text:
#   Unit "cca.jl" of the Diagonalization.jl Package for Julia language
#
#   MIT License
#   Copyright (c) 2019-2021,
#   Marco Congedo, CNRS, Grenoble, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements Maximum Covariance Analysys (MCA) and
#   Canonical Correlation Analysis(CCA) filters, that is,
#   filters based on the singular value decomposition.

"""
```julia
(1)
function mca(Cxy :: Mat;
             eVar      :: TeVaro = ○,
             eVarMeth  :: Function = searchsortedfirst,
             simple    :: Bool = false)

(2)
function mca(X::Mat, Y::Mat;
             dims     :: Into = ○,
             meanX    :: Tmean = 0,
             meanY    :: Tmean = 0,
             wXY      :: Tw = ○,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false)

(3)
function mca(𝐗::VecMat, 𝐘::VecMat;
             dims     :: Into = ○,
             meanX    :: Into = 0,
             meanY    :: Into = 0,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false)
```

Return a [LinearFilter](@ref) object:

**(1) Maximum covariance analysis** with real or complex covariance matrix
`Cxy` of dimension ``n_x⋅n_y`` as input.

Differently from [PCA](@ref) and [Whitening](@ref), `Cxy`
is a generic `Matrix` object since it is not symmetric/Hermitian,
left alone square.

`eVar` and `eVarMeth` are keyword optional arguments for defining the
[subspace dimension](@ref) ``p`` using the `.arev` vector given by Eq. [mca.2].

The default values are:
- `eVar=0.999`
- `eVarMeth=searchsortedfirst`

If `simple` is set to `true`, ``p`` is set equal to ``min(n_x, n_y)``
and only the fields `.F` and `.iF` are written in the constructed object.
This option is provided for low-level work when you don't need to define
a subspace dimension or you want to define it by your own methods.

**(2) Maximum covariance analysis**
with real or complex data matrices `X` and `Y` as input.

`dims`, `meanX`, `meanY` and `wXY` are optional keyword arguments to
regulate the estimation of the cross-covariance matrix ``C_{xy}``.
Particularly (see [covariance matrix estimations](@ref)),
- `meanX` is the `mean` argument for data matrix `X`.
- `meanY` is the `mean` argument for data matrix `Y`.
- `wXY` is the weight argument for estimating a weighted cross-covariance matrix
  ``C_{XY}``.

Once the cross-covariance matrix estimated, method (1) is invoked
with optional keyword arguments `eVar`, `eVarMeth` and `simple`.
See method (1) for details.

**(3) Maximum covariance analysis**
with two vectors of real or complex data matrices
`𝐗` and `𝐘` as input. `𝐗` and `𝐘` must hold the same number of matrices
and corresponding pairs of matrices therein must comprise the same
number of samples.

`dims`, `meanX` and `meanY` are optional keyword arguments to
regulate the estimation of the cross-covariance matrices for
all pairs of corresponding data matrices in `𝐗` and `𝐘`.
See method (2) and [covariance matrix estimations](@ref).

The arithmetic mean of these cross-covariance matrices
is computed and method (1) is invoked with optional
keyword arguments `eVar`, `eVarMeth` and `simple`.
See method (1) for details.

**See also:** [PCA](@ref), [CCA](@ref), [gMCA](@ref), [mAJD](@ref).

**Examples:**

```julia
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
n, t=10, 100
X=genDataMatrix(n, t)
Y=genDataMatrix(n, t)
Cx=Symmetric((X*X')/t)
Cy=Symmetric((Y*Y')/t)
Cxy=(X*Y')/t
mC=mca(Cxy, simple=true)
@test Cxy≈mC.F[1]*mC.D*mC.F[2]'
D=mC.F[1]'Cxy*mC.F[2]
@test norm(D-Diagonal(D))+1. ≈ 1.

# Method (1) complex
Xc=genDataMatrix(ComplexF64, n, t)
Yc=genDataMatrix(ComplexF64, n, t)
Cxc=Symmetric((Xc*Xc')/t)
Cyc=Symmetric((Yc*Yc')/t)
Cxyc=(Xc*Yc')/t
mCc=mca(Cxyc, simple=true)
@test Cxyc≈mCc.F[1]*mCc.D*mCc.F[2]'
Dc=mCc.F[1]'Cxyc*mCc.F[2]
@test norm(Dc-Diagonal(Dc))+1. ≈ 1.

# Method (2) real
mXY=mca(X, Y, simple=true)
D=mXY.F[1]'*Cxy*mXY.F[2]
@test norm(D-Diagonal(D))+1≈1.
@test mXY==mC

# Method (2) complex
mXYc=mca(Xc, Yc, simple=true)
Dc=mXYc.F[1]'*Cxyc*mXYc.F[2]
@test norm(Dc-Diagonal(Dc))+1. ≈ 1.
@test mXYc==mCc

# Method (3) real
# maximum covariance analysis of the average covariance and cross-covariance
k=10
Xset=[genDataMatrix(n, t) for i=1:k]
Yset=[genDataMatrix(n, t) for i=1:k]

m=mca(Xset, Yset)

# ... selecting subspace dimension allowing an explained variance = 0.5
m=mca(Xset, Yset; eVar=0.5)

# ... subtracting the mean from the matrices in Xset and Yset
m=mca(Xset, Yset; meanX=nothing, meanY=nothing, eVar=0.5)

# mca on the average of the covariance and cross-covariance matrices
# computed along dims 1
m=mca(Xset, Yset; dims=1, eVar=0.5)

# name of the filter
m.name

using Plots
# plot regularized accumulated eigenvalues
plot(m.arev)

# plot the original cross-covariance matrix and the rotated
# cross-covariance matrix
 Cmax=maximum(abs.(Cxy));
 h1 = heatmap(Cxy, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="Cxy");
 D=mC.F[1]'*Cxy*mC.F[2];
 Dmax=maximum(abs.(D));
 h2 = heatmap(D, clim=(0, Dmax), yflip=true, c=:amp, title="F[1]'*Cxy*F[2]");
 📈=plot(h1, h2, size=(700,300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigMCA.png")
```

 ![Figure MCA](assets/FigMCA.png)

```julia
# Method (3) complex
# maximum covariance analysis of the average covariance and cross-covariance

k=10
Xcset=[genDataMatrix(ComplexF64, n, t) for i=1:k]
Ycset=[genDataMatrix(ComplexF64, n, t) for i=1:k]

mc=mca(Xcset, Ycset)

# ... selecting subspace dimension allowing an explained variance = 0.5
mc=mca(Xcset, Ycset; eVar=0.5)

# ... subtracting the mean from the matrices in Xset and Yset
mc=mca(Xcset, Ycset; meanX=nothing, meanY=nothing, eVar=0.5)

# mca on the average of the covariance and cross-covariance matrices
# computed along dims 1
mc=mca(Xcset, Ycset; dims=1, eVar=0.5)

# name of the filter
mc.name
```

"""
function mca(Cxy :: Mat;
             eVar      :: TeVaro = ○,
             eVarMeth  :: Function = searchsortedfirst,
             simple    :: Bool = false)

  args=("Maximum Covariance Anaysis", false)
  M=Matrix
  U, λ, V = svd(Cxy; full=true)
  λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

  # typecast V, U', V' U1' and U2' as `Matrix` since they are of the Adjoint type
  simple ? LF([U, M(V)], [M(U'), M(V')], Diagonal(λ), ○, ○, ○, args...) :
  begin
     eVar, D, U1, U2, p, arev=_ssdxy!(eVar, λ, U, M(V), _minDim(Cxy), eVarMeth)
     LF([U1, U2], [M(U1'), M(U2')], D, eVar, λ, arev, args...)
  end
end

# MCA with data as input (separated data): API
function mca(X::Mat, Y::Mat;
             dims     :: Into = ○,
             meanX    :: Tmean = 0,
             meanY    :: Tmean = 0,
             wXY      :: Tw = ○,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false)

   dims===○ && (dims=_set_dims(X, Y))
   _check_data(X, Y, dims, SCM, meanX, meanY, wXY)===○ && return

   return mca(_cov(X, Y, dims, meanX, meanY, wXY);
              eVar=eVar, eVarMeth=eVarMeth, simple=simple)
end

function mca(𝐗::VecMat, 𝐘::VecMat;
             dims     :: Into = ○,
             meanX    :: Into = 0,
             meanY    :: Into = 0,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false)

   dims===○ && (dims=_set_dims(𝐗, 𝐘))
   _check_data(𝐗, 𝐘, dims, SCM, meanX, meanY, ○)===○ && return

   𝐂=_cov(𝐗, 𝐘; dims=dims, meanX=meanX, meanY=meanY)

   mca(PosDefManifold.fVec(mean, Vector{Matrix}(𝐂)),
       eVar=eVar, eVarMeth=eVarMeth, simple=simple)
end

"""
```julia
(1)
function cca(Cx :: SorH, Cy :: SorH, Cxy :: Mat;
             eVarCx   :: TeVaro=○,
             eVarCy   :: TeVaro=○,
             eVar     :: TeVaro=○,
             eVarMeth :: Function=searchsortedfirst,
             simple   :: Bool=false)

(2)
function cca(X::Mat, Y::Mat;
             covEst   :: StatsBase.CovarianceEstimator=SCM,
             dims     :: Into = ○,
             meanX    :: Tmean = 0,
             meanY    :: Tmean = 0,
             wX       :: Tw = ○,
             wY       :: Tw = ○,
             wXY      :: Tw = ○,
          eVarCx   :: TeVaro = ○,
          eVarCy   :: TeVaro = ○,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false)

(3)
function cca(𝐗::VecMat, 𝐘::VecMat;
             covEst   :: StatsBase.CovarianceEstimator=SCM,
             dims     :: Into = ○,
             meanX    :: Into = 0,
             meanY    :: Into = 0,
          eVarCx   :: TeVaro = ○,
          eVarCy   :: TeVaro = ○,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false,
       metric   :: Metric = Euclidean,
       wCx      :: Vector = [],
       wCy      :: Vector = [],
       ✓w       :: Bool = true,
       initCx   :: SorHo = nothing,
       initCy   :: SorHo = nothing,
       tol      :: Real = 0.,
       verbose  :: Bool = false)
```

Return a [LinearFilter](@ref) object:

**(1) Canonical correlation analysis**
with, as input, real or complex:
- covariance matrix `Cx` of dimension ``n_x⋅n_x``,
- covariance matrix `Cy` of dimension ``n_y⋅n_y`` and
- cross-covariance matrix `Cxy` of dimension ``n_x⋅n_y``.

`Cx` and `Cy` must be flagged as `Symmetric`, if real or `Hermitian`,
if real or complex, see [data input](@ref). Instead `Cxy`
is a generic `Matrix` object since it is not symmetric/Hermitian,
left alone square.

`eVarCx`, `eVarCy`, `eVar` and `eVarMeth` are keyword optional arguments
for defining the [subspace dimension](@ref) ``p``. Particularly:
- `eVarCx` and `eVarCy` are used to determine a subspace dimension
  in the whitening of `Cx` and `Cy`, i.e., in the first step of the two-step
  procedure descrived above in the **solution** section.
- `eVar` is used to determine the final subspace dimension using the `.arev`
  vector given by Eq. [cca.5] here above.
- `eVarMeth` applies to `eVarCx`, `eVarCy` and `eVar`.

The default values are:
- `eVarCx=0.999`
- `eVarCy=0.999`
- `eVar=0.999`
- `eVarMeth=searchsortedfirst`

If `simple` is set to `true`, ``p`` is set equal to ``min(n_x, n_y)``
and only the fields `.F` and `.iF` are written in the constructed object.
This option is provided for low-level work when you don't need to define
a subspace dimension or you want to define it by your own methods.

**(2) Canonical correlation analysis**
with real or complex data matrices `X` and `Y` as input.

`covEst`, `dims`, `meanX`, `meanY`,  `wX`, `wY` and `wXY` are optional
keyword arguments to regulate the estimation of the
covariance matrices of `X` and `Y` and their cross-covariance matrix ``C_{xy}``.
Particularly (See [covariance matrix estimations](@ref)),
- `meanX` is the `mean` argument for data matrix `X`.
- `meanY` is the `mean` argument for data matrix `Y`.
- `wX` is the `w` argument for estimating a weighted covariance matrix ``C_x``.
- `wY` is the `w` argument for estimating a weighted covariance matrix ``C_y``.
- `wXY` is the `w` argument for estimating a weighted
   cross-covariance matrix ``C_{XY}``.
- `covEst` applies only to the estimation of ``C_x`` and ``C_y``.

Once ``C_x``, ``C_y`` and ``C_{xy}`` estimated, method (1) is invoked
with optional keyword arguments `eVar`, `eVarCx`, `eVarCy`,
`eVarMeth` and `simple`. See method (1) for details.

**(3) Canonical correlation analysis**
with two vectors of real or complex data matrices
`𝐗` and `𝐘` as input. `𝐗` and `𝐘` must hold the same number of matrices
and corresponding pairs of matrices therein must comprise the same
number of samples.

`covEst`, `dims`, `meanX` and `meanY` are optional keyword arguments to
regulate the estimation of the covariance matrices of all matrices in
`𝐗` and `𝐘` and the cross-covariance matrices for
all pairs of their corresponding data matrices.
See method (2) and [covariance matrix estimations](@ref).

A mean of these covariance matrices is computed.
For the cross-covariance matrices the arithmetic mean is used.
For the covariance matrices of `𝐗` and `𝐘`,
optional keywords arguments `metric`, `wCx`, `wCy`, `✓w`,
`initCx`, `initCy`, `tol` and `verbose` are used to allow
non-Euclidean mean estimations. Particularly
(see [mean covariance matrix estimations](@ref)),
- `wCx` are the weights for the covariance matrices of `𝐗`,
- `wCy` are the weights for the covariance matrices of `𝐘`,
- `initCx` is the initialization for the mean of the covariance matrices of `𝐗`,
- `initCy` is the initialization for the mean of the covariance matrices of `𝐘`.
By default, the arithmetic mean is computed.

Once the mean covariance and cross-covariance matrices are estimated,
method (1) is invoked with optional
keyword arguments `eVarCx`, `eVarCy`, `eVar`, `eVarMeth` and `simple`.
See method (1) for details.

**See also:** [Whitening](@ref), [MCA](@ref), [gCCA](@ref), [mAJD](@ref).

**Examples:**

```julia
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
n, t=10, 100
X=genDataMatrix(n, t)
Y=genDataMatrix(n, t)
Cx=Symmetric((X*X')/t)
Cy=Symmetric((Y*Y')/t)
Cxy=(X*Y')/t
cC=cca(Cx, Cy, Cxy, simple=true)
@test cC.F[1]'*Cx*cC.F[1]≈I
@test cC.F[2]'*Cy*cC.F[2]≈I
D=cC.F[1]'*Cxy*cC.F[2]
@test norm(D-Diagonal(D))+1. ≈ 1.

# Method (1) complex
Xc=genDataMatrix(ComplexF64, n, t)
Yc=genDataMatrix(ComplexF64, n, t)
Cxc=Hermitian((Xc*Xc')/t)
Cyc=Hermitian((Yc*Yc')/t)
Cxyc=(Xc*Yc')/t
cCc=cca(Cxc, Cyc, Cxyc, simple=true)
@test cCc.F[1]'*Cxc*cCc.F[1]≈I
@test cCc.F[2]'*Cyc*cCc.F[2]≈I
Dc=cCc.F[1]'*Cxyc*cCc.F[2]
@test norm(Dc-Diagonal(Dc))+1. ≈ 1.

# Method (2) real
cXY=cca(X, Y, simple=true)
@test cXY.F[1]'*Cx*cXY.F[1]≈I
@test cXY.F[2]'*Cy*cXY.F[2]≈I
D=cXY.F[1]'*Cxy*cXY.F[2]
@test norm(D-Diagonal(D))+1. ≈ 1.
@test cXY==cC

# Method (2) complex
cXYc=cca(Xc, Yc, simple=true)
@test cXYc.F[1]'*Cxc*cXYc.F[1]≈I
@test cXYc.F[2]'*Cyc*cXYc.F[2]≈I
Dc=cXYc.F[1]'*Cxyc*cXYc.F[2]
@test norm(Dc-Diagonal(Dc))+1. ≈ 1.
@test cXYc==cCc

# Method (3) real
# canonical correlation analysis of the average covariance and cross-covariance
k=10
Xset=[genDataMatrix(n, t) for i=1:k]
Yset=[genDataMatrix(n, t) for i=1:k]

c=cca(Xset, Yset)

# ... selecting subspace dimension allowing an explained variance = 0.9
c=cca(Xset, Yset; eVar=0.9)

# ... subtracting the mean from the matrices in Xset and Yset
c=cca(Xset, Yset; meanX=nothing, meanY=nothing, eVar=0.9)

# cca on the average of the covariance and cross-covariance matrices
# computed along dims 1
c=cca(Xset, Yset; dims=1, eVar=0.9)

# name of the filter
c.name

using Plots
# plot regularized accumulated eigenvalues
plot(c.arev)

# plot the original covariance and cross-covariance matrices
# and their transformed counterpart
 CxyMax=maximum(abs.(Cxy));
 h1 = heatmap(Cxy, clim=(-CxyMax, CxyMax), title="Cxy", yflip=true, c=:bluesreds);
 D=cC.F[1]'*Cxy*cC.F[2];
 Dmax=maximum(abs.(D));
 h2 = heatmap(D, clim=(0, Dmax), title="F1'CxyF2", yflip=true, c=:amp);
 CxMax=maximum(abs.(Cx));
 h3 = heatmap(Cx, clim=(-CxMax, CxMax), title="Cx", yflip=true, c=:bluesreds);
 h4 = heatmap(cC.F[1]'*Cx*cC.F[1], clim=(0, 1), title="F1'CxF1", yflip=true, c=:amp);
 CyMax=maximum(abs.(Cy));
 h5 = heatmap(Cy, clim=(-CyMax, CyMax), title="Cy", yflip=true, c=:bluesreds);
 h6 = heatmap(cC.F[2]'*Cy*cC.F[2], clim=(0, 1), title="F2'CyF2", yflip=true, c=:amp);
 📈=plot(h3, h5, h1, h4, h6, h2, size=(800,400))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigCCA.png")

```

![Figure CCA](assets/FigCCA.png)

```julia
# Method (3) complex
# canonical correlation analysis of the average covariance and cross-covariance
k=10
Xcset=[genDataMatrix(ComplexF64, n, t) for i=1:k]
Ycset=[genDataMatrix(ComplexF64, n, t) for i=1:k]

cc=cca(Xcset, Ycset)

# ... selecting subspace dimension allowing an explained variance = 0.9
cc=cca(Xcset, Ycset; eVar=0.9)

# ... subtracting the mean from the matrices in Xset and Yset
cc=cca(Xcset, Ycset; meanX=nothing, meanY=nothing, eVar=0.9)

# cca on the average of the covariance and cross-covariance matrices
# computed along dims 1
cc=cca(Xcset, Ycset; dims=1, eVar=0.9)

# name of the filter
cc.name
```

"""
function cca(Cx :: SorH, Cy :: SorH, Cxy :: Mat;
             eVarCx   :: TeVaro=○,
             eVarCy   :: TeVaro=○,
             eVar     :: TeVaro=○,
             eVarMeth :: Function=searchsortedfirst,
             simple   :: Bool=false)

   size(Cx, 1)==size(Cxy, 1) || throw(ArgumentError(📌*", cca function: Matrices `Cx` and `Cxy` must have the same number of columns"))
   size(Cy, 2)==size(Cxy, 2) || throw(ArgumentError(📌*", cca function: Matrices `Cy` and `Cxy` must have the same number of rows"))
   if eVar isa Int
      eVarCx isa Int && eVar>eVarCx && throw(ArgumentError(📌*", cca function: `eVar` cannot be larger than `eVarCx`"))
      eVarCy isa Int && eVar>eVarCy && throw(ArgumentError(📌*", cca function: `eVar` cannot be larger than `eVarCy`"))
   end

   args=("Canonical Correlation Analysis", false)

   x=whitening(Cx; eVar=eVarCx, eVarMeth=eVarMeth)
   y=whitening(Cy; eVar=eVarCy, eVarMeth=eVarMeth)

   # alert user if eVar passed as an integer exceeds the minimum whitening dim
   whiteDim=min(size(x.F, 2), size(y.F, 2))
   if eVar isa Int && eVar>whiteDim
      @warn(📌*", cca function: the whitening step reduced the rank to $(whiteDim); `eVar` has been lowered to this value.")
      eVar=whiteDim
   end

   m=mca(x.F'*Cxy*y.F; eVar=eVar, eVarMeth=eVarMeth)

   if simple
      LF([x.F*m.F[1], y.F*m.F[2]], [m.iF[1]*x.iF, m.iF[2]*y.iF], m.D, ○, ○, ○, args...)
   else
      e, D, U, V, p, arev=_ssdxy!(eVar, diag(m.D), m.F[1], m.F[2], _minDim(Cxy), eVarMeth)
      LF([x.F*U, y.F*V], [U'*x.iF, V'*y.iF], D, e, diag(m.D), arev, args...)
   end
end

function cca(X::Mat, Y::Mat;
             covEst   :: StatsBase.CovarianceEstimator=SCM,
             dims     :: Into = ○,
             meanX    :: Tmean = 0,
             meanY    :: Tmean = 0,
             wX       :: Tw = ○,
             wY       :: Tw = ○,
             wXY      :: Tw = ○,
          eVarCx   :: TeVaro = ○,
          eVarCy   :: TeVaro = ○,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false)

   dims===○ && (dims=_set_dims(X, Y))
   _check_data(X, Y, dims, covEst, meanX, meanY, wXY)===○ && return

   cca(_cov(X, covEst, dims, meanX, wX),
       _cov(Y, covEst, dims, meanY, wY),
       _cov(X, Y, dims, 0, 0, wXY);
       eVarCx=eVarCx, eVarCy=eVarCy, eVar=eVar, eVarMeth=eVarMeth,
       simple=simple)
end

function cca(𝐗::VecMat, 𝐘::VecMat;
             covEst   :: StatsBase.CovarianceEstimator=SCM,
             dims     :: Into = ○,
             meanX    :: Into = 0,
             meanY    :: Into = 0,
          eVarCx   :: TeVaro = ○,
          eVarCy   :: TeVaro = ○,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false,
       metric   :: Metric = Euclidean,
       wCx      :: Vector = [],
       wCy      :: Vector = [],
       ✓w       :: Bool = true,
       initCx   :: SorHo = nothing,
       initCy   :: SorHo = nothing,
       tol      :: Real = 0.,
       verbose  :: Bool = false)

   Metric==VonNeumann && throw(ArgumentError(📌*", cca function: A solution for the mean is not available for the Von Neumann metric. Use another metric as `metric` argument"))
   dims===○ && (dims=_set_dims(𝐗, 𝐘))
   _check_data(𝐗, 𝐘, dims, covEst, meanX, meanY, ○)===○ && return

   𝐂x= _cov(𝐗; covEst=covEst, dims = dims, meanX = meanX)
   𝐂y= _cov(𝐘; covEst=covEst, dims = dims, meanX = meanX)
   𝐂xy=_cov(𝐗, 𝐘; dims=dims, meanX=0, meanY=0)

   cca(mean(metric, 𝐂x;
            w = wCx, ✓w = ✓w,
            init = initCx===○ ? ○ : Hermitian(initCx), #just init here when you upfate PodDefManifold
            tol = tol, verbose = verbose),
       mean(metric, 𝐂y;
            w = wCy, ✓w = ✓w,
            init = initCy===○ ? ○ : Hermitian(initCy), #just init here when you upfate PodDefManifold
            tol = tol, verbose = verbose),
       PosDefManifold.fVec(mean, Vector{Matrix}(𝐂xy)), # multi-threaded Euclidean mean
       eVarCx=eVarCx, eVarCy=eVarCy, eVar=eVar, eVarMeth=eVarMeth,
       simple=simple)
end


text:
# Truncated Inverse Gaussian distribution

TruncatedInverseGaussian(μ::Float64, λ::Float64, a::Float64, b::Float64) =
    Distributions.Truncated(Distributions.InverseGaussian(μ, λ), a, b)
TruncatedInverseGaussian(μ::Real, λ::Real, a::Real, b::Real) =
    TruncatedInverseGaussian(Float64(μ), Float64(λ), Float64(a), Float64(b))

# ### statistics

function Distributions.minimum(
    d::Distributions.Truncated{Distributions.InverseGaussian{T},Distributions.Continuous},
) where {T<:Real}
    min(d.lower, Distributions.minimum(d.untruncated))
end
function Distributions.maximum(
    d::Distributions.Truncated{Distributions.InverseGaussian{T},Distributions.Continuous},
) where {T<:Real}
    max(d.lower, Distributions.maximum(d.untruncated))
end

## sampling

### sampler for upper truncated inverse gaussian

function Distributions.rand(
    rng::Distributions.AbstractRNG,
    d::Distributions.Truncated{Distributions.InverseGaussian{T},Distributions.Continuous},
) where {T<:Real}
    mu = d.untruncated.μ
    lambda = d.untruncated.λ
    R = d.upper
    if d.lower == 0.0
        #upper Truncated
        X = R + 1
        if (mu > R)
            alpha = 0.0
            while (Random.rand(rng) > alpha)
                ## X = R + 1
                ## while (X > R) {
                ##     X = 1.0 / rgamma(1, 0.5, rate=0.5)
                ## }
                E = [Random.randexp(rng) for i = 1:2]
                while (E[1]^2 > 2 * E[2] / R)
                    E = [Random.randexp(rng) for i = 1:2]
                end
                X = R / (1 + R * E[1])^2
                alpha = exp(-0.5 * (1 / mu)^2 * X)
            end
        else
            while (X > R)
                lambda = 1.0
                Y = Random.randn(rng)^2
                X =
                    mu + 0.5 * mu^2 / lambda * Y -
                    0.5 * mu / lambda * sqrt(4 * mu * lambda * Y + (mu * Y)^2)
                if (Random.randexp(rng) > mu / (mu + X))
                    X = mu^2 / X
                end
            end
        end
        return X
    else
        error("The sampler is available only for upper TruncatedInverseGaussian.")
    end
end

ext: jl


text:
using Test, NaturallyUnitful

@testset "tests" begin
    @test natural(1u"m") == 5.067730759202785e6u"eV^-1"

    @test unnatural(u"m", 5.067730759202785e6u"eV^-1") == 1.0u"m"

    @test natural(1e8u"m/s") == 0.33356409519815206

    @test unnatural(u"m/s", 0.33356409519815206) == 1e8u"m/s"

    @test natural(1u"kg", base=u"GeV") ==  5.609588650020685e26u"GeV"

    @test natural(1u"m", base=u"GeV") ==  5.067730759202785e15u"GeV^-1"

    @test unnatural(u"K", 1u"eV") == 11604.522060401006u"K"
end


text:
#= NOTE: If running tests incrementally, import the following packages:
using Revise, Test
=#

import Cytof5.Util: similarity_FM, meanabsdiff, sumabsdiff

@testset "Util similarity_FM" begin
  I = 3
  J = 20
  K1 = 5
  K2 = 4

  W1 = rand(I, K1)
  W1 ./= sum(W1, dims=2)

  Z1 = Int.(rand(J, K1) .> .5)
  Z2 = Int.(rand(J, K2) .> .5)

  @time s = similarity_FM(Z1, W1, Z2)
  @assert s <= J * I

  @test true
end


text:
@compile reference_vadd """
__global__ void reference_vadd(const float *a, const float *b, float *c)
{
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    c[i] = a[i] + b[i];
}
"""

@target ptx function kernel_vadd(a::CuDeviceArray{Float32}, b::CuDeviceArray{Float32},
                                 c::CuDeviceArray{Float32})
    i = blockId_x() + (threadId_x()-1) * numBlocks_x()
    c[i] = a[i] + b[i]

    return nothing
end


text:
# DARNN.jl
#
# Layers for a DA-RNN network
# Based on paper Quin et. al. "A Dual-Stage Attention-Based Recurrent Neural
# Network for Time Series Prediction" https://arxiv.org/abs/1704.02971

# struct that stores layers and some extra information
mutable struct DARNNCell{A, B, C, D, E, F, W, X, Y, Z}
    # Encoder part
	encoder_lstm::A
	encoder_attn::B
    # Decoder part
	decoder_lstm::C
	decoder_attn::D
	decoder_fc::E
	decoder_fc_final::F
    # Index for original data etc
	encodersize::W
	decodersize::X
	orig_idx::Y
	poollength::Z
end

# Initialization function to obtain the right size of the hidden state
# after a `reset!` of the model.
function darnn_init(m::DARNNCell,x)
	m.encoder_lstm(x[:,1,1,:])
	m.decoder_lstm(x[m.orig_idx,1,1,:]')
	return nothing
end
Flux.Zygote.@nograd darnn_init

# Dispatch on `trainable` and `reset!` to make standard Flux scripts work
Flux.trainable(m::DARNNCell) = (m.encoder_lstm, m.encoder_attn, m.decoder_lstm,
    m.decoder_attn, m.decoder_fc, m.decoder_fc_final)
Flux.reset!(m::DARNNCell) = Flux.reset!.((m.encoder_lstm, m.decoder_lstm))

"""
    DARNN(inp, encodersize, decodersize, poollength, orig_idx)

Creates a DA-RNN layer based on the architecture described in
[Qin et. al.](https://arxiv.org/abs/1704.02971), as implemented for PyTorch
[here](https://github.com/Seanny123/da-rnn/blob/master/modules.py). `inp` specifies the
number of input features. `encodersize` defines the number of LSTM encoder layers, and
`decodersize` defines the number of LSTM decoder layers. `poolsize` gives the
length of the window for the pooled input data, and `orig_idx` defines the array
index where the original time series is stored in the input data,

Data is expected as array with dimensions `features x poolsize x 1 x data`, i.e.
for 1000 data points containing 31 features that have been windowed over 6
timesteps, `DARNN` expects an input size of `(31, 6, 1, 1000)`.

Takes the keyword arguments `initW` and `initb` for the initialization of the
weight vector and bias of the linear layers.
"""
function DARNN(inp::Integer, encodersize::Integer, decodersize::Integer, poollength::Integer, orig_idx::Integer;
	initW = Flux.glorot_uniform, initb = Flux.zeros)

	# Encoder part
	encoder_lstm = LSTM(inp, encodersize)
	encoder_attn = Chain(Dense(2*encodersize + poollength, poollength, initW=initW, initb=initb),
	                    a -> tanh.(a),
	                    Dense(poollength,1,initW=initW, initb=initb))
	# Decoder part
	decoder_lstm = LSTM(1, decodersize)
	decoder_attn = Chain(Dense(2*decodersize + encodersize, encodersize, initW=initW, initb=initb),
	                    a -> tanh.(a),
	                    Dense(encodersize, 1, initW=initW, initb=initb))
	decoder_fc = Dense(encodersize + 1, 1)
	decoder_fc_final = Dense(decodersize + encodersize, 1, initW=initW, initb=initb)

	return DARNNCell(encoder_lstm, encoder_attn, decoder_lstm, decoder_attn, decoder_fc,
	 		  decoder_fc_final, encodersize, decodersize, orig_idx, poollength)
end

# model output
function (m::DARNNCell)(x)
	# Initialize after reset
	size(m.encoder_lstm.state[1],2) == 1 && darnn_init(m,x)
	input_data = dropdims(x; dims=3)
	input_encoded = darnn_encoder(m, input_data)
	context = darnn_decoder(m, input_encoded, input_data)
	return m.decoder_fc_final( cat(m.decoder_lstm.state[1], context, dims=1))
end

function darnn_encoder(m::DARNNCell, input_data)
	input_encoded = Flux.Zygote.Buffer(similar(input_data, m.encodersize, m.poollength,
	 				size(input_data,3)))
	@inbounds for t in 1:m.poollength
	      hidden = m.encoder_lstm.state[1]
	      cell = m.encoder_lstm.state[2]

	      x = cat(repeat(hidden, inner=(1,1,size(input_data,1))),
	              repeat(cell, inner=(1,1,size(input_data,1))),
	              permutedims(input_data,[2,3,1]), dims=1) |>
	           a -> m.encoder_attn(reshape(a, (:, size(input_data,1)*size(input_data,3) ) ))
	      attn_weights = Flux.softmax( reshape(x, (size(input_data,1), size(input_data,3))  ) )
	      weighted_input = attn_weights .* input_data[:,t,:]
	      _ = m.encoder_lstm(weighted_input)

	      input_encoded[:,t,:] = Flux.unsqueeze(m.encoder_lstm.state[1],2)
	end
	return copy(input_encoded)
end

function darnn_decoder(m::DARNNCell, input_encoded, input_data)
	context = similar(input_data, m.encodersize, size(input_data,3))
	@inbounds for t in 1:m.poollength
	      hidden = m.decoder_lstm.state[1]
	      cell = m.decoder_lstm.state[2]
	      x = cat(permutedims(repeat(hidden, inner=(1,1,m.poollength)), [1,3,2]),
	              permutedims(repeat(cell, inner=(1,1,m.poollength)), [1,3,2]),
	              input_encoded, dims=1) |>
	          a -> m.decoder_attn(reshape(a, (2*m.decodersize + m.encodersize,:))) |>
	          a -> Flux.softmax(reshape(a, (m.poollength,:)))
	      context = dropdims(NNlib.batched_mul(input_encoded, Flux.unsqueeze(x,2) ), dims=2)
	      ỹ = m.decoder_fc(cat(context, input_data[m.orig_idx,t,:]', dims=1))
	      _ = m.decoder_lstm(ỹ)
	end
	return context
end

# pretty printing
function Base.show(io::IO, l::DARNNCell)
	print(io, "DARNN(", size(l.encoder_lstm.cell.Wi,2))
	print(io, ", ", l.encodersize)
	print(io, ", ", Int(size(l.decoder_lstm.cell.Wi,1)/4))
	print(io, ", ", l.poollength)
	print(io, ", ", l.orig_idx)
	print(io, ")")
end

max_stars_repo_name: AndyTGardner/FluxArchitectures

max_issues_repo_name: AndyTGardner/FluxArchitectures

max_forks_repo_name: AndyTGardner/FluxArchitectures


text:
abstract type ReservoirGeometry end

surface_area(g::ReservoirGeometry, z) =
    error("Please implement a surface_area function for type ", typeof(g))

struct Extrusion <: ReservoirGeometry
    A::Float64
end

surface_area(g::Extrusion, z) = g.A
volume(g::Extrusion, z) = g.A * z

struct DiscreteCurve <: ReservoirGeometry
    n::Int
    z::Vector{Float64}
    S::Vector{Float64}
end

surface_area(g::DiscreteCurve, z) = _linear_interpolation(g.n, g.z, g.S, z)

function volume(g::DiscreteCurve, z)
    V = 0.0
    for i = 2:g.n
        V += (S[i - 1] + S[i]) * (z[i] - z[i - 1])
    end
    return V / 2.0
end


text:
using Base.Test
using Lora

fnames = (
  :value,
  :loglikelihood,
  :logprior,
  :logtarget,
  :gradloglikelihood,
  :gradlogprior,
  :gradlogtarget,
  :tensorloglikelihood,
  :tensorlogprior,
  :tensorlogtarget,
  :dtensorloglikelihood,
  :dtensorlogprior,
  :dtensorlogtarget,
  :diagnosticvalues
)
filepath = dirname(@__FILE__)
filesuffix = "csv"
filenames = Array(AbstractString, 14)
for i in 1:14
  filenames[i] = joinpath(filepath, string(fnames[i])*"."*filesuffix)
end

println("    Testing BasicContParamIOStream constructors and close method...")

iostreamsize = ()
iostreamn = 4
iostream = BasicContParamIOStream(iostreamsize, iostreamn, filepath=filepath)

@test isa(iostream.value, IOStream)
for i in 2:14
  @test iostream.(fnames[i]) == nothing
end
@test iostream.size == iostreamsize
@test iostream.n == iostreamn
@test length(iostream.diagnostickeys) == 0

close(iostream)
rm(filenames[1])

iostreamsize = (2,)
iostreamn = 4
iostream = BasicContParamIOStream(
  iostreamsize,
  iostreamn,
  monitor=[true; fill(false, 5); true; fill(false, 6)],
  filepath=filepath,
  diagnostickeys=[:accept]
)

for i in (1, 7, 14)
  @test isa(iostream.(fnames[i]), IOStream)
end
for i in [2:6; 8:13]
  @test iostream.(fnames[i]) == nothing
end
@test iostream.size == iostreamsize
@test iostream.n == iostreamn
@test length(iostream.diagnostickeys) == 1

close(iostream)
for i in (1, 7, 14); rm(filenames[i]); end

iostreamsize = (3,)
iostreamn = 7
iostream = BasicContParamIOStream(iostreamsize, iostreamn, [:value, :logtarget], filepath=filepath, diagnostickeys=[:accept])

for i in (1, 4, 14)
  @test isa(iostream.(fnames[i]), IOStream)
end
for i in [2:3; 5:13]
  @test iostream.(fnames[i]) == nothing
end
@test iostream.size == iostreamsize
@test iostream.n == iostreamn
@test length(iostream.diagnostickeys) == 1

close(iostream)
for i in [1, 4, 14]; rm(filenames[i]); end

println("    Testing BasicContParamIOStream IO methods...")

println("      Interaction with BasicContUnvParameterState...")

nstatev = Float64[5.70, 1.44, -1.21, 5.67]
iostreamsize = ()
iostreamn = length(nstatev)

iostream = BasicContParamIOStream(iostreamsize, iostreamn, filepath=filepath)
for v in nstatev
  write(iostream, BasicContUnvParameterState(v))
end

close(iostream)

iostream = BasicContParamIOStream(iostreamsize, iostreamn, filepath=filepath, mode="r")
nstate = read(iostream, Float64)

@test isa(nstate, ContUnvMarkovChain{Float64})
@test nstate.value == nstatev
for i in 2:13
  @test length(nstate.(fnames[i])) == 0
end
@test length(nstate.diagnostickeys) == 0
@test size(nstate.diagnosticvalues) == (0, 0)
@test nstate.n == iostream.n

close(iostream)
rm(filenames[1])

println("      Interaction with ContUnvMarkovChain...")

nstatev = Float32[1.93, 98.46, -3.61, -0.99, 74.52, 9.90]
nstated = Any[false, true, true, false, true, false]'
iostreamsize = ()
iostreamn = length(nstatev)

iostream = BasicContParamIOStream(iostreamsize, iostreamn, [:value], filepath=filepath, diagnostickeys=[:accept])
nstatein = ContUnvMarkovChain(iostreamn, [:value], [:accept], Float32)
nstatein.value = nstatev
nstatein.diagnosticvalues = nstated
write(iostream, nstatein)

close(iostream)

iostream = BasicContParamIOStream(iostreamsize, iostreamn, [:value], filepath=filepath, diagnostickeys=[:accept], mode="r")
nstateout = read(iostream, Float32)

@test isa(nstateout, ContUnvMarkovChain{Float32})
@test nstateout.value == nstatein.value
for i in 2:13
  @test length(nstateout.(fnames[i])) == 0
end
@test length(nstateout.diagnostickeys) == 1
@test nstateout.diagnosticvalues == nstatein.diagnosticvalues
@test nstateout.n == nstatein.n

close(iostream)
for i in (1, 14); rm(filenames[i]); end

println("      Interaction with BasicContMuvParameterState...")

nstatev = Float64[1.33 2.44 3.14 -0.82; 7.21 -9.75 -5.26 -0.63]
nstategll = Float64[3.13 -12.10 13.11 -0.99; 9.91 -5.25 -8.15 -9.69]
nstated = Any[false, true, true, false]'
iostreamsize = (size(nstatev, 1),)
iostreamn = size(nstatev, 2)

iostream = BasicContParamIOStream(
  iostreamsize, iostreamn, [:value, :gradloglikelihood], filepath=filepath, diagnostickeys=[:accept]
)
for i in 1:iostreamn
  state = BasicContMuvParameterState(nstatev[:, i], Symbol[], [:accept], [nstated[i]])
  state.gradloglikelihood = nstategll[:, i]
  write(iostream, state)
end

close(iostream)

iostream = BasicContParamIOStream(
  iostreamsize, iostreamn, [:value, :gradloglikelihood], filepath=filepath, diagnostickeys=[:accept], mode="r"
)
nstate = read(iostream, Float64)

@test isa(nstate, ContMuvMarkovChain{Float64})
@test nstate.value == nstatev
@test nstate.gradloglikelihood == nstategll
for i in [2:4; 6:13]
  @test length(nstate.(fnames[i])) == 0
end
@test length(nstate.diagnostickeys) == 1
@test nstate.diagnosticvalues == nstated
@test nstate.size == iostream.size[1]
@test nstate.n == iostream.n

close(iostream)
for i in (1, 5, 14); rm(filenames[i]); end

println("      Interaction with ContMuvMarkovChain...")

nstatev = Float32[-1.85 -0.09 0.36; -0.45 -0.85 1.91]
nstatell = Float32[-1.30, -1.65, -0.18]
nstatelt = Float32[-0.44, 0.72, -0.21]
nstated = Any[false true; true false; true true]'
iostreamsize = (size(nstatev, 1),)
iostreamn = size(nstatev, 2)

iostream = BasicContParamIOStream(
  iostreamsize, iostreamn, [:value, :loglikelihood, :logtarget], filepath=filepath, diagnostickeys=[:accept]
)
nstatein = ContMuvMarkovChain(iostreamsize[1], iostreamn, [:value, :loglikelihood, :logtarget], [:accept], Float32)
nstatein.value = nstatev
nstatein.loglikelihood = nstatell
nstatein.logtarget = nstatelt
nstatein.diagnosticvalues = nstated
write(iostream, nstatein)

close(iostream)

iostream = BasicContParamIOStream(
  iostreamsize, iostreamn, [:value, :loglikelihood, :logtarget], filepath=filepath, diagnostickeys=[:accept], mode="r"
)
nstateout = read(iostream, Float32)

@test isa(nstateout, ContMuvMarkovChain{Float32})
@test nstateout.value == nstatein.value
@test nstateout.loglikelihood == nstatein.loglikelihood
@test nstateout.logtarget == nstatein.logtarget
for i in [3; 5:13]
  @test length(nstateout.(fnames[i])) == 0
end
@test length(nstateout.diagnostickeys) == 1
@test nstateout.diagnosticvalues == nstatein.diagnosticvalues
@test nstateout.n == nstatein.n

close(iostream)
for i in (1, 2, 4, 14); rm(filenames[i]); end


text:
using Pkg
Pkg.activate(@__DIR__)

# parse command line arguments if given
if length(ARGS) > 0
    subdir = ARGS[1]
# otherwise prompt user to specify
else
    print("Please enter sub directory to run RDE in: ")
    subdir = readline()
end
# input validation
while !isdir(joinpath(@__DIR__, subdir))
    print("Invalid directory $subdir. Please enter sub directory to run RDE in: ")
    global subdir = readline()
end

using PyCall
pushfirst!(PyVector(pyimport("sys")["path"]), joinpath(@__DIR__, subdir))

import FrankWolfe
include("custom_oralces.jl")
include(joinpath(@__DIR__, subdir, "config_birkhoff.jl"))  # load indices, rates, max_iter
cd(subdir)

# Get the Python side of RDE
rde = pyimport("rde")

for idx in indices

    # Load data sample and distortion functional
    x, fname = rde.get_data_sample(idx)
    f, df, node, pred = rde.get_distortion(x)

    # Setup LMO
    lmo = FrankWolfe.BirkhoffPolytopeLMO()

    # helper functions for prototype ordering vectors
    lin_p = convert(Vector{eltype(x)}, LinRange(0.0, 1.0, length(x)))

    function k_to_p(k)
        pk = zeros(eltype(x), length(x))
        pk[end-k+1:end] .= 1.0
        return pk
    end

    # Wrap objective and gradient functions
    function func(S, k)
        if (S isa Matrix) && !(S isa Matrix{eltype(x)})
            S = convert(Matrix{eltype(x)}, S)
        end
        return f(S*k_to_p(k))
    end

    function grad!(storage, S, k)
        if (S isa Matrix) && !(S isa Matrix{eltype(x)})
            S = convert(Matrix{eltype(x)}, S)
        end
        pk = k_to_p(k)
        g = vec(df(S*pk))
        if any(isnan, g) || any(isnan, pk)
            @info "Warning: Numerical instabilities, skipped rate k=$k"
        else
            # storage += g*transpose(pk), but faster with BLAS?
            BLAS.ger!(convert(eltype(x), 1.0), g, pk, storage)
        end
        return storage
    end

    f_stoch = FrankWolfe.StochasticObjective(func, grad!, 1:length(x), zeros(eltype(x), length(x), length(x)))

    # Run FrankWolfe
    println("Running sample $idx")
    S0 = Matrix{eltype(x)}(I, length(x), length(x))

    @time S, V, primal, dual_gap = FrankWolfe.stochastic_frank_wolfe(
        f_stoch,
        lmo,
        S0,
        ;fw_arguments...
    )
    # reset adaptive step size if necessary
    if fw_arguments.line_search isa FrankWolfe.MonotonousNonConvexStepSize
        fw_arguments.line_search.factor = 0
    end

    # Store results
    all_s = zeros(eltype(x), (length(rates), length(x)))
    for (ridx, rate) in enumerate(rates)
        all_s[ridx, :] = S*k_to_p(rate)
        rde.store_single_result(all_s[ridx,:], idx, fname, rate)
    end

    # Store multiple rate results
    rde.store_collected_results(all_s, idx, node, pred, fname, rates, nothing, S, S*lin_p)

end

max_stars_repo_name: jmaces/fw-rde

max_issues_repo_name: ZIB-IOL/fw-rde

max_forks_repo_name: ZIB-IOL/fw-rde


text:
##################################################
# Test star DF
##################################################

"""
    _HdHTest(a,b,c,d,x)

Computes the `H(a,b,c,d;x)` function used in the distribution function and its derivative.
Returns a tuple `(H, dH/dx, d2H/dx2)`.

# Arguments
- `x::Float64`: `x` parameter.
"""
function _HdHTest(a::Float64, b::Float64, c::Float64, d::Float64, x::Float64)
    if (x <= 1)
        pref = 1/(GAMMA_ca*GAMMA_ad)
        HG   = _₂F₁(a+b,1+a-c,a+d,x)#tabHGInt[1](x) #
        HGp  = _₂F₁(a+b+1,a-c+2,a+d+1,x)#tabHGInt[3](x) #
        HGpp = _₂F₁(a+b+2,a-c+3,a+d+2,x)#tabHGInt[5](x) #
        H = pref*x^a*HG
        dH = pref*x^(a-1)*(a*HG + (a+b)*(1+a-c)*x*HGp/(a+d))
        d2H = (pref*x^(a-2)*(a*(a-1)*HG + (2*a)*(a+b)*(1+a-c)*x*HGp/(a+d)
               + ((a+b)*(1+a-c)/(a+d))*((a+b+1)*(2+a-c)/(a+d+1)) *x^2*HGpp))
        return H, dH, d2H
    else
        pref = 1/(GAMMA_db*GAMMA_bc)
        HG   = _₂F₁(a+b,1+b-d,b+c,1/x)#tabHGInt[2](1/x) #
        HGp  = _₂F₁(a+b+1,b-d+2,b+c+1,1/x)#tabHGInt[4](1/x) #
        HGpp = _₂F₁(a+b+2,b-d+3,b+c+2,1/x)#tabHGInt[6](1/x) #
        H = pref*x^(-b)*HG
        dH = pref*x^(-b-1)*((-b)*HG - (1/x)*HGp*(a+b)*(1+b-d)/(b+c))
        d2H = (pref*x^(-b-2)*(b*(b+1)*HG + (2*b+2)*(a+b)*(1+b-d)*(1/x)*HGp/(b+c)
               + ((a+b)*(1+b-d)/(b+c))*((a+b+1)*(2+b-d)/(b+c+1)) *(1/x^2)*HGpp))
        return H, dH, d2H
    end
end

"""
    _tFdF(tE,tL,[q=qCalc])

Gradients of the non-dimensional distribution function of a Plummer sphere.
Returns a tuple `(dtF/dtE, dtF/dtL, d2tF/dtE2, d2tF/dtEtL, d2tF/dtL2)`.

# Arguments
- `tE::Float64`: Reduced energy `tE = E/E0`.
- `tL::Float64`: Reduced energy `tL = L/L0`.
- `q ::Float64`: Anisotropy parameter.
"""
function _tdFTest(tE::Float64, tL::Float64, q::Float64=qCalc)
    if (tE <= 0.0 || tL <= 0.0) # If E or L are negative, the DF vanishes
        dE = 0.0
        dL = 0.0
        dEE = 0.0
        dEL = 0.0
        dLL = 0.0
        return dE, dL, dEE, dEL, dLL
    else
        x = tL^2/(2*tE)
        if (q == 0.0) # Isotropic case
            dE = 3/(PI^3) * (2*tE)^(5/2)
            dL = 0.0
            dEE = 15/(PI^3) * (2*tE)^(3/2)
            dEL = 0.0
            dLL = 0.0
            return dE, dL, dEE, dEL, dLL
        else
            if (q == 2.0)
                if (x <= 1.0)
                    pref = 18/(2*PI)^3
                    dE = pref * (2*tE-tL^2)^(1/2)
                    dL = -pref*tL * (2*tE-tL^2)^(1/2)
                    dEE = pref* (2*tE-tL^2)^(-1/2)
                    dEL = -pref*tL * (2*tE-tL^2)^(-1/2)
                    dLL = -pref*((2*tE-tL^2)^(1/2) - tL^2*(2*tE-tL^2)^(-1/2))
                    # here L ???
                    return dE, dL, dEE, dEL, dLL
                else
                    dE = 0.0
                    dL = 0.0
                    dEE = 0.0
                    dEL = 0.0
                    dLL = 0.0
                    return dE, dL, dEE, dEL, dLL
                end
            else
                pref = 3*GAMMA_6q/(2*(2*PI)^(5/2))
                pref1 = pref*tE^(5/2-q)
                pref2 = pref*tE^(3/2-q)
                H, dH, d2H = _HdHTest(0.0,q/2,9/2-q,1.0,x)
                dE = pref1*((7/2-q)*H - x*dH)
                dL = pref1*(tL*dH)
                dEE = pref2*((7/2-q)*(5/2-q)*H - (5-2*q)*x*dH + x^2*d2H)
                dEL = pref2*tL*((5/2-q)*dH - x*d2H)
                dLL = pref2*(tE*dH +tL^2*d2H)
                return dE, dL, dEE, dEL, dLL
            end
        end
    end
end

"""
    _dTestF(E,L,[q=qCalc])

Gradients of the distribution function of a Plummer sphere.
Normalized to the total mass of the cluster `_M`.
Returns a tuple `(dF/dE, dF/dL, d2F/dE2, d2F/dEL, d2F/dL2)`.

# Arguments
- `E::Float64`: Energy.
- `L::Float64`: Angular momentum.
- `q ::Float64`: Anisotropy parameter.
"""
function _dFTest(E::Float64, L::Float64, q::Float64=qCalc)
    tE = _tE(E)
    tL = _tL(L)
    tdE, tdL, tdEE, tdEL, tdLL = _tdFTest(tE,tL,q)
    dE = (_M*_F0/_E0)*tdE
    dL = (_M*_F0/_L0)*tdL
    dEE = (_M*_F0/(_E0^2))*tdEE
    dEL = (_M*_F0/(_E0*_L0))*tdEL
    dLL = (_M*_F0/(_L0^2))*tdLL
    return dE, dL, dEE, dEL, dLL
end

##################################################
# Compute numerical-gradient flux and its divergence
##################################################

"""
    divflux_NR_num(Jr,L,m_field,[eps=10^(-5)*_L0,nbK=nbK_default,nbAvr=nbAvr_default,m_test=m_field])

Computes the divergence of the flux in action space.

# Remarks
- `dF/dt = -div(flux)`.

# Arguments
- `Jr::Float64`: Radial action parameter.
- `L::Float64`: Angular momentum parameter.
- `m_field::Float64`: Mass of field star.
- `eps::Float64`: Finite difference used in numerical gradient.
- `nbK::Int64`: Sampling number of Rosenbluth integrals.
- `nbAvr::Int64`: Sampling number of orbit average.
- `m_test::Float64`: Mass of test star.
"""
function divflux_NR_num(Jr::Float64, L::Float64, m_field::Float64, eps::Float64=10^(-5)*_L0 ,
                    nbK::Int64=nbK_default, nbAvr::Int64=nbAvr_default,
                    m_test::Float64=m_field)

    Jr_p = Jr+eps
    Jr_m = Jr-eps

    L_p = L+eps
    L_m = L-eps

    E = _E_from_Jr(Jr,L)
    E_Lp = _E_from_Jr(Jr,L_p,nbAvr)
    E_Lm = _E_from_Jr(Jr,L_m,nbAvr)
    E_p = _E_from_Jr(Jr_p,L,nbAvr)
    E_m = _E_from_Jr(Jr_m,L,nbAvr)
    E_p_Lp = _E_from_Jr(Jr_p,L_p,nbAvr)
    E_p_Lm = _E_from_Jr(Jr_p,L_m,nbAvr)
    E_m_Lp = _E_from_Jr(Jr_m,L_p,nbAvr)
    E_m_Lm = _E_from_Jr(Jr_m,L_m,nbAvr)

    # Distribution function
    Ftot = _F(E,L)
    dFtotdE, dFtotdL, d2FtotdE2, d2FtotdEL, d2FtotdL2 = _dFTest(E,L)

    # Compute dJr/dE and dJr/dL
    sp, sa = radius_s_bounds_action(Jr,L,nbAvr)
    a = sma_eff(sp,sa)
    ecc = ecc_eff(sp,sa)

    dJrdE = 0.0
    dJrdL = 0.0
    d2JrdE2 = 0.0
    d2JrdEL = 0.0
    d2JrdL2 = 0.0

    for iu=1:nbAvr
        uloc = -1+2*(iu-0.5)/nbAvr
        rloc = r_eff_anomaly(uloc,a,ecc)
        trloc = rloc/_b
        sloc = sqrt(trloc^2+1)

        jac_loc = dth1du(uloc,sp,sa)
        djacdE_loc, djacdL_loc, dsdL_loc = djac_and_ds(uloc,sp,sa,L)

        dJrdE += jac_loc
        dJrdL += jac_loc/rloc^2
        d2JrdE2 += djacdE_loc
        d2JrdEL += djacdL_loc
        d2JrdL2 += (jac_loc/trloc^2+L/trloc^4*
                (djacdL_loc*trloc^2-2*jac_loc*sloc*dsdL_loc))

    end

    dJrdE *= (2/nbAvr)*(1/_Omega0)/PI
    dJrdL *= (2/nbAvr)*(-L/PI)*(1/_Omega0)
    d2JrdE2 *= (2/nbAvr)*(1/_Omega0)/PI
    d2JrdEL *= (2/nbAvr)*(1/_Omega0)/PI
    d2JrdL2 *= (2/nbAvr)*(-1/PI)/_L0

    # Compute DFtot/DJr and DFtot/DL
    dEdL = -dJrdL/dJrdE
    dEdJr = 1/dJrdE

    d2EdJr2 = -d2JrdE2/(dJrdE)^3
    d2EdJrL = -(d2JrdE2*dEdL+d2JrdEL)/(dJrdE)^2
    d2EdL2 = -((dEdL*d2JrdE2+d2JrdEL)*dEdL+dEdL*d2JrdEL+d2JrdL2)/(dJrdE)

    DFtotDJr = (1/dJrdE)*dFtotdE
    DFtotDL = dEdL*dFtotdE + dFtotdL

    D2FtotDJr2 = d2EdJr2*dFtotdE+(dEdJr)^2*d2FtotdE2
    D2FtotDJrL = d2EdJrL*dFtotdE + dEdJr*(dEdL*d2FtotdE2+d2FtotdEL)
    D2FtotDL2 = d2EdL2*dFtotdE + dEdL*(dEdL*d2FtotdE2+d2FtotdEL)+dEdL*d2FtotdEL+d2FtotdL2

    # Diffusion coefficients
    DJr, DL, DJrJr, DJrL, DLL = avr_action_coefficients(Jr,L,m_field,nbK,nbAvr,m_test)
    DJr_p, _, DJrJr_p, DJrL_Jrp, _ = avr_action_coefficients(Jr_p,L,m_field,nbK,nbAvr,m_test)
    DJr_m, _, DJrJr_m, DJrL_Jrm, _ = avr_action_coefficients(Jr_m,L,m_field,nbK,nbAvr,m_test)
    _, DL_p, _, DJrL_Lp, DLL_p = avr_action_coefficients(Jr,L_p,m_field,nbK,nbAvr,m_test)
    _, DL_m, _, DJrL_Lm, DLL_m = avr_action_coefficients(Jr,L_m,m_field,nbK,nbAvr,m_test)
    _, _, _, DJrL_pp, _ = avr_action_coefficients(Jr_p,L_p,m_field,nbK,nbAvr,m_test)
    _, _, _, DJrL_mm, _ = avr_action_coefficients(Jr_m,L_m,m_field,nbK,nbAvr,m_test)
    _, _, _, DJrL_pm, _ = avr_action_coefficients(Jr_p,L_m,m_field,nbK,nbAvr,m_test)
    _, _, _, DJrL_mp, _ = avr_action_coefficients(Jr_m,L_p,m_field,nbK,nbAvr,m_test)

    # d(DJrF)/dJr
    termJr1 = (2*L*DJr_p-2*L*DJr_m)/(2*eps*_L0)*Ftot
    termJr2 = 2*L*DJr*DFtotDJr
    termJr = termJr1 + termJr2

    # d(DLF)/dL
    termL1 = (2*L_p*DL_p-2*L_m*DL_m)/(2*eps)*Ftot
    termL2 = 2*L*DL*DFtotDL
    termL = termL1 + termL2

    # d2(DJrJrF)/dJr2
    termJrJr1 = (2*L*DJrJr_p+2*L*DJrJr_m-2*2*L*DJrJr)/((eps)^2)*Ftot
    termJrJr2 = 2*(2*L*DJrJr_p-2*L*DJrJr_m)/(2*eps)*DFtotDJr
    termJrJr3 = 2*L*DJrJr*D2FtotDJr2
    termJrJr = termJrJr1 + termJrJr2 + termJrJr3

    # d2(DLLF)/dL2
    termLL1 = (2*L_p*DLL_p+2*L_m*DLL_m-2*2*L*DLL)/((eps)^2)*Ftot
    termLL2 = 2*(2*L_p*DLL_p-2*L_m*DLL_m)/(2*eps)*DFtotDL
    termLL3 = 2*L*DLL*D2FtotDL2
    termLL = termLL1 + termLL2 + termLL3

    # d2(DJrLF)/dJrL
    termJrL1 = (2*L_p*DJrL_pp+2*L_m*DJrL_mm-2*L_m*DJrL_pm-2*L_p*DJrL_mp)/(4*(eps)^2)*Ftot
    termJrL2 = (2*L*DJrL_Jrp-2*L*DJrL_Jrm)/(2*eps)*DFtotDL
    termJrL3 = (2*L_p*DJrL_Lp-2*L_m*DJrL_Lm)/(2*eps)*DFtotDJr
    termJrL4 = 2*L*DJrL*D2FtotDJrL
    termJrL = termJrL1 + termJrL2 + termJrL3 + termJrL4

    return termJr+termL-0.5*termJrJr-0.5*termLL-termJrL
end

"""
    flux_NR_num(Jr,L,m_field,[eps=10^(-5)*_L0,nbK=nbK_default,nbAvr=nbAvr_default,m_test=m_field])

Computes the flux in action space.

# Arguments
- `Jr::Float64`: Radial action parameter.
- `L::Float64`: Angular momentum parameter.
- `m_field::Float64`: Mass of field star.
- `eps::Float64`: Finite difference used in numerical gradient.
- `nbK::Int64`: Sampling number of Rosenbluth integrals.
- `nbAvr::Int64`: Sampling number of orbit average.
- `m_test::Float64`: Mass of test star.
"""
function flux_NR_num(Jr::Float64, L::Float64, m_field::Float64, eps::Float64=10^(-5)*_L0,
                    nbK::Int64=nbK_default, nbAvr::Int64=nbAvr_default,
                    m_test::Float64=m_field)

        Jr_p = Jr+eps
        Jr_m = Jr-eps

        L_p = L+eps
        L_m = L-eps

        E = _E_from_Jr(Jr,L)

        Ftot, dFtotdE, dFtotdL = _FdF(E,L)

        # Compute dJr/dE and dJr/dL
        sp, sa = radius_s_bounds_action(Jr,L,nbAvr)
        a = sma_eff(sp,sa)
        ecc = ecc_eff(sp,sa)

        dJrdE = 0.0
        dJrdL = 0.0

        for iu=1:nbAvr
            uloc = -1+2*(iu-0.5)/nbAvr
            rloc = r_eff_anomaly(uloc,a,ecc)

            jac_loc = dth1du(uloc,sp,sa)
            dJrdE += jac_loc
            dJrdL += jac_loc/rloc^2
        end

        dJrdE *= (2/nbAvr)*(1/_Omega0)/PI
        dJrdL *= (2/nbAvr)*(-L/PI)*(1/_Omega0)

        # Compute DFtot/DJr and DFtot/DL
        DFtotDJr = (1/dJrdE)*dFtotdE

        dEdL = -dJrdL/dJrdE
        DFtotDL = dEdL*dFtotdE + dFtotdL

        # Compute diffusion coefficients
        DJr, DL, DJrJr, DJrL, DLL = avr_action_coefficients(Jr,L,m_field,nbK,nbAvr,m_test) #
        _, _, DJrJr_Jrp, DJrL_Jrp, _ = avr_action_coefficients(Jr_p,L,m_field,nbK,nbAvr,m_test) #
        _, _, DJrJr_Jrm, DJrL_Jrm, _ = avr_action_coefficients(Jr_m,L,m_field,nbK,nbAvr,m_test) #
        _, _, _, DJrL_Lp, DLL_Lp = avr_action_coefficients(Jr,L_p,m_field,nbK,nbAvr,m_test) #
        _, _, _, DJrL_Lm, DLL_Lm = avr_action_coefficients(Jr,L_m,m_field,nbK,nbAvr,m_test) #

        # Flux Jr
        termJr1 = 2*L*DJr*Ftot
        termJr2 = -0.5*((2*L*DJrJr_Jrp-2*L*DJrJr_Jrm)/(2*eps))*Ftot
        termJr3 = -L*DJrJr*DFtotDJr
        termJr4 = -0.5*((2*L_p*DJrL_Lp-2*L_m*DJrL_Lm)/(2*eps))*Ftot
        termJr5 = -L*DJrL*DFtotDL

        FJr = termJr1+termJr2+termJr3+termJr4+termJr5

        # Flux Jr
        termL1 = 2*L*DL*Ftot
        termL2 = -0.5*((2*L*DJrL_Jrp-2*L*DJrL_Jrm)/(2*eps))*Ftot
        termL3 = -L*DJrL*DFtotDJr
        termL4 = -0.5*((2*L_p*DLL_Lp-2*L_m*DLL_Lm)/(2*eps))*Ftot
        termL5 = -L*DLL*DFtotDL

        FL = termL1+termL2+termL3+termL4+termL5

        return FJr, FL

end

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
struct RadauIIA3Tableau{T,T2}
  T11::T
  T12::T
  T21::T
  T22::T
  TI11::T
  TI12::T
  TI21::T
  TI22::T
  c1::T2
  c2::T2
  α::T
  β::T
end

function RadauIIA3Tableau(T,T2)
  T11 = T(0.10540925533894596)
  T12 = T(-0.29814239699997197)
  T21 = T(0.9486832980505138)
  T22 = T(0.0)
  TI11 = T(0.0)
  TI12 = T(1.0540925533894598)
  TI21 = T(-3.3541019662496843)
  TI22 = T(0.3726779962499649)
  c1 = T2(1/3)
  c2 = T2(1.0)
  α  = T(2.0)
  β  = T(-sqrt(2))
  RadauIIA3Tableau{T,T2}(T11, T12, T21, T22,
  TI11, TI12, TI21, TI22,
  c1, c2,   α,    β)
end

struct RadauIIA5Tableau{T,T2}
  T11::T
  T12::T
  T13::T
  T21::T
  T22::T
  T23::T
  T31::T
  #T32::T = 1
  #T33::T = 0
  TI11::T
  TI12::T
  TI13::T
  TI21::T
  TI22::T
  TI23::T
  TI31::T
  TI32::T
  TI33::T
  c1::T2
  c2::T2
  #c3::T2 = 1
  γ::T
  α::T
  β::T
  e1::T
  e2::T
  e3::T
end

# inv(T) * inv(A) * T = [γ  0  0
#                        0  α -β
#                        0  β  α]
function RadauIIA5Tableau(T,T2)
  T11  = convert(T, 9.1232394870892942792e-02)
  T12  = convert(T, -0.14125529502095420843e0)
  T13  = convert(T, -3.0029194105147424492e-02)
  T21  = convert(T, 0.24171793270710701896e0)
  T22  = convert(T, 0.20412935229379993199e0)
  T23  = convert(T, 0.38294211275726193779e0)
  T31  = convert(T, 0.96604818261509293619e0)
  TI11 = convert(T, 4.3255798900631553510e0)
  TI12 = convert(T, 0.33919925181580986954e0)
  TI13 = convert(T, 0.54177053993587487119e0)
  TI21 = convert(T, -4.1787185915519047273e0)
  TI22 = convert(T, -0.32768282076106238708e0)
  TI23 = convert(T, 0.47662355450055045196e0)
  TI31 = convert(T, -0.50287263494578687595e0)
  TI32 = convert(T, 2.5719269498556054292e0)
  TI33 = convert(T, -0.59603920482822492497e0)

  sqrt6 = sqrt(6)

  c1 = convert(T2, (4-sqrt6)/10)
  c2 = convert(T2, (4+sqrt6)/10)

  cbrt9 = cbrt(9)

  γ′ = convert(T, (6.0+cbrt9*(cbrt9-1))/30  ) # eigval of `A`
  α′ = convert(T, (12.0-cbrt9*(cbrt9-1))/60 ) # eigval of `A`
  β′ = convert(T, cbrt9*(cbrt9+1)*sqrt(3)/60) # eigval of `A`
  scale = α′^2 + β′^2
  γ  = inv(γ′)  # eigval of `inv(A)`
  α  = α′/scale # eigval of `inv(A)`
  β  = β′/scale # eigval of `inv(A)`

  e1 = convert(T, -(13+7*sqrt6)/3)
  e2 = convert(T, (-13+7*sqrt6)/3)
  e3 = convert(T, -1/3)
  RadauIIA5Tableau{T,T2}(
  T11 , T12 , T13 , T21 , T22 , T23 , T31 , #= T32 = 1 =# #= T33 = 0 =#
  TI11, TI12, TI13, TI21, TI22, TI23, TI31, TI32, TI33,
  c1,   c2, #= c3 = 1 =#
  γ,    α,    β,
  e1,   e2,   e3)
end


text:
"""
    foreachdump(processdump, input::IO)

`processdump` is a function with the following signature

    processdump(processir)

where `processir` is a function with the following signature

    processir(processline) -> name::AbstractString

where `name` is the LLVM level name of the function and `processline` is a
function with the following signature

    processline(line::AbstractString)

which is called for each `line` of the IR.
"""
function foreachdump(processdump, input::IO)
    firstline = Ref{String}()
    while true
        line = readline(input; keep = true)
        if startswith(line, "*** ")
            firstline[] = line
            break
        end
        @warn "Ignoreing:" line
        eof(input) && return
    end
    while true
        processdump() do processline
            processline(firstline[])
            name = nothing
            while true
                line = readline(input; keep = true)
                if name === nothing
                    name = parsename(line)
                end
                if startswith(line, "*** ")
                    firstline[] = line
                    break
                end
                processline(line)
                eof(input) && break
            end
            return name
        end
        eof(input) && return
    end
end

function parsename(line)
    m = match(r"^define .*@([^(]+)", line)
    m === nothing && return nothing
    return m[1]
end

asfilename(::Nothing) = "UNKNOWN"
function asfilename(name::AbstractString)
    name = replace(name, r"^[^a-zA-Z0-9]+" => "")
    name = replace(name, r"[^a-zA-Z0-9]+$" => "")
    name = replace(name, r"[^a-zA-Z0-9]+" => "_")
    return name
end

function tryparse_passname_line(ln::AbstractString)
    m = match(r"^\*\*\* IR Dump After (.*) \*\*\*", ln)
    if m !== nothing
        return m[1]
    end
    return nothing
end

function passof(code::AbstractString)
    for line in eachline(IOBuffer(code))
        name = tryparse_passname_line(line)
        name === nothing || return name
    end
    return nothing
end

stripoff_passname(code::AbstractString) = stripoff_passname(IOBuffer(code))
function stripoff_passname(input::IO)
    output = IOBuffer()
    for line in eachline(input; keep = true)
        if tryparse_passname_line(line) === nothing
            print(output, line)
        end
    end
    return String(take!(output))
end

extract_function(code::AbstractString) = extract_function(IOBuffer(code))
function extract_function(input::IO)
    output = IOBuffer()
    while !eof(input)
        line = readline(input; keep = true)
        if parsename(line) !== nothing
            print(output, line)
            break
        end
    end
    while !eof(input)
        line = readline(input; keep = true)
        print(output, line)
        if match(r"^}\s*$", line) !== nothing
            break
        end
    end
    return String(take!(output))
end

cleanup_ir_line(line::AbstractString) = replace(line, r" *addrspace\([0-9]+\)" => "")

# TODO: Use `opt`
cleanup_ir(code::AbstractString) = cleanup_ir(IOBuffer(code))
function cleanup_ir(input::IO)
    output = IOBuffer()
    for line in eachline(input; keep = true)

        addrspacecast = match(r"^( *%[^ ]* += *)addrspacecast .*(%[^ ]*) +", line)
        if addrspacecast !== nothing
            print(output, addrspacecast[1], addrspacecast[2])
            print(output, " ; ", line)
            continue
        end

        print(output, cleanup_ir_line(line))
    end
    return String(take!(output))
end

IRDumps.parsefile(filepath) = open(IRDumps.parse, filepath)
function IRDumps.parse(io::IO)
    irs = FunctionCollection()
    counter = Ref(0)
    buffer = IOBuffer()
    foreachdump(io) do foreachline
        fnname = foreachline() do ln
            write(buffer, ln)
        end
        code = String(take!(buffer))
        id = counter[] += 1
        ir = IRDump(id, fnname, passof(code), code)
        push!(irs, ir)
    end
    return irs
end


text:
################helper functions for the tables and regressions

function regmbb(X::Vector{TXcol},y; bssamples, bswidth, β=nothing) where TXcol<:AbstractVector

  #draw the coefficients for the bootstrap
  function bsstatfunc(sim)
    y = sim[1]
    X = reduce(hcat, sim[2:end])

    #try
      b = cholesky!(Symmetric(X'*X))\(X'*y)
    #=catch
      println("X[1:5,:]: $(X[1:5,:])
        y[1:5]: $(y[1:5])")
      throw("cholesky factorization failed")

      #return bsstatfunc(sim)
    end=#
    return Dict{Symbol, Vector{Float64}}(
      :beta => b,
      #:betamr => ρ * std(m)/std(r)
      )
  end

  function bsaggfunc(bs)
    bsmat = reduce(hcat,bs)
    aggfunc = (;se=(v->std(skipmissing(v))).(eachrow(bsmat)), p=nonparametricp.(eachrow(bsmat)))
  end

  return mbbootstrap([y, X...,], bsstatfunc, bsaggfunc, Vector{Float64}, B=bssamples, w=bswidth)
end

#helper funciton to handle the likely scenario of an X matrix instead of a vector
regmbb(X::Matrix, y; bssamples, bswidth) = regmbb(X |> eachcol |> collect, y; bssamples, bswidth)

function regnullar1(X::AbstractMatrix,y; β, nullar1simulations, integrated)

  function nullar1statfunc(Xy)
    Xs = @view Xy[:, 1:(end-1)]
    ys = @view Xy[:, end]
    return Dict(:beta => cholesky!(Symmetric(Xs'*Xs))\(Xs'*ys),)
  end

  stats::Dict{Symbol, Vector{Float64}} = Dict(:beta=>β)
  nullar1aggfunc(M,statname) = (;
    se=std(M, dims=1),
    p=nonparametricnullp(M, stats[statname]))

  integrated && return comparearima11undernull(X,y, nullar1statfunc,
    nullar1aggfunc, Nsimulations=nullar1simulations)
  return comparear1undernull(X,y, nullar1statfunc,
    nullar1aggfunc, Nsimulations=nullar1simulations)
end

#####################

#the below function and info is for formatting numbers
function llformat(β;scalefactor, decimals, prefix, suffix)
  valstring = num2str(β, decimals; scalefactor)
  return "$(prefix)$(valstring)$(suffix)"
end

llbetaformat =(; scalefactor=1.0, decimals=2, prefix="", suffix="\\:\\:")
llpercentageformat =(; scalefactor=100.0, decimals=2, prefix="", suffix="\\:\\:")
llseformat = (;scalefactor=1.0, decimals=2, prefix="\\mathit{(", suffix=")}")
llNformat = (;scalefactor=1.0, decimals=0, prefix="", suffix="")

#put all of the field information here
llfields = Dict{Symbol, Any}(
  :measure=>Dict(
    :L1g_WXmomentum121s=>(;
      id=:L1g_WXmomentum121s,
      label="\$\\%gMOM_{t-1}\$",
      format=llpercentageformat,
      se=(;
        id=:L1g_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1p_WXmomentum121s=>(;
      id=:L1p_WXmomentum121s,
      label="\$\\%Par_{t-1}\$",
      format=llpercentageformat,
      se=(;
        id=:L1p_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1ag_WXmomentum121s=>(;
      id=:L1ag_WXmomentum121s,
      label="\\%gMOM (active)",
      format=llpercentageformat,
      se=(;
        id=:L1ag_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1cgg_WXmomentum121s=>(;
      id=:L1cgg_WXmomentum121s,
      label="\\%gMOM (returns)",
      format=llpercentageformat,
      se=(;
        id=:L1cgg_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1GMT_WXmomentum121s=>(;
      id=:L1GMT_WXmomentum121s,
      label="\$GMT_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:L1GMT_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1nGMT_WXmomentum121s=>(;
      id=:L1nGMT_WXmomentum121s,
      label="\$GMT_{t-1}^{net}\$",
      format=llbetaformat,
      se=(;
        id=:L1nGMT_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :dtXL1lMT_WXmomentum121s=>(;
      id=:dtXL1lMT_WXmomentum121s,
      label="dtXL1lMT",
      format=llbetaformat,
      se=(;
        id=:dtXL1lMT_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1GAN_WXmomentum121s=>(;
      id=:L1GAN_WXmomentum121s,
      label="GAN",
      format=llbetaformat,
      se=(;
        id=:L1GAN_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1GEON_WXmomentum121s=>(;
      id=:L1GEON_WXmomentum121s,
      label="GEON",
      format=llbetaformat,
      se=(;
        id=:L1GEON_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1GDEON_WXmomentum121s=>(;
      id=:L1GDEON_WXmomentum121s,
      label="GDEON",
      format=llbetaformat,
      se=(;
        id=:L1GDEON_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1M_WXmomentum121s=>(;
      id=:L1M_WXmomentum121s,
      label="\$M^{L+S}_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:L1M_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1MT_WXmomentum121s=>(;
      id=:L1MT_WXmomentum121s,
      label="\$M^{tot}_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:L1MT_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1AN_WXmomentum121s=>(;
      id=:L1AN_WXmomentum121s,
      label="\$MF^{all}_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:L1AN_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L1DEON_WXmomentum121s=>(;
      id=:L1DEON_WXmomentum121s,
      label="\$MF^{dom. eq.}_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:L1DEON_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :dtXL1M_WXmomentum121s=>(;
      id=:dtXL1M_WXmomentum121s,
      label="\$dtM^{L+S}_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:dtXL1M_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :dtXL1MT_WXmomentum121s=>(;
      id=:dtXL1MT_WXmomentum121s,
      label="\$dtM^{tot}_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:dtXL1MT_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :dtXL1AN_WXmomentum121s=>(;
      id=:dtXL1AN_WXmomentum121s,
      label="\$dtMF^{all}_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:dtXL1AN_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :dtXL1DEON_WXmomentum121s=>(;
      id=:dtXL1DEON_WXmomentum121s,
      label="\$dtMF^{dom. eq.}_{t-1}\$",
      format=llbetaformat,
      se=(;
        id=:dtXL1DEON_WXmomentum121s_se,
        label="",
        format=llseformat)),
      :L12M_WXmomentum121s=>(;
        id=:L12M_WXmomentum121s,
        label="\$M^{L+S}_{t-12}\$",
        format=llbetaformat,
        se=(;
          id=:L12M_WXmomentum121s_se,
          label="",
          format=llseformat)),
      :L12MT_WXmomentum121s=>(;
        id=:L12MT_WXmomentum121s,
        label="\$M^{tot}_{t-12}\$",
        format=llbetaformat,
        se=(;
          id=:L12MT_WXmomentum121s_se,
          label="",
          format=llseformat)),
      :L12AN_WXmomentum121s=>(;
        id=:L12AN_WXmomentum121s,
        label="\$MF^{all}_{t-12}\$",
        format=llbetaformat,
        se=(;
          id=:L12AN_WXmomentum121s_se,
          label="",
          format=llseformat)),
      :L12DEON_WXmomentum121s=>(;
        id=:L12DEON_WXmomentum121s,
        label="\$MF^{dom. eq.}_{t-12}\$",
        format=llbetaformat,
        se=(;
          id=:L12DEON_WXmomentum121s_se,
          label="",
          format=llseformat)),
      :dtXL12M_WXmomentum121s=>(;
        id=:dtXL12M_WXmomentum121s,
        label="\$dtM^{L+S}_{t-12}\$",
        format=llbetaformat,
        se=(;
          id=:dtXL12M_WXmomentum121s_se,
          label="",
          format=llseformat)),
      :dtXL12MT_WXmomentum121s=>(;
        id=:dtXL12MT_WXmomentum121s,
        label="\$dtM^{tot}_{t-12}\$",
        format=llbetaformat,
        se=(;
          id=:dtXL12MT_WXmomentum121s_se,
          label="",
          format=llseformat)),
      :dtXL2AN_WXmomentum121s=>(;
        id=:dtXL12AN_WXmomentum121s,
        label="\$dtMF^{all}_{t-12}\$",
        format=llbetaformat,
        se=(;
          id=:dtXL12AN_WXmomentum121s_se,
          label="",
          format=llseformat)),
      :dtXL12DEON_WXmomentum121s=>(;
        id=:dtXL12DEON_WXmomentum121s,
        label="\$dtMF^{dom. eq.}_{t-12}\$",
        format=llbetaformat,
        se=(;
          id=:dtXL12DEON_WXmomentum121s_se,
          label="",
          format=llseformat)),
    :L1F_ff5_mktrf=>(;
      id=:L1F_ff5_mktrf,
      label="\$\\%Mkt_{t-1}\$",
      format=llpercentageformat,
      se=(;
        id=:L1F_ff5_mktrf_se,
        label="",
        format=llseformat)),
    :L1ng_WXmomentum121s=>(;
      id=:L1ng_WXmomentum121s,
      label="\$\\%gMOM_{t-1}-Mkt_{t-1}\$",
      format=llpercentageformat,
      se=(;
        id=:L1ng_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L2GMT_WXmomentum121s=>(;
      id=:L2GMT_WXmomentum121s,
      label="\$GMT_{t-2}\$",
      format=llbetaformat,
      se=(;
        id=:L2GMT_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L2nGMT_WXmomentum121s=>(;
      id=:L2nGMT_WXmomentum121s,
      label="\$GMT_{t-2}^{net}\$",
      format=llbetaformat,
      se=(;
        id=:L2nGMT_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L2p_WXmomentum121s=>(;
      id=:L2p_WXmomentum121s,
      label="\$\\%Par_{t-2}\$",
      format=llpercentageformat,
      se=(;
        id=:L2p_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L2g_WXmomentum121s=>(;
      id=:L2g_WXmomentum121s,
      label="\$\\%gMOM_{t-2}\$",
      format=llpercentageformat,
      se=(;
        id=:L2g_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L2ng_WXmomentum121s=>(;
      id=:L2ng_WXmomentum121s,
      label="\$\\%gMOM_{t-2}-Mkt_{t-2}\$",
      format=llpercentageformat,
      se=(;
        id=:L2ng_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L2GEON_WXmomentum121s=>(;
      id=:L2GEON_WXmomentum121s,
      label="\$\\frac{A_{t-1}-A_{t-2}}{AUM^{eq}_{t-2}}\$",
      format=llbetaformat,
      se=(;
        id=:L2GEON_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L2F_ff5_mktrf=>(;
      id=:L2F_ff5_mktrf,
      label="\$\\%Mkt_{t-2}\$",
      format=llpercentageformat,
      se=(;
        id=:L2F_ff5_mktrf_se,
        label="",
        format=llseformat)),
    :L3g_WXmomentum121s=>(;
      id=:L3g_WXmomentum121s,
      label="\$\\%gMOM_{t-3}\$",
      format=llpercentageformat,
      se=(;
        id=:L3g_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :L3F_ff5_mktrf=>(;
      id=:L3F_ff5_mktrf,
      label="\$Mkt_{t-3}\$",
      format=llpercentageformat,
      se=(;
        id=:L3F_ff5_mktrf_se,
        label="",
        format=llseformat)),

    ######for contemporaneous return relationships
    :R_WXmomentum121s=>(;
      id=:R_WXmomentum121s,
      label="\$mom_{t}\$",
      format=llbetaformat,
      se=(;
        id=:R_WXmomentum121s_se,
        label="",
        format=llseformat)),
    :WXmomentum121s_retlevgross=>(;
      id=:WXmomentum121s_retlevgross,
      label="\$mom_{t}^{L+S}\$",
      format=llbetaformat,
      se=(;
        id=:WXmomentum121s_retlevgross_se,
        label="",
        format=llseformat)),
    :WXmomentum121s_long=>(;
      id=:WXmomentum121s_long,
      label="\$mom_{t}^{L}\$",
      format=llbetaformat,
      se=(;
        id=:WXmomentum121s_long_se,
        label="",
        format=llseformat)),
    :WXmomentum121s_short=>(;
      id=:WXmomentum121s_short,
      label="\$mom_{t}^{S}\$",
      format=llbetaformat,
      se=(;
        id=:WXmomentum121s_short_se,
        label="",
        format=llseformat)),
    :F_ff3m_umd=>(;
      id=:F_ff3m_umd,
      label="\$umd_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff3m_umd_se,
        label="",
        format=llseformat)),
    :F_ff3m_mktrf=>(;
      id=:F_ff3m_mktrf,
      label="\$mktrf_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff3m_mktrf_se,
        label="",
        format=llseformat)),
    :F_ff3m_smb=>(;
      id=:F_ff3m_smb,
      label="\$smb_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff3m_smb_se,
        label="",
        format=llseformat)),
    :F_ff3m_hml=>(;
      id=:F_ff3m_hml,
      label="\$hml_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff3m_hml_se,
        label="",
        format=llseformat)),
    :F_ff5_mktrf=>(;
      id=:F_ff5_mktrf,
      label="\$MktRf_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff5_mktrf_se,
        label="",
        format=llseformat)),
    :F_ff5_smb=>(;
      id=:F_ff5_smb,
      label="\$smb_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff5_smb_se,
        label="",
        format=llseformat)),
    :F_ff5_hml=>(;
      id=:F_ff5_hml,
      label="\$hml_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff5_hml_se,
        label="",
        format=llseformat)),
    :F_ff5_rmw=>(;
      id=:F_ff5_rmw,
      label="\$rmw_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff5_rmw_se,
        label="",
        format=llseformat)),
    :F_ff5_cma=>(;
      id=:F_ff5_cma,
      label="\$cma_{t}\$",
      format=llbetaformat,
      se=(;
        id=:F_ff5_cma_se,
        label="",
        format=llseformat)),
  ),
  :description=>Dict(
    :Nmonths=>(;
      id=:Nmonths,
      label="\\midrule Months",
      format=llNformat),
    :Nweeks=>(;
      id=:Nweeks,
      label="\\midrule Weeks",
      format=llNformat),
  ),
)

#start of a table spec
tablelt112_3672spec = (;
  tablename="lllt112_3672",
  Fmeasures=[
    [:L1MT_WXmomentum121s],
    [:L12MT_WXmomentum121s],
    [:L1MT_WXmomentum121s],
    [:L12MT_WXmomentum121s],
    [:L1MT_WXmomentum121s],
    [:L12MT_WXmomentum121s],
    ],
    #[:L1g_WXmomentum121s, :L1F_ff5_mktrf]],
  colnames=[["12-month", "36-month", "72-month"], 1:6 .|> i->"($i)",],
  widthcolnames= [[2,2,2],ones(Int,6),],
  Freturns=[[:FW12R_WXmomentum121s for i ∈ 1:2]...;[:FW36R_WXmomentum121s for i ∈ 1:2]...;[:FW72R_WXmomentum121s for i ∈ 1:2]...;],
  FN=:Nmonths,
  Fse=:se,
  alignmentstring="l | rr | rr | rr"
)

#start of a table spec
tableltdt112_3672spec = (;
  tablename="llltdt112_3672",
  Fmeasures=[
    [:dtXL1MT_WXmomentum121s],
    [:dtXL12MT_WXmomentum121s],
    [:dtXL1MT_WXmomentum121s],
    [:dtXL12MT_WXmomentum121s],
    [:dtXL1MT_WXmomentum121s],
    [:dtXL12MT_WXmomentum121s],
    ],
    #[:L1g_WXmomentum121s, :L1F_ff5_mktrf]],
  colnames=[["12-month", "36-month", "72-month"], 1:6 .|> i->"($i)",],
  widthcolnames= [[2,2,2],ones(Int,6),],
  Freturns=[[:dtXFW12R_WXmomentum121s for i ∈ 1:2]...;[:dtXFW36R_WXmomentum121s for i ∈ 1:2]...;[:dtXFW72R_WXmomentum121s for i ∈ 1:2]...;],
  FN=:Nmonths,
  Fse=:se,
  alignmentstring="l | rr | rr | rr"
)

#start of a table spec
tablest11spec = (;
  tablename="llst11",
  Fmeasures=[
    [:L1F_ff5_mktrf],
    [:L1g_WXmomentum121s],
    [:L1GMT_WXmomentum121s],
    [:L1nGMT_WXmomentum121s],
    [:L1ng_WXmomentum121s],
    [:L1GMT_WXmomentum121s, :L1F_ff5_mktrf]],
    #[:L1g_WXmomentum121s, :L1F_ff5_mktrf]],
  Freturns=[[:FW1R_WXmomentum121s for i ∈ 1:6]...;],
  FN=:Nweeks,
  Fse=:se,
  alignmentstring="l | rrr | rrr"
)

tablest21spec = (;
  tablename="llst21",
  Fmeasures=[
    [:L2F_ff5_mktrf],
    [:L2GMT_WXmomentum121s],
    [:L2g_WXmomentum121s],
    [:L2nGMT_WXmomentum121s],
    [:L2ng_WXmomentum121s],
    [:L2GMT_WXmomentum121s, :L2F_ff5_mktrf]],
    #[:L2g_WXmomentum121s, :L2F_ff5_mktrf]],
  Freturns=[[:FW1R_WXmomentum121s for i ∈ 1:6]...;],
  #=[[:FW1R_WXmomentum121s; :FW2R_WXmomentum121s; :FW3R_WXmomentum121s];
    [:FW1R_WXmomentum121s; :FW2R_WXmomentum121s; :FW3R_WXmomentum121s];],=#
  #[[:FW1R_WXmomentum121s for i ∈ 1:6]...;],
  FN=:Nweeks,
  Fse=:se,
  alignmentstring="l | rrr | rrr"
)

tablestpspec = (;
  tablename="llstp",
  Fmeasures=[
    [:L1p_WXmomentum121s],
    [:L1p_WXmomentum121s, :L1F_ff5_mktrf],
    [:L1p_WXmomentum121s],
    [:L1p_WXmomentum121s, :L1F_ff5_mktrf],
    [:L2p_WXmomentum121s],
    [:L2p_WXmomentum121s, :L2F_ff5_mktrf],
    [:L2p_WXmomentum121s],
    [:L2p_WXmomentum121s, :L2F_ff5_mktrf]],
  colnames=[["1-week return","2-week return", "1-week return","2-week return"], 1:8 .|> i->"($i)",],
  widthcolnames= [[2,2,2,2],ones(Int,8),],
    #[:L2g_WXmomentum121s, :L2F_ff5_mktrf]],
  Freturns=[[:FW1R_WXmomentum121s for i ∈ 1:2]...;[:FW2R_WXmomentum121s for i ∈ 1:2]...;
    [:FW1R_WXmomentum121s for i ∈ 1:2]...;[:FW2R_WXmomentum121s for i ∈ 1:2]...;],
  #=[[:FW1R_WXmomentum121s; :FW2R_WXmomentum121s; :FW3R_WXmomentum121s];
    [:FW1R_WXmomentum121s; :FW2R_WXmomentum121s; :FW3R_WXmomentum121s];],=#
  #[[:FW1R_WXmomentum121s for i ∈ 1:6]...;],
  FN=:Nweeks,
  Fse=:se,
  alignmentstring="l | rr | rr | rr | rr"
)

table00spec = (;
  tablename="ll00",
  Fmeasures=[
    [:R_WXmomentum121s],
    [:WXmomentum121s_long],
    [:WXmomentum121s_short],
    #[:WXmomentum121s_retlevgross], #maybe consider breaking this all up into two tables
    [:F_ff3m_umd],
    [:F_ff3m_mktrf],
    [:F_ff3m_smb],
    [:F_ff3m_hml],],
    #[:F_ff3m_umd, :F_ff3m_mktrf, :F_ff3m_smb, :F_ff3m_hml, ]],
  Freturns=[[:g_WXmomentum121s for i ∈ 1:7]...;],
  #=[[:FW1R_WXmomentum121s; :FW2R_WXmomentum121s; :FW3R_WXmomentum121s];
    [:FW1R_WXmomentum121s; :FW2R_WXmomentum121s; :FW3R_WXmomentum121s];],=#
  #[[:FW1R_WXmomentum121s for i ∈ 1:6]...;],
  FN=:Nweeks,
  Fse=:se,
  alignmentstring="l | rrr | rrrr"
)

function lltables(;analysispath = PARAM[:analysispath],)

  bssamples = PARAM[:tabllbsar1samples]
  nullar1simulations = PARAM[:tabllbsar1samples]

  measure00 = CSV.File(
    "$analysispath\\$(PARAM[:tabll00measurefilename])" *
    "\\$(PARAM[:tabll00fundrid])\\$(PARAM[:tabll00fundrid])_leadlag_series.csv"
    ) |> DataFrame
  ll00bswidth = nrow(measure00)^(1/3) |> ceil |> Int
  lltable(measure00, (X,y; β)->regmbb(X,y, bswidth=ll00bswidth; bssamples).beta.se;
    table00spec...)

  measurelt112_3672 = CSV.File(
    "$analysispath\\$(PARAM[:tablllt112_3672measurefilename])" *
    "\\$(PARAM[:tablllt112_3672fundrid])\\$(PARAM[:tablllt112_3672fundrid])_leadlag_series.csv"
    ) |> DataFrame
  lt112_3672bswidth = nrow(measurelt112_3672)^(1/4) |> ceil |> Int
  lltable(measurelt112_3672, (X,y; β)->regnullar1(X,y;β,nullar1simulations,integrated=true).beta.se;
    tablelt112_3672spec...)

  measureltdt112_3672 = CSV.File(
    "$analysispath\\$(PARAM[:tabllltdt112_3672measurefilename])" *
    "\\$(PARAM[:tabllltdt112_3672fundrid])\\$(PARAM[:tabllltdt112_3672fundrid])_leadlag_series.csv"
    ) |> DataFrame
  ltdt112_3672bswidth = nrow(measureltdt112_3672)^(1/4) |> ceil |> Int
  lltable(measureltdt112_3672, (X,y; β)->regnullar1(X,y;β,nullar1simulations,integrated=false).beta.se;
    tableltdt112_3672spec...)

  measurest11 = CSV.File(
    "$analysispath\\$(PARAM[:tabllst11measurefilename])" *
    "\\$(PARAM[:tabllst11fundrid])\\$(PARAM[:tabllst11fundrid])_leadlag_series.csv"
    ) |> DataFrame
  st11bswidth = nrow(measurest11)^(1/3) |> ceil |> Int
  lltable(measurest11, (X,y; β)->regmbb(X,y, bswidth=st11bswidth; bssamples).beta.se;
    tablest11spec...)

  measurest21 = CSV.File(
    "$analysispath\\$(PARAM[:tabllst21measurefilename])" *
    "\\$(PARAM[:tabllst21fundrid])\\$(PARAM[:tabllst21fundrid])_leadlag_series.csv"
    ) |> DataFrame
  st21bswidth = nrow(measurest21)^(1/3) |> ceil |> Int
  lltable(measurest21, (X,y; β)->regmbb(X,y, bswidth=st21bswidth; bssamples).beta.se;
    tablest21spec...)

  measurestp = CSV.File(
    "$analysispath\\$(PARAM[:tabllstpmeasurefilename])" *
    "\\$(PARAM[:tabllstpfundrid])\\$(PARAM[:tabllstpfundrid])_leadlag_series.csv"
    ) |> DataFrame
  stpbswidth = nrow(measurestp)^(1/3) |> ceil |> Int
  lltable(measurestp, (X,y; β)->regmbb(X,y, bswidth=stpbswidth; bssamples).beta.se;
    tablestpspec...)

end

#note the se function can be anything that takes X and y, so it cna return multpile stats
#this is not meant to be high performance
function llregression(measure,se::Tse; FX,Fy) where Tse<:Function
  #construct the info for each row

  cleanmeasure = view(measure, completecases(measure, [FX; Fy]),:)
  #form X from the measure DF + an intercept column
  X = hcat(cleanmeasure[:, FX] |> Matrix{Float64}, ones(nrow(cleanmeasure)))
  y = cleanmeasure[:, Fy] |> Vector{Float64}

  #compute the estimate and standard error
  β = svd(X)\y
  @assert (FMLM(X,y).β .≈ β) |> all

  se = se(X, y; β)

  N = nrow(cleanmeasure)

  t = β ./ se

  return (; FX, β, t, se, N)
end

#create a table for the lead/lag
function lltable(measure, se::Tse;
    tablename,
    Fmeasures,
    Freturns::TFreturns,
    FN,
    Fse,
    tablepath = PARAM[:tablepath],
    colnames=[1:length(Fmeasures) .|> i->"($i)",],
    widthcolnames= [ones(Int,length(Fmeasures)),],
    alignmentcolnames = [
      (rc-> "c").(colnamerow) for (r,colnamerow) ∈ enumerate(colnames)],
    alignmentstring = string("l",["r" for i ∈ 1:length(Fmeasures)]...),

    #construct the info for each row
    #first extract lightly constructed field objects IN ORDER
    datafields = reduce(vcat, (reduce(vcat, Fmeasures) |> unique) .|> (
      (s)->[llfields[:measure][s], llfields[:measure][s][Fse]])),
    descriptionfields = [
      llfields[:description][FN]
    ],

    ) where {TFreturns, Tse<:Function}

  #get the row ids IN ORDER
  Frowids = (f->f.id).([datafields; descriptionfields])

  #create an index of row information. order is not preserved.
  Frowindex = Dict(f.id=>f for f ∈ [datafields; descriptionfields])

  #double check alignment
  @assert (id->Frowindex[id].id≡id).(Frowids) |> all

  #extract the row labels, using the order preserved in Frowid
  descrownames = (id->Frowindex[id].label).(Frowids)

  #do this in case I want to use a different parsing routine etc
  RHS = Fmeasures
  lhs::Vector{Symbol} = TFreturns <: Symbol ? [Freturns for r ∈ 1:size(RHS,1)] : Freturns

  @assert length(lhs) ≡ size(RHS,1)

  #NOTE: each ROW in the dataframe corresponds to a COLUMN in the table
  #cols preserve order, though I will try not to rely on it
  regressioncols = DataFrame([id=>String[] for id ∈ Frowids])
  for (FX, Fy) ∈ zip(RHS, lhs)
    reg = llregression(measure, se; FX, Fy)

    regressioncol = Dict{Symbol, Any}()
    for (i,Fx) ∈ enumerate(FX)

      #pull the field information from the GLOBAL VARIABLE llfields
      estimatefield=llfields[:measure][Fx]
      @assert estimatefield.id ≡ Fx #double check the field is correct
      regressioncol[estimatefield.id] = llformat(reg.β[i]; estimatefield.format...)

      sefield = llfields[:measure][Fx][Fse]
      @assert sefield.id ≡ Symbol(Fx, :_, Fse)
      regressioncol[sefield.id] = llformat(reg.t[i]; sefield.format...)
    end

    #compute N, and eventually any other descriptive fields
    Nfield= llfields[:description][FN]
    @assert Nfield.id ≡ FN
    regressioncol[Nfield.id] = llformat(reg.N; Nfield.format...)

    #make sure we captured the expected fields
    #note- below is UNORDERED
    @assert issetequal(keys(regressioncol)|>collect,
      [FX; FX .|> (s)->llfields[:measure][s][Fse].id; Nfield.id]) "
        Unexpected and/or missing fields detected in lltable!:
        keys(regressioncol)|>collect: $(keys(regressioncol)|>collect)
        [...]: $([FX; FX .|> (s)->llfields[:measure][s][Fse].id; Nfield])
        "

    #all other columns and rows are blank
    blankfields = Dict(f => "" for f ∈ setdiff(propertynames(regressioncols), keys(regressioncol)))

    #form the regression col
    push!(regressioncols, merge(regressioncol, blankfields))
  end

  Ncols = length(Fmeasures)
  #desccontent = [Vector{String}(undef, Ncols) for i ∈ 1:Nrows]

  @assert (propertynames(regressioncols) .≡ Frowids) |> all
  desccontent = eachcol(regressioncols) |> collect .|> Vector

  t = textable(;colnames, descrownames, desccontent,
    widthcolnames, alignmentstring, alignmentcolnames)

  writetextable(t, path=tablepath, outname="$tablename.tex")
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Test
using Random
using LogicCircuits
using ProbabilisticCircuits

@testset "Expectation Maximization algorithm test" begin
    Random.seed!(1337)
    data = dataset(twenty_datasets("nltcs"); do_shuffle=false, batch_size=-1)
    train_x = train(data)

    pc1 = zoo_psdd("nltcs.clt.psdd")
    pc2 = zoo_psdd("nltcs.10split.psdd")
    pcs = [pc1, pc2]
    mf = train_mixture(pcs, convert(XBatches, train_x), 1.0, 3)

    @test mf isa FlatMixtureWithFlow
    @test mf.origin isa FlatMixture
    @test components(mf) == pcs
    @test component_weights(mf) ≈ [0.41093382427863423, 0.5890661757213657] atol=1.0e-9
    @test log_likelihood(mf, convert(XBatches, train_x)) / num_examples(train_x) ≈ -6.475868465920385 atol=1.0e-9
end


text:
using TransformVariables

export asparams

"""
`asparams` build on `TransformVariables.as` to construct bijections to the
parameter space of a given parameterized measure. Because this is only possible
for continuous parameter spaces, we allow constraints to assign values to any
subset of the parameters.

--------

    asparams(::Type{<:ParameterizedMeasure}, ::Val{::Symbol})

Return a transformation for a particular parameter of a given parameterized
measure. For example,

```
julia> asparams(Normal, Val(:σ))
asℝ₊
```

-----

    asparams(::Type{<: ParameterizedMeasure{N}}, constraints::NamedTuple) where {N}

Return a transformation for a given parameterized measure subject to the named tuple
`constraints`. For example,

```
julia> asparams(Binomial{(:p,)}, (n=10,))
TransformVariables.TransformTuple{NamedTuple{(:p,), Tuple{TransformVariables.ScaledShiftedLogistic{Float64}}}}((p = as𝕀,), 1)
```

------------

    aspararams(::ParameterizedMeasure)

Return a transformation with no constraints. For example,

```
julia> asparams(Normal{(:μ,:σ)})
TransformVariables.TransformTuple{NamedTuple{(:μ, :σ), Tuple{TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}}}}((μ = asℝ, σ = asℝ₊), 2)
```
"""
function asparams end

asparams(μ::ParameterizedMeasure, v::Val) = asparams(constructor(μ), v)
asparams(μ, s::Symbol) = asparams(μ, Val(s))

asparams(M::Type{A}) where {A<:AbstractMeasure} = asparams(M, NamedTuple())

function asparams(::Type{M}, constraints::NamedTuple{N}) where {N, M<: ParameterizedMeasure}
    # @show M
    thekeys = paramnames(M, constraints)
    t1 = NamedTuple{thekeys}(asparams(M, Val(k)) for k in thekeys)
    t2 = NamedTuple{N}(map(asConst, values(constraints)))
    C = constructorof(M)
    # @show C
    # @show constraints
    # @show transforms
    # Make sure we end up with a consistent ordering
    ordered_transforms = params(C(merge(t1, t2)))
    return TV.as(ordered_transforms)
end

# TODO: Make this work for Affine measures
# function asparams(::Type{M}, constraints::NamedTuple{N}) where {N, M<: Affine}
#     ...
# end

asparams(μ::ParameterizedMeasure, nt::NamedTuple=NamedTuple()) = asparams(constructor(μ), nt)

TV.as(::Half) = asℝ₊


text:
using Distributed
addprocs(24)
@everywhere using ProgressMeter

@everywhere using Lasso
using NPZ
@everywhere using Statistics
using PyCall
@everywhere using JLD
@everywhere using SparseArrays

py"""
import scipy.sparse as sp

def spload(x):
    return sp.load_npz(x)

def spsave(x):
    return sp.save_npz("Results/results.npz", x)
"""

PyObject(S::SparseMatrixCSC) =
    pyimport("scipy.sparse").csc_matrix((S.nzval, S.rowval .- 1, S.colptr .- 1), shape=size(S))

# G = py"spload"("Data/Gcsc.npz")
# G = SparseMatrixCSC(G.shape[1], G.shape[2], G.indptr .+ 1, G.indices .+ 1, G.data)

G = npzread("Data/G.npy")

@everywhere G = $G

@everywhere function fit_model(G,dv,i)
    if isfile("Results/Outputs/output$i.jld") == false
        dn = std(dv)
        dv ./= dn
        path = fit(LassoPath, G, dv, α=0.95, λminratio = 1e-4, standardize=false, intercept=false)
        #λmin = path.λ[argmin(aicc(path))]
        # m = coef(path; select=MinAICc())
        m = coef(path; select=MinCVmse(path))
        #save("Results/Outputs/output$i.jld", Dict("i"=>i,"reg"=>λmin, "dn"=>dn, "m"=>m))
        save("Results/Outputs/output$i.jld", Dict("i"=>i, "dn"=>dn, "m"=>m))
    end
end

@everywhere fit_function(dv,i) = fit_model(G,dv,i)

d = npzread("Data/wvtdata.npy")
d = [(i, d[i,:]) for i = 1:size(d,1)];

p = Progress(length(d), dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)
progress_pmap(d, progress=p) do x
    fit_function(x[2], x[1])
end

results = [load("Results/Outputs/output$i.jld") for i = 1:size(d,1)]

sres = reduce(hcat, [x["m"].*x["dn"] for x in results])
py"spsave"(PyObject(sres))

# save("output.jld", "results", results)

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
include("ZeroHalfCuts.jl")

module ZeroHalfOnIP

using JuMP
using CPLEX
using LightGraphs
using SimpleWeightedGraphs
using SparseArrays

using ..ZeroHalfChvatGom
using ..readLP

function writeDOTFile(g::SimpleWeightedDiGraph, filename)
    open(filename, "w") do f
        w = weights(g)
        write(f, "digraph G {\n")
        write(f, "\toverlap=scale;\n")
      	write(f, "\tsplines=true;\n")

        n = nv(g)
        for e in edges(g)
            #println("(",src(e),",",dst(e),")" )
            println(f, "\t",src(e)," -> ", dst(e), "[",
                    "label = \"", w[src(e),dst(e)], "\"];")
        end
        write(f, "}\n")
     end
end

function findGraph(A,consLB, consUB)
    # go through the constraints. Any ... = 0 constraint will give rise to a node in the graph
    (m,n) = size(A)
    Atranspose = sparse(A')
    nNodes = 0
    # each variable corresponds to an arc.
    # the head an tail variables define each arc
    head = zeros(Int64, n)
    tail = zeros(Int64, n)

    rows = rowvals(Atranspose)
    vals = nonzeros(Atranspose)
    for i=1:m
        if consLB[i] == consUB[i] && consLB[i] != 1
            #println("Constraint $i ... ")
            nNodes += 1
            # go through  row i of A (the same as column i in A^T )
            for j in nzrange(Atranspose, i)
                varID = rows[j]
                val = vals[j]
                if val != 1 && val != 0 && val != -1
                    throw("Unexpected value in flow assumed conservation constraint")
                else
                    if val < 0
                        if tail[varID] != 0
                            throw("constraint $i: arc $varID already has a tail")
                        end
                        tail[varID] = nNodes
                    elseif val > 0
                        if head[varID] != 0
                            throw("constraint $i: arc $varID already has a head")
                        end
                        head[varID] = nNodes
                    end
                end
            end
        end
    end
    g = SimpleWeightedDiGraph(nNodes+2)
    for i=1:n
        if tail[i] == 0
            throw("arc $i has no tail")
        elseif head[i] == 0
            throw("arc $i has no head")
        else
            add_edge!(g, tail[i], head[i], i)
        end
    end
    writeDOTFile(g, "graph.dot")
end

function test(lpFileName, maxIter=1, strength=3, printCuts=false)
    A, obj, varLB, varUB, consLB, consUB, varType = readLP.loadLP(lpFileName)
    findGraph(A,consLB, consUB)
    zhOrigSubSys = ZeroHalfChvatGom.findZeroHalfOrigSubSystem(A, consLB, consUB, varType, false, 1)
    (m,n) = size(A)
    model = Model(solver = CplexSolver(CPX_PARAM_SCRIND=0))
    @variable(model, varLB[j] <= x[j=1:n] <= varUB[j] )

    @objective(model, Min, sum( obj[j] * x[j] for j=1:n) )
    @constraint(model, consLB .<= A*x .<= consUB )
    solve(model)

    println("Objective value: ", getobjectivevalue(model))
    #println("x = ", getvalue(x))
    xVal = getvalue(x)
    for i=1:length(xVal)
        if xVal[i] > 0.001
            println("x[$i] = $(xVal[i])")
        end
    end
    #writeLP(model,"test.lp")

    done = false
    iter = 1
    while !done && iter <= maxIter
        zhsParams = ZeroHalfChvatGom.ZeroHalfSepParams(strength,1)
        cuts = ZeroHalfChvatGom.ZeroHalfSepTranslateVars(zhOrigSubSys.ASub, zhOrigSubSys.bSub, xVal, 0.0001, zhsParams, zhOrigSubSys.varsSorted)
        for i = 1:length(cuts)
            #println("*** Cut $i***")
            cut = cuts[i]
            vars = cut[1]
            coeffs = cut[2]
            rhs = cut[3]
            if printCuts
                for j=1:length(vars)
                    preString = ""
                    if j > 1
                        preString = " +"
                    end
                    if coeffs[j] == 1.0
                        print("$preString x$(vars[j])")
                    elseif coeffs[j] == -1.0
                        print(" - x$(vars[j])")
                    else
                        print("$preString $(coeffs[j])*x$(vars[j])")
                    end
                end
                println(" <= $rhs")
            end
            @constraint(model, sum(coeffs[j] * x[vars[j]] for j=1:length(vars)) <= rhs )
        end
        if length(cuts) > 0
            solve(model)

            println("Objective value (iter $iter): ", getobjectivevalue(model))
            #println("x = ", getvalue(x))
            xVal = getvalue(x)
            iter += 1
        else
            done = true
        end
    end
end

ZeroHalfOnIP.test("LP_and_blocks/Misc/wt008_028_2_mip.lp")
#ZeroHalfOnIP.test("LP_and_blocks/Misc/DagIntFlow_A-n32.lp")

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export datepickerinput

"""
    datepickerinput(;kwargs...)

A DatePickerInput component.

Keyword arguments:
- `id` (String; required): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `disabled` (Bool; optional): Specify whether or not the input should be disabled
- `hideLabel` (Bool; optional): Specify if the label should be hidden
- `iconDescription` (String; optional): The description of the calendar icon.
- `invalid` (Bool; optional): Specify whether or not the input should be invalid
- `invalidText` (String; optional): Specify the text to be rendered when the input is invalid
- `labelText` (String; required): Provide the text that will be read by a screen reader when visiting this control
- `placeholder` (String; optional): Specify the placeholder text
- `size` (a value equal to: 'sm', 'xl'; optional): Specify the size of the Date Picker Input. Currently supports either sm or xl as an option.
- `warn` (Bool; optional): Specify whether the control is currently in warning state
- `warnText` (String; optional): Provide the text that is displayed when the control is in warning state
- `value` (String; optional): Current value of the date
- `datePickerType` (a value equal to: 'range', 'single', 'simple'; optional): The type of the date picker
"""
function datepickerinput(; kwargs...)
        available_props = Symbol[:id, :disabled, :hideLabel, :iconDescription, :invalid, :invalidText, :labelText, :placeholder, :size, :warn, :warnText, :value, :datePickerType]
        wild_props = Symbol[]
        return Component("datepickerinput", "DatePickerInput", "dash_carbon_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components