text:
#   Unit "enlr.jl" of the PosDefManifoldML Package for Julia language

#   MIT License
#   Copyright (c) 2019,
#   Marco Congedo, CNRS, Grenoble, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements the elastic net logistic regression (ENLR)
#   machine learning classifier, including the lasso LR (default) and Ridge LR
#   as specific instances.

"""
```
abstract type ENLRmodel<:TSmodel end
```

Abstract type for **Elastic Net Logistic Rgression (ENLR)**
machine learning models. See [MLmodel](@ref).
"""
abstract type ENLRmodel<:TSmodel end

"""
```
mutable struct ENLR <: ENLRmodel
	metric      :: Metric = Fisher;
	alpha       :: Real = 1.0
	standardize :: Bool
	intercept   :: Bool
	meanISR     :: Union{ℍVector, Nothing}
	vecRange    :: UnitRange
	featDim     :: Int
	# GLMNet Models
	path        :: GLMNet.GLMNetPath
	cvλ         :: GLMNet.GLMNetCrossValidation
	best        :: GLMNet.GLMNetPath
end
```
ENLR machine learning models are incapsulated in this
mutable structure. Fields:

`.metric`, of type
[Metric](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#Metric::Enumerated-type-1),
is the metric that will be adopted to compute the mean used
as base-point for tangent space projection. By default the
Fisher metric is adopted. See [mdm.jl](@ref)
for the available metrics. If the data used to train the model
are not positive definite matrices, but Euclidean feature vectors,
the `.metric` field has no use. In order to use metrics you need to install the
*PosDefManifold* package.

`.alpha` is the hyperparameter in ``[0, 1]`` trading-off
the **elestic-net model**. ``α=0`` requests a pure **ridge** model and
``α=1`` a pure **lasso** model. By default, ``α=1`` is specified (lasso model).
This argument is usually passed as parameter to
the [`fit`](@ref) function, defaulting therein to ``α=1`` too.
See the examples here below.

All other fields do not correspond to arguments passed
upon creation of the model by the default creator.
Instead, they are filled later when a model is created by the
[`fit`](@ref) function:

For the content of fields `standardize`, `intercept`, `meanISR` and `vecRange`,
please see the documentation of the [`fit`](@ref) function.

if the data used to train the model are positive definite matrices,
`.featDim` is the length of the vectorized tangent vectors.
This is given by ``n(n+1)÷2`` (integer division), where ``n``
is the dimension of the original PD matrices on which the model is applied
once they are mapped onto the tangent space.
If feature vectors are used to train the model, `.featDim` is the length
of these vectors. If for fitting the model you have provided an optional
keyword argument `vecRange`, `.featDim` will be reduced accordingly.

`.path` is an instance of the following `GLMNetPath`
structure of the
[GLMNet.jl](https://github.com/JuliaStats/GLMNet.jl) package.
It holds the regularization path
that is created when the [`fit`](@ref) function is invoked
with optional keyword parameter `fitType` = `:path` or = `:all`:

```
struct GLMNetPath{F<:Distribution}
    family::F                        # Binomial()
    a0::Vector{Float64}              # intercept values for each solution
    betas::CompressedPredictorMatrix # coefficient values for each solution
    null_dev::Float64                # Null deviance of the model
    dev_ratio::Vector{Float64}       # R^2 values for each solution
    lambda::Vector{Float64}          # lambda values for each solution
    npasses::Int                     # actual number of passes over the
                                     # data for all lamda values
end
```

`.cvλ` is an instance of the following `GLMNetCrossValidation`
structure of the [GLMNet.jl](https://github.com/JuliaStats/GLMNet.jl) package.
It holds information about the
cross-validation used for estimating the optimal lambda
hyperparameter by the [`fit`](@ref) function when this is invoked
with optional keyword parameter `fitType` = `:best` (default) or = `:all`:

```
struct GLMNetCrossValidation
    path::GLMNetPath            # the cv path
    nfolds::Int                 # the number of folds for the cv
    lambda::Vector{Float64}     # lambda values for each solution
    meanloss::Vector{Float64}   # mean loss for each solution
    stdloss::Vector{Float64}    # standard deviation of the mean losses
end
```

`.best` is an instance of the `GLMNetPath`
structure (see above). It holds the model with the optimal
lambda parameter found by cross-validation
that is created by default
when the [`fit`](@ref) function is invoked.

**Examples**:
```
# Note: creating models with the default creator is possible,
# but not useful in general.

using PosDefManifoldML, PosDefManifold

# create an empty lasso model
m = ENLR(Fisher)

# since the Fisher metric is the default metric,
# this is equivalent to
m = ENLR()

# create an empty ridge model using the logEuclidean metric
m = ENLR(logEuclidean; alpha=0)

# Empty models can be passed as first argument of the `fit` function
# to fit a model. For instance, this will fit a ridge model of the same
# kind of `m` and put the fitted model in `m1`:
m1=fit(m, PTr, yTr)

# in general you don't need this machinery for fitting a model,
# since you can specify a model by creating one on the fly:
m2=fit(ENLR(logEuclidean; alpha=0), PTr, yTr)

# which is equivalent to
m2=fit(ENLR(logEuclidean), PTr, yTr; alpha=0)

# note that, albeit model `m` has been created as a ridge model,
# you have passed `m` and overwritten the `alpha` hyperparameter.
# The metric, instead, cannot be overwritten.

```
"""
mutable struct ENLR <: ENLRmodel
    metric   :: Metric
    alpha    :: Real
    standardize
    intercept
    meanISR
	vecRange
    featDim
    path
    cvλ
    best
    function ENLR(metric :: Metric=Fisher;
               alpha     :: Real = 1.0,
               standardize = nothing,
               intercept   = nothing,
               meanISR     = nothing,
			   vecRange    = nothing,
               featDim     = nothing,
               path        = nothing,
               cvλ         = nothing,
               best        = nothing)
        new(metric, alpha, standardize, intercept,
            meanISR, vecRange, featDim, path, cvλ, best)
    end
end

"""
```
function fit(model	:: ENLRmodel,
             𝐏Tr	 :: Union{ℍVector, Matrix{Float64}},
             yTr	:: IntVector;
	# parameters for projection onto the tangent space
	w		:: Union{Symbol, Tuple, Vector} = [],
	meanISR	:: Union{ℍ, Nothing} = nothing,
	meanInit:: Union{ℍ, Nothing} = nothing,
	vecRange:: UnitRange = 𝐏Tr isa ℍVector ? (1:size(𝐏Tr[1], 2)) : (1:size(𝐏Tr, 2)),
	fitType	:: Symbol = :best,
	verbose	:: Bool = true,
	⏩	   :: Bool = true,
	# arguments for `GLMNet.glmnet` function
	alpha			:: Real = model.alpha,
	weights			:: Vector{Float64} = ones(Float64, length(yTr)),
	intercept		:: Bool = true,
	standardize		:: Bool = true,
	penalty_factor	:: Vector{Float64} = ones(Float64, _getDim(𝐏Tr, vecRange)),
	constraints		:: Matrix{Float64} = [x for x in (-Inf, Inf), y in 1:_getDim(𝐏Tr, vecRange)],
	offsets			:: Union{Vector{Float64}, Nothing} = nothing,
	dfmax			:: Int = _getDim(𝐏Tr, vecRange),
	pmax			:: Int = min(dfmax*2+20, _getDim(𝐏Tr, vecRange)),
	nlambda			:: Int = 100,
	lambda_min_ratio:: Real = (length(yTr) < _getDim(𝐏Tr, vecRange) ? 1e-2 : 1e-4),
	lambda			:: Vector{Float64} = Float64[],
	tol				:: Real = 1e-5,
	maxit			:: Int = 1000000,
	algorithm		:: Symbol = :newtonraphson,
	# selection method
	λSelMeth	:: Symbol = :sd1,
	# arguments for `GLMNet.glmnetcv` function
	nfolds		:: Int = min(10, div(size(yTr, 1), 3)),
	folds		:: Vector{Int} =
	begin
		n, r = divrem(size(yTr, 1), nfolds)
		shuffle!([repeat(1:nfolds, outer=n); 1:r])
	end,
	parallel 	:: Bool=true)
```

Create and fit an [`ENLR`](@ref) machine learning model,
with training data `𝐏Tr`, of type
[ℍVector](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#%E2%84%8DVector-type-1),
and corresponding labels `yTr`, of type [IntVector](@ref).
Return the fitted model(s) as an instance of the [`ENLR`](@ref) structure.

As for all ML models acting in the tangent space,
fitting an ENLR model involves computing a mean of all the
matrices in `𝐏Tr`, mapping all matrices onto the tangent space
after parallel transporting them at the identity matrix
and vectorizing them using the
[vecP](https://marco-congedo.github.io/PosDefManifold.jl/dev/riemannianGeometry/#PosDefManifold.vecP)
operation. Once this is done, the elastic net logistic regression is fitted.

The mean is computed according to the `.metric` field
of the `model`, with optional weights `w`.
The `.metric` field of the `model` is passed to the [`tsMap`](@ref) function.
By default the metric is the Fisher metric. See the examples
here below to see how to change metric.
See [mdm.jl](@ref) or check out directly the documentation
of [PosDefManifold.jl](https://marco-congedo.github.io/PosDefManifold.jl/dev/)
for the available metrics.

**Optional keyword arguments**

By default, uniform weights will be given to all observations
for computing the mean to pass in the tangent space.
This is equivalent to passing as argument `w=:uniform` (or `w=:u`).
You can also pass as argument:

- `w=:balanced` (or simply `w=:b`). If the two classes are unbalanced,
  the weights should be inversely proportional to the number of examples
  for each class, in such a way that each class contributes equally
  to the computation of the mean.
  This is equivalent of passing `w=tsWeights(yTr)`. See the
  [`tsWeights`](@ref) function for details.
- `w=v`, where `v` is a user defined vector of non-negative weights for the
  observations, thus, `v` must contain the same number of elements as `yTr`.
  For example, `w=[1.0, 1.0, 2.0, 2.0, ...., 1.0]`
- `w=t`, where `t` is a 2-tuple of real weights, one weight for each class,
  for example `w=(0.5, 1.5)`.
  This is equivalent to passing `w=tsWeights(yTr; classWeights=collect(t))`,
  see the [`tsWeights`](@ref) function for details.

If `meanISR` is passed as argument, the mean is not computed,
instead this matrix is the inverse square root (ISR) of the mean
used for projecting the matrices in the tangent space (see [`tsMap`](@ref)).
Passed or computed, the inverse square root (ISR) of the mean
will be written in the `.meanISR` field of the created [`ENLR`](@ref) structure.
If `meanISRis` is not provided and the `.metric` field of the `model`
is Fisher, logdet0 or Wasserstein, the tolerance of the iterative algorithm
used to compute the mean is set to the argument passed as `tol` (default 1e-5).
Also, in this case a particular initialization for those iterative algorithms
can be provided as an `Hermitian` matrix with argument `meanInit`.

This function also allows to fit a model passing as
training data `𝐏Tr` directly a matrix of feature vectors,
where each feature vector is a row of the matrix.
In this case the `metric` of the ENLR model and argument `meanISR` are not used.
Therefore, the `.meanISR` field of the created [`ENLR`](@ref) structure
will be set to `nothing`.

If a `UnitRange` is passed with optional keyword argument `vecRange`,
then if `𝐏Tr` is a vector of `Hermitian` matrices, the vectorization
of those matrices once they are projected onto the tangent space
concerns only the rows (or columns) given in the specified range,
else if `𝐏Tr` is a matrix with feature vectors arranged in its rows, then
only the columns of `𝐏Tr` given in the specified range will be used.

If `fitType` = `:best` (default), a cross-validation procedure is run
to find the best lambda hyperparameter for the given training data.
This finds a single model that is written into the `.best` field
of the [`ENLR`](@ref) structure that will be created.

If `fitType` = `:path`, the regularization path for several values of
the lambda hyperparameter is found for the given training data.
This creates several models, which are written into the
`.path` field of the [`ENLR`](@ref) structure that will be created,
none of which is optimal, in the cross-validation sense, for the
given training data.

If `fitType` = `:all`, both the above fits are performed and all fields
of the [`ENLR`](@ref) structure that will be created will be filled in.

If `verbose` is true (default), information is printed in the REPL.
This option is included to allow repeated calls to this function
without crowding the REPL.

The `⏩` argument (true by default) is passed to the [`tsMap`](@ref)
function for projecting the matrices in `𝐏Tr` onto the tangent space
and to the `GLMNet.glmnetcv` function to run inner cross-validation
to find the `best` model using multi-threading.

The remaining optional keyword arguments, are

- the arguments passed to the `GLMNet.glmnet` function for fitting the models.
  Those are always used.

- the `λSelMeth` argument and the arguments passed to the `GLMNet.glmnetcv`
  function for finding the best lambda hyperparamater by cross-validation.
  Those are used only if `fitType` = `:path` or = `:all`.

**Optional keyword arguments for fitting the model(s) using GLMNet**

`alpha`: the hyperparameter in ``[0, 1]`` to trade-off
an elestic-net model. ``α=0`` requests a pure *ridge* model and
``α=1`` a pure *lasso* model. This defaults to 1.0,
which specifies a lasso model, unless the input [`ENLR`](@ref) `model`
has another value in the `alpha` field, in which case this value
is used. If argument `alpha` is passed here, it will overwrite
the `alpha` field of the input `model`.

`weights`: a vector of weights for each matrix (or feature vectors)
of the same size as `yTr`.
It defaults to 1 for all matrices.

`intercept`: whether to fit an intercept term.
The intercept is always unpenalized. Defaults to true.

`standardize`: if true (default), GLMNet standardize the predictors
(presumably this amounts to transform to unit variance) so that they are
in the same units. This is a common choice for regularized
regression models.

`penalty_factor`: a vector of length ``n(n+1)/2``, where ``n``
is the dimension of the original PD matrices on which the model is applied,
of penalties for each predictor in the tangent vectors.
This defaults to all ones, which weights each predictor equally.
To specify that a predictor should be unpenalized,
set the corresponding entry to zero.

`constraints`: an ``[n(n+1)/2]`` x ``2`` matrix specifying lower bounds
(first column) and upper bounds (second column) on each predictor.
By default, this is [-Inf Inf] for each predictor (each element
of tangent vectors).

`offset`: see documentation of original GLMNet package [🎓](@ref).

`dfmax`: The maximum number of predictors in the largest model.

`pmax`: The maximum number of predictors in any model.

`nlambda`: The number of values of ``λ`` along the path to consider.

`lambda_min_ratio`: The smallest ``λ`` value to consider,
as a ratio of the value of ``λ`` that gives the null model
(*i.e.*, the model with only an intercept).
If the number of observations exceeds the number of variables,
this defaults to 0.0001, otherwise 0.01.

`lambda`: The ``λ`` values to consider for fitting.
By default, this is determined
from `nlambda` and `lambda_min_ratio`.

`tol`: Is the convergence criterion for both the computation
of a mean for projecting onto the tangent space
(if the metric requires an iterative algorithm)
and for the GLMNet fitting algorithm. Defaults to 1e-5.
In order to speed up computations, you may try to set a lower `tol`;
The convergence will be faster but more coarse,
with a possible drop of classification accuracy,
depending on the signal-to-noise ratio of the input features.

`maxit`: The maximum number of iterations of the cyclic coordinate
descent algorithm. If convergence is not achieved, a warning is returned.

`algorithm`: the algorithm used to find the regularization path.
Possible values are `:newtonraphson` (default) and
`:modifiednewtonraphson`.

For further informations on those arguments, refer to the
resources on the GLMNet package [🎓](@ref).

**Optional Keyword arguments for finding the best model by cv**

`λSelMeth` = `:sd1` (default), the best model is defined as the one
allowing the highest `cvλ.meanloss` within one standard deviation of the
minimum, otherwise it is defined as the one allowing the minimum `cvλ.meanloss`.
Note that in selecting a model, the model with only the intercept term,
if it exists, is ignored. See [`ENLRmodel`](@ref) for a description
of the `.cvλ` field of the model structure.

Arguments `nfolds` and `folds` are passed to the
`GLMNet.glmnetcv` function along with the `⏩` argument.
Please refer to the resources on GLMNet
for details [🎓](@ref).

**See**: [notation & nomenclature](@ref), [the ℍVector type](@ref).

**See also**: [`predict`](@ref), [`cvAcc`](@ref).

**Tutorial**: [Example using the ENLR model](@ref).

**Examples**
```
using PosDefManifoldML, PosDefManifold

# generate some data
PTr, PTe, yTr, yTe=gen2ClassData(10, 30, 40, 60, 80, 0.1)

# Fit an ENLR lasso model and find the best model by cross-validation:
m=fit(ENLR(), PTr, yTr)

# ... balancing the weights for tangent space mapping
m=fit(ENLR(), PTr, yTr; w=tsWeights(yTr))

# ... using the log-Eucidean metric for tangent space projection
m=fit(ENLR(logEuclidean), PTr, yTr)

# Fit an ENLR ridge model and find the best model by cv:
m=fit(ENLR(Fisher), PTr, yTr; alpha=0)

# Fit an ENLR elastic-net model (α=0.9) and find the best model by cv:
m=fit(ENLR(Fisher), PTr, yTr; alpha=0.9)

# Fit an ENLR lasso model and its regularization path:
m=fit(ENLR(), PTr, yTr; fitType=:path)

# Fit an ENLR lasso model, its regularization path
# and the best model found by cv:
m=fit(ENLR(), PTr, yTr; fitType=:all)

```

"""
function fit(model  :: ENLRmodel,
               𝐏Tr  :: Union{ℍVector, Matrix{Float64}},
               yTr  :: IntVector;
		   # parameters for projection onto the tangent space
           w        	:: Union{Symbol, Tuple, Vector} = [],
           meanISR  	:: Union{ℍ, Nothing} = nothing,
		   meanInit 	:: Union{ℍ, Nothing} = nothing,
           vecRange 	:: UnitRange = 𝐏Tr isa ℍVector ? (1:size(𝐏Tr[1], 2)) : (1:size(𝐏Tr, 2)),
           fitType  	:: Symbol = :best,
		   verbose  	:: Bool = true,
           ⏩      	   :: Bool = true,
           # arguments for `GLMNet.glmnet` function
           alpha            :: Real = model.alpha,
           weights          :: Vector{Float64} = ones(Float64, length(yTr)),
           intercept        :: Bool = true,
		   standardize  	:: Bool = true,
           penalty_factor   :: Vector{Float64} = ones(Float64, _getDim(𝐏Tr, vecRange)),
           constraints      :: Matrix{Float64} = [x for x in (-Inf, Inf), y in 1:_getDim(𝐏Tr, vecRange)],
           offsets          :: Union{Vector{Float64}, Nothing} = nothing,
           dfmax            :: Int = _getDim(𝐏Tr, vecRange),
           pmax             :: Int = min(dfmax*2+20, _getDim(𝐏Tr, vecRange)),
           nlambda          :: Int = 100,
           lambda_min_ratio :: Real = (length(yTr) < _getDim(𝐏Tr, vecRange) ? 1e-2 : 1e-4),
           lambda           :: Vector{Float64} = Float64[],
           tol              :: Real = 1e-5,
           maxit            :: Int = 1000000,
           algorithm        :: Symbol = :newtonraphson,
           # selection method
           λSelMeth :: Symbol = :sd1,
           # arguments for `GLMNet.glmnetcv` function
           nfolds   :: Int = min(10, div(size(yTr, 1), 3)),
           folds    :: Vector{Int} =
           begin
               n, r = divrem(size(yTr, 1), nfolds)
               shuffle!([repeat(1:nfolds, outer=n); 1:r])
           end,
           parallel ::Bool=true)

    ⌚=now() # get the time in ms
    ℳ=deepcopy(model) # output model

	# overwrite fields in `ℳ` if the user has passed them here as arguments,
	# otherwise use as arguments the values in the fields of `ℳ`, e.g., the default
	if alpha ≠ 1.0 ℳ.alpha = alpha else alpha = ℳ.alpha end

    # check w argument and get weights for input matrices
    (w=_getWeights(w, yTr, "fit ("*_modelStr(ℳ)*" model)")) == nothing && return

    # other checks
    𝐏Tr isa ℍVector ? nObs=length(𝐏Tr) : nObs=size(𝐏Tr, 1)
    !_check_fit(ℳ, nObs, length(yTr), length(w), length(weights), "ENLR") && return

	# project data onto the tangent space or just copy the features if 𝐏Tr is a matrix
	X=_getFeat_fit!(ℳ, 𝐏Tr, meanISR, meanInit, tol, w, vecRange, true, verbose, ⏩)

    # convert labels in GLMNet format
    y = convert(Matrix{Float64}, [(yTr.==1) (yTr.==2)])

    # write some fields in output model struct
    ℳ.standardize = standardize
    ℳ.intercept   = intercept
    ℳ.featDim     = size(X, 2)
	ℳ.vecRange    = vecRange

    # collect the argumenst for `glmnet` function excluding the `lambda` argument
    fitArgs_λ = (alpha            = alpha,
                 weights          = weights,
                 standardize      = standardize,
                 intercept        = intercept,
                 penalty_factor   = penalty_factor,
                 constraints      = constraints,
                 offsets          = offsets,
                 dfmax            = dfmax,
                 pmax             = pmax,
                 nlambda          = nlambda,
                 lambda_min_ratio = lambda_min_ratio,
                 tol              = tol,
                 maxit            = maxit,
                 algorithm        = algorithm)

    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
    if fitType ∈(:path, :all)
        # fit the regularization path
        verbose && println("Fitting "*_modelStr(ℳ)*" reg. path...")
        ℳ.path = glmnet(X, y, Binomial();
                         lambda = lambda,
                         fitArgs_λ...) # glmnet Args but lambda
    end
    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
    if fitType ∈(:best, :all)
        verbose && println("Fitting best "*_modelStr(ℳ)*" model...")
        ℳ.cvλ=glmnetcv(X, y, Binomial();
                        nfolds   = nfolds,
                        folds    = folds,
                        parallel = parallel,
                        lambda   = lambda,
                        fitArgs_λ...) # glmnet Args but lambda

        # Never consider the model with only the intercept (0 degrees of freedom)
        l, i=length(ℳ.cvλ.lambda), max(argmin(ℳ.cvλ.meanloss), 1+intercept)

        # if bestλsel==:sd1 select the highest model with mean loss withinh 1sd of the minimum
        # otherwise the model with the smallest mean loss.
        thr=ℳ.cvλ.meanloss[i]+ℳ.cvλ.stdloss[i]
        λSelMeth==:sd1 ? (while i<l ℳ.cvλ.meanloss[i+1]<=thr ? i+=1 : break end) : nothing

        # fit the best model (only for the optimal lambda)
        ℳ.best = glmnet(X, y, Binomial();
                         lambda  = [ℳ.cvλ.path.lambda[i]],
                         fitArgs_λ...) # glmnet Args but lambda
    end
    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

    verbose && println(defaultFont, "Done in ", now()-⌚,".")
    return ℳ
end

"""
```
function predict(model   :: ENLRmodel,
		𝐏Te		:: Union{ℍVector, Matrix{Float64}},
		what		:: Symbol = :labels,
		fitType		:: Symbol = :best,
		onWhich		:: Int    = Int(fitType==:best);
	transfer	:: Union{ℍ, Nothing} = nothing,
	verbose		:: Bool = true,
	⏩		:: Bool = true)
```

Given an [`ENLR`](@ref) `model` trained (fitted) on 2 classes
and a testing set of ``k`` positive definite matrices `𝐏Te` of type
[ℍVector](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#%E2%84%8DVector-type-1),

if `what` is `:labels` or `:l` (default), return
the predicted **class labels** for each matrix in `𝐏Te`,
as an [IntVector](@ref).
Those labels are '1' for class 1 and '2' for class 2;

if `what` is `:probabilities` or `:p`, return the predicted **probabilities**
for each matrix in `𝐏Te` to belong to each classe, as a ``k``-vector
of ``z`` vectors holding reals in ``[0, 1]`` (probabilities).
The 'probabilities' are obtained passing to a
[softmax function](https://en.wikipedia.org/wiki/Softmax_function)
the output of the ENLR model and zero;

if `what` is `:f` or `:functions`, return the **output function** of the model,
which is the raw output of the ENLR model.

If `fitType` = `:best` (default), the best model that has been found by
cross-validation is used for prediction.

If `fitType` = `:path`,

- if `onWhich` is a valid serial number for a model in the `model.path`,
then this model is used for prediction,

- if `onWhich` is zero, all model in the `model.path` will be used for
predictions, thus the output will be multiplied by the number of
models in `model.path`.

Argumet `onWhich` has no effect if `fitType` = `:best`.

!!! note "Nota Bene"
    By default, the [`fit`](@ref) function fits only the `best` model.
    If you want to use the `fitType` = `:path` option you need to invoke
    the fit function with optional keyword argument `fitType`=`:path` or
    `fitType`=`:all`. See the [`fit`](@ref) function for details.

Optional keyword argument `transfer` can be used to specify the principal
inverse square root (ISR) of a new mean to be used as base point for
projecting the matrices in `𝐏Te` onto the tangent space.
By default `transfer` is equal to nothing,
implying that the base point will be the mean used to fit the model.
Passing a new mean ISR allows the *adaptation* first described in
Barachant et *al.*(2013). Typically `transfer` is the ISR
of the mean of the matrices in `𝐏Te` or of a subset of them.
Notice that this actually performs *transfer learning* by parallel
transporting both the training and test data to the identity matrix
as defined in Zanini et *al.*(2018) and later taken up in
Rodrigues et *al.*(2019)[🎓](@ref).

If `verbose` is true (default), information is printed in the REPL.
This option is included to allow repeated calls to this function
without crowding the REPL.

If ⏩ = true (default) and `𝐏Te` is an ℍVector type, the projection onto the
tangent space is multi-threaded.

**See**: [notation & nomenclature](@ref), [the ℍVector type](@ref).

**See also**: [`fit`](@ref), [`cvAcc`](@ref), [`predictErr`](@ref).

**Examples**
```
using PosDefManifoldML, PosDefManifold

# generate some data
PTr, PTe, yTr, yTe=gen2ClassData(10, 30, 40, 60, 80)

# fit an ENLR lasso model and find the best model by cv
m=fit(ENLR(Fisher), PTr, yTr)

# predict labels from the best model
yPred=predict(m, PTe, :l)
# prediction error
predErr=predictErr(yTe, yPred)

# predict probabilities from the best model
predict(m, PTe, :p)

# output functions from the best model
predict(m, PTe, :f)

# fit a regularization path for an ENLR lasso model
m=fit(ENLR(Fisher), PTr, yTr; fitType=:path)

# predict labels using a specific model
yPred=predict(m, PTe, :l, :path, 10)

# predict labels for all models
yPred=predict(m, PTe, :l, :path, 0)
# prediction error for all models
predErr=[predictErr(yTe, yPred[:, i]) for i=1:size(yPred, 2)]

# predict probabilities from a specific model
predict(m, PTe, :p, :path, 12)

# predict probabilities from all models
predict(m, PTe, :p, :path, 0)

# output functions from specific model
predict(m, PTe, :f, :path, 3)

# output functions for all models
predict(m, PTe, :f, :path, 0)

```

"""
function predict(model   :: ENLRmodel,
                 𝐏Te     :: Union{ℍVector, Matrix{Float64}},
                 what    :: Symbol = :labels,
                 fitType :: Symbol = :best,
                 onWhich :: Int    = Int(fitType==:best);
			transfer   :: Union{ℍ, Nothing} = nothing,
            verbose    :: Bool = true,
            ⏩        :: Bool = true)

    ⌚=now()

    # checks
    if !_whatIsValid(what, "predict ("*_modelStr(model)*")") return end
    if !_fitTypeIsValid(fitType, "predict ("*_modelStr(model)*")") return end
    if fitType==:best && model.best==nothing @error 📌*", predict function: the best model has not been fitted; run the `fit`function with keyword argument `fitType=:best` or `fitType=:all`"; return end
    if fitType==:path && model.path==nothing @error 📌*", predict function: the regularization path has not been fitted; run the `fit`function with keyword argument `fitType=:path` or `fitType=:all`"; return end
    if !_ENLRonWhichIsValid(model, fitType, onWhich, "predict ("*_modelStr(model)*")") return end

    # projection onto the tangent space
	X=_getFeat_Predict!(model, 𝐏Te, transfer, model.vecRange, true, verbose, ⏩)

    # prediction
    verbose && println("Predicting "*_ENLRonWhichStr(model, fitType, onWhich)*"...")
    if 		fitType==:best
        	path=model.best
        	onWhich=1
    elseif  fitType==:path
        	path=model.path
    end

    onWhich==0 ? π=GLMNet.predict(path, X) : π=GLMNet.predict(path, X, onWhich)

    k, l=size(π, 1), length(path.lambda)
    if     	what == :functions || what == :f
        	🃏=π
    elseif 	what == :labels || what == :l
        	onWhich==0 ? 🃏=[π[i, j]<0 ? 1 : 2 for i=1:k, j=1:l] : 🃏=[y<0 ? 1 : 2 for y ∈ π]
    elseif 	what == :probabilities || what == :p
        	onWhich==0 ? 🃏=[softmax([-π[i, j], 0]) for i=1:k, j=1:l] : 🃏=[softmax([-y, 0]) for y ∈ π]
    end

    verbose && println(defaultFont, "Done in ", now()-⌚,".")
    verbose && println(titleFont, "\nPredicted ",_what2Str(what),":", defaultFont)
    return 🃏
end

# ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, M::ENLR)
    if M.path==nothing
        if M.best==nothing
            println(io, greyFont, "\n↯ ENLR GLMNet machine learning model")
            println(io, "⭒  ⭒    ⭒       ⭒          ⭒")
            println(io, ".metric : ", string(M.metric))
            println(io, ".alpha  : ", "$(round(M.alpha, digits=3))", defaultFont)
            println(io, "Unfitted model")
            return
        end
    end

    println(io, titleFont, "\n↯ GLMNet ENLR machine learning model")
    println(io, separatorFont, "  ", _modelStr(M))
    println(io, separatorFont, "⭒  ⭒    ⭒       ⭒          ⭒", defaultFont)
    println(io, "type    : PD Tangent Space model")
    println(io, "features: tangent vectors of length $(M.featDim)")
    println(io, "classes : 2")
    println(io, separatorFont, "Fields  : ")
	# # #
	println(io, greyFont, " Tangent Space Parametrization", defaultFont)
    println(io, separatorFont," .metric      ", defaultFont, string(M.metric))
	if M.meanISR == nothing
        println(io, greyFont, " .meanISR      not created")
    else
        n=size(M.meanISR, 1)
        println(io, separatorFont," .meanISR     ", defaultFont, "$(n)x$(n) Hermitian matrix")
    end
    println(io, separatorFont," .vecRange    ", defaultFont, "$(M.vecRange)")
    println(io, separatorFont," .featDim     ", defaultFont, "$(M.featDim)")
	# # #
	println(io, greyFont, " ENLR Parametrization", defaultFont)
    println(io, separatorFont," .alpha       ", defaultFont, "$(round(M.alpha, digits=3))")
    println(io, separatorFont," .intercept   ", defaultFont, string(M.intercept))
	println(io, separatorFont," .standardize ", defaultFont, string(M.standardize))

    if M.path==nothing
        println(io, greyFont," .path struct `GLMNetPath`:")
        println(io, "       not created ")
    else
        println(io, separatorFont," .path", defaultFont," struct `GLMNetPath`:")
        println(io, titleFont,"       .family, .a0, .betas, .null_dev, ")
        println(io, titleFont,"       .dev_ratio, .lambda, .npasses")
    end
    if M.cvλ==nothing
        println(io, greyFont," .cvλ  struct `GLMNetCrossValidation`:")
        println(io, "       not created ")
    else
        println(io, separatorFont," .cvλ", defaultFont,"  struct `GLMNetCrossValidation`:")
        println(io, titleFont,"       .path, .nfolds, .lambda, ")
        println(io, titleFont,"       .meanloss, stdloss")
    end
    if M.best==nothing
        println(io, greyFont," .best struct `GLMNetPath`:")
        println(io, "       not created ")
    else
        println(io, separatorFont," .best", defaultFont," struct `GLMNetPath`:")
        println(io, titleFont,"       .family, .a0, .betas, .null_dev, ")
        println(io, titleFont,"       .dev_ratio, .lambda, .npasses")
    end
end


text:
basis_blades(sig, grade=1) = [Blade{sig}(1, i) for i ∈ bits_of_grade(grade, dimension(sig))]

basis(sig; grade=1) = basis_blades(sig, grade)
basis(dim::Integer; grade=1) = basis(Tuple(fill(1, dim)); grade=1)
basis(osig::OffsetSignature{sig,indices}; grade=1) where {sig,indices} = OffsetArray(basis_blades(osig, grade), indices)

# return expression containing variable assignments for each basis blade for the given signature.
function generate_blades(combos, sig)
	bvs = collect(basis(sig)) # needed because `OffsetArrays` don't work with, e.g., `powerset`
	vars = Symbol[]
	vals = Blade{sig}[]
	for (ordered_bvs, ublade) ∈ zip(combos(bvs), combos(1:dimension(sig)))
		varname = basis_blade_label(sig, ublade)
		isempty(varname) && continue
		push!(vars, Symbol(varname))
		push!(vals, prod(ordered_bvs))
	end
	quote
		@info "Defined basis blades $(join($vars, ", "))"
		$([ :($(esc(k)) = $(v)) for (k, v) ∈ zip(vars, vals) ]...)
		nothing
	end
end
generate_blades(sig) = generate_blades(x -> [[i] for i ∈ x], sig)

parse_sig(n::Integer) = Tuple(ones(n))
parse_sig(sig::MetricSignature) = sig
parse_sig(sig::Tuple{Vararg{<:Integer}}) = sig
parse_sig(sig::NamedTuple) = sig
parse_sig(labels::Tuple{Vararg{Symbol}}) = NamedTuple{labels}(ones(Int, length(labels)))
parse_sig(sig::AbstractString) = Tuple(get(Dict('+'=>1, '-'=>-1, '0'=>0), i) do i
	throw(ArgumentError("Invalid character in signature string: $i"))
end for i ∈ sig)
parse_sig(sig) = throw(ArgumentError("could not interpret $sig as metric signature"))

function parse_sig(sig::Tuple)
	# signature specified by series of symbols/exprs `label[=square]`
	labels = Symbol[]
	squares = []
	for bv ∈ sig
		if bv isa Symbol
			label, sq = bv, 1
		elseif bv isa Expr && bv.head == :(=)
			label, sq = bv.args
		else
			throw(ArgumentError("Invalid term '$bv' in signature specification"))
		end
		push!(labels, label)
		push!(squares, sq)
	end
	NamedTuple{tuple(labels...)}(squares)
end

# for a single argument `@basis X`, interpret `X` as an expression which evaluates to a signature.
# for many arguments, `@basis x1 x2 x3...`, interpret `(x1, x2, x3)` as a signature.
parse_sig(::Module, args...) = parse_sig(args)
parse_sig(m::Module, args) = let a = @eval m $args
	parse_sig(a)
end

"""
	@basis sig
	@basis x y z...

Populate the local namespace with basis blades of all grades for
the geometric algebra with signature `sig`, which may be provided
as a string, tuple, named tuple or `<:AbstractMetricSignature`.
If multiple arguments are given, each is interpreted as specifying
the label and square (defaulting to `1`) of a basis vector.
For an `n`-dimensional signature, `2^n` basis blades are defined.

See also `@basisall`.

Examples
===
```jldoctest
julia> @basis "+++"
[ Info: Defined basis blades v, v1, v2, v3, v12, v13, v23, v123

julia> @basis (+1,0,-1)
[ Info: Defined basis blades v, v1, v2, v3, v12, v13, v23, v123

julia> @basis x y
[ Info: Defined basis blades x, y, xy

julia> @basis dt=-1 dx=1 dy
[ Info: Defined basis blades dt, dx, dy, dtdx, dtdy, dxdy, dtdxdy

julia> dt^2 + dxdy*dy
MixedMultivector{⟨dt-,dx+,dy+⟩, Vector{Int64}}:
 -1
 1 dx
```
"""
macro basis(a, args...)
	generate_blades(powerset, parse_sig(__module__, a, args...))
end

"""
	@basisall

Similarly to `@basis`, populate the local namespace with basis blades,
but including all permutations of the basis blades.

Note that more than `2^n` variables are defined for a signature of dimension `n`!

Example
===
```jldoctest
julia> @basisall x y z
[ Info: Defined basis blades x, y, z, xy, yx, xz, zx, yz, zy, xyz, xzy, yxz, yzx, zxy, zyx

julia> zyx
Grade-3 Blade{⟨x+,y+,z+⟩, 3, 0b111, Int64}:
 -1 xyz
```
"""
macro basisall(a, args...)
	generate_blades(parse_sig(__module__, a, args...)) do bvs
		Iterators.flatten(permutations.(powerset(bvs)))
	end
end

ext: jl


text:
"""
    seek_connected(f, N::Int, head::DataType, tail::Tuple) -> (N, remainder)

Recursively scan a tuple of `DataType` (split into its `head` and
`tail`) to compute the length `N` of the uninterrupted sequence
of adjacent operations where the predicate `f` is true.
Additionally the `remainder` of the tuple (without that sequence)
is also returned.
"""
@inline function seek_connected(f, N::Int, head::Type{<:Operation}, tail::Tuple)
    if f(head)
        seek_connected(f, N+1, first(tail), Base.tail(tail))
    else
        N, (head, tail...)
    end
end

@inline function seek_connected(f, N::Int, head::Type{<:Operation}, tail::Tuple{})
    if f(head)
        N+1, ()
    else
        N, (head,)
    end
end

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

@inline supports_lazy(head::Type{<:Operation}, tail::Tuple{}) = false
@inline supports_lazy(head::Type{<:Operation}, tail::Tuple) =
    supports_lazy(head) && supports_lazy(first(tail))
@inline uses_affinemap(head::Type{<:Operation}, tail::Tuple{}) = false
@inline uses_affinemap(head::Type{<:Operation}, tail::Tuple) =
    uses_affinemap(head) && uses_affinemap(first(tail))

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

function augment_impl(var_offset::Int, op_offset::Int, pipeline::Tuple, args...)
    augment_impl(var_offset, op_offset, first(pipeline), Base.tail(pipeline), args...)
end

function augment_impl(var_offset::Int, op_offset::Int, pipeline::Tuple{}, args...)
    :($(Symbol(:img_, var_offset)))
end

function augment_impl(var_offset::Int, op_offset::Int, head::DataType, tail::NTuple{N,DataType}, avoid_eager = false) where N
    var_in  = Symbol(:img_, var_offset)
    var_out = Symbol(:img_, var_offset+1)
    if supports_lazy(head, tail)
        # If reached there are at least two adjacent lazy operations
        num_affine, after_affine = uses_affinemap(head, tail) ? seek_connected(uses_affinemap, 0, head, tail) : (0, nothing)
        num_special, _ = seek_connected(x->(supports_permute(x)||supports_view(x)||supports_stepview(x)), 0, head, tail)
        num_lazy, after_lazy = seek_connected(supports_lazy, 0, head, tail)
        if num_special >= num_affine
            quote
                $var_out = unroll_applylazy($(Expr(:tuple, (:(pipeline[$i]) for i in op_offset:op_offset+num_lazy-1)...)), $var_in)
                $(augment_impl(var_offset+1, op_offset+num_lazy, after_lazy, avoid_eager))
            end
        else
            quote
                $var_out = unroll_applyaffine($(Expr(:tuple, (:(pipeline[$i]) for i in op_offset:op_offset+num_affine-1)...)), $var_in)
                $(augment_impl(var_offset+1, op_offset+num_affine, after_affine, avoid_eager))
            end
        end
    else
        # At most "head" is lazy (i.e. tail[1] is surely not).
        # Unless "avoid_eager==true" we prefer using "applyeager" in
        # this case because there is no neighbour synergy and we
        # assume "applyeager" is more efficient.
        if !supports_lazy(head) || (!avoid_eager && (length(tail) == 0 || supports_eager(head)))
            quote
                $var_out = applyeager(pipeline[$op_offset], $var_in)
                $(augment_impl(var_offset+1, op_offset+1, tail, avoid_eager))
            end
        else
            quote
                $var_out = applylazy(pipeline[$op_offset], $var_in)
                $(augment_impl(var_offset+1, op_offset+1, tail, avoid_eager))
            end
        end
    end
end

function augment_impl(varname::Symbol, pipeline::NTuple{N,DataType}, args...) where N
    quote
        img_1 = $varname
        $(augment_impl(1, 1, pipeline, args...))
    end
end

# --------------------------------------------------------------------
# The following two methods are just for user inspection and
# debugging purposes as they show what kind of code the pipeline
# generates. They are not called internally by the package itself.
#
# Example:
#
#   julia> AugmentorBase.augment_impl(Rotate(45) |> Scale(0.9) |> CacheIntermediate() |> FlipX() |> FlipY())
#   quote
#       img_1 = input_image
#       begin
#           img_2 = unroll_applyaffine((pipeline[1], pipeline[2]), img_1)
#           begin
#               img_3 = applyeager(pipeline[3], img_2)
#               begin
#                   img_4 = unroll_applylazy((pipeline[4], pipeline[5]), img_3)
#                   img_4
#               end
#           end
#       end
#   end

function augment_impl(pipeline::AbstractPipeline; avoid_eager = false)
    augment_impl(:input_image, map(typeof, operations(pipeline)), avoid_eager)
end

augment_impl(op::Operation; kw...) = augment_impl((op,); kw...)


text:
using Gurobi
using JuMP
using YAML
using Configurations

include("make_instance.jl")
include("load_instance.jl")

model = Model(Gurobi.Optimizer)
set_optimizer_attribute(model, "threads", 8)

@option struct Configuration
    instance_file::String
    instance_name::String
end

function load_config(file_name)
    data = YAML.load_file(file_name, dicttype=Dict{String, Any})
    from_dict(Configuration, data)
end

struct ModelVariables
    x
    y
    z
end

function build_model!(model, A, Γ, c, Φ, Π, Σ, κ, t, Λ)
    n, m = size(A)
    @variable(model, 0 <= x[1:n], Int)
    @variable(model, z[1:n], Bin)
    @variable(model, y[i = 1:n, j = 1:m], Bin)

    @constraint(model, c0, Γ' * z + (c .+ Φ)' * x <= t)
    @constraint(model, c1[i = 1:n],  x[i] <= Σ[i] * z[i])
    @constraint(model, c2, x .>= Π)
    @constraint(model, c3[i  = 1:n], κ * x[i] <=  Λ[i])
    @constraint(model, c4[i = 1:n], A[i, :]' * y[i, :] <= κ * x[i])

    @objective(model, Max, A[:]' * y[:])
    ModelVariables(x, y, z)
end

config = load_config("config.yml")

instance = load_istance(config.instance_file, config.instance_name)

A, Λ = make_instance(instance.clients, instance.stations)
vars = build_model!(model, A, instance.Γ, instance.c, instance.ϕ, instance.Π, instance.Σ, instance.k, instance.T, Λ)
optimize!(model)
if termination_status(model) == OPTIMAL
    println(value.(vars.x))
    println(value.(vars.z))
end


text:
function FieldHermiteInterpolator()
    return FieldHermiteInterpolator(())
end

function add_sample_point(obj::FieldHermiteInterpolator, arg0::FieldElement, arg1::Vector{Vector{FieldElement}})
    return jcall(obj, "addSamplePoint", void, (FieldElement, Vector{Vector{FieldElement}}), arg0, arg1)
end

function derivatives(obj::FieldHermiteInterpolator, arg0::FieldElement, arg1::jint)
    return jcall(obj, "derivatives", Vector{Vector{FieldElement}}, (FieldElement, jint), arg0, arg1)
end

function value(obj::FieldHermiteInterpolator, arg0::FieldElement)
    return jcall(obj, "value", Vector{FieldElement}, (FieldElement,), arg0)
end


text:
using NamedDims
using NamedDims: names
using Test
using Statistics

@testset "sum" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test sum(nda) == 100
    @test sum(nda; dims=:x) == [40 60] == sum(nda; dims=1)

    @test names(sum(nda; dims=:x)) == (:x, :y) == names(sum(nda; dims=1))
end

@testset "mean" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test mean(nda) == 25
    @test mean(nda; dims=:x) == [20 30] == mean(nda; dims=1)

    @test names(mean(nda; dims=:x)) == (:x, :y) == names(mean(nda; dims=1))
end

@testset "mapslices" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test mapslices(join, nda; dims=:x) == ["1030" "2040"] == mapslices(join, nda; dims=1)
    @test mapslices(join, nda; dims=:y) == reshape(["1020", "3040"], Val(2)) == mapslices(join, nda; dims=2)

    @test names(mapslices(join, nda; dims=:y)) == (:x, :y) == names(mapslices(join, nda; dims=2))
end

@testset "mapreduce" begin
    nda = NamedDimsArray([10 20; 31 40], (:x, :y))
    @test mapreduce(isodd, |, nda; dims=:x) == [true false] == mapreduce(isodd, |, nda; dims=1)
    @test mapreduce(isodd, |, nda; dims=:y) == [false true]' == mapreduce(isodd, |, nda; dims=2)

    @test names(mapreduce(isodd, |, nda; dims=:y)) == (:x, :y) == names(mapreduce(isodd, |, nda; dims=2))
end

@testset "zero" begin
    nda = NamedDimsArray([10 20; 31 40], (:x, :y))
    @test zero(nda) == [0 0; 0 0]
    @test names(zero(nda)) == (:x, :y)
end


text:
using Test
using OpenQuantumSystems
using SparseArrays, LinearAlgebra

@testset "superoperators" begin

    b = GenericBasis([2])
    A = DenseOperator(b, b, [1 2; 3 4])
    B = DenseOperator(b, b, [5 6; 7 8])
    C = Commutator(A)
    C_ref = DenseSuperOperator(
        (b, b),
        (b, b),
        [0 2 -3 0; 3 3 0 -3; -2 0 -3 2; 0 -2 3 0],
    )
    @test Commutator(A) == C_ref
    @test Commutator(A) * B == A * B - B * A

end # testset

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems


text:
# test FluxModelsPose2Pose2

using LinearAlgebra
using Flux
using RoME

import RoME: FluxModelsPose2Pose2

##

include(joinpath((@__DIR__), "LoadPyNNTxt.jl"))

allModels = []
for i in 0:99
# /home/dehann/data/racecar/results/conductor/models/retrained_network_weights0
  push!(allModels, loadTfModelIntoFlux(ENV["HOME"]*"/data/racecar/results/conductor/models/retrained_network_weights$i") )
end

##

# start with a basic factor graph
fg = generateCanonicalFG_ZeroPose(varType=Pose2)

mvnNaive = MvNormal(zeros(3), diagm([1.0;1.0;0.01]))

addVariable!(fg, :x1, Pose2)

jvd = zeros(25,4)
pp = FluxModelsPose2Pose2(allModels, jvd, mvnNaive, 0.5)

##

addFactor!(fg, [:x0;:x1], pp)

pts = approxConv(fg, :x0x1f1, :x1)

#

##  Special async task to add Neural odo to fg when data becomes available.

# # latest pose
# allVars = ls(fec.slam.dfg, r"x\d") |> sortDFG
# # does it have FluxModelsPose2Pose2 factor?
# allFluxFct = ls(fec.slam.dfg, FluxModelsPose2Pose2)
# mask = (x->intersect(ls(fec.slam.dfg, x), allFluxFct ) |> length).(allVars) .== 0
# varsWithout = allVars[mask][end-5:end]
#
# # what are the command values from previous pose
# let prevPs = prevPs, i=i
# prevPs = varsWithout[1]
# i = 2
# for i in 2:length(varsWithout)
#
# ps = varsWithout[i]
# theVar = getVariable(fec.slam.dfg, ps)
# # skip if no entry yet
# # hasDataEntry(theVar, :JOYSTICK_CMD_VALS) ? nothing : continue
# cmdData = fetchDataElement(theVar, fec.datastore, :JOYSTICK_CMD_VALS)
# throttle = (x->x[3].axis[2]).(cmdData)
# steering = (x->x[3].axis[4]).(cmdData)
# DT = ( getTimestamp(theVar) - getTimestamp(getVariable(fec.slam.dfg, prevPs)) ).value*1e-3
# xj = getPPE(fec.slam.dfg, ps).suggested
# xi = getPPE(fec.slam.dfg, prevPs).suggested
# biRw = TU.R(-xi[3])
# biV = biRw * (xj[1:2] - xi[1:2]) / DT
# biVrep = repeat(biV', length(throttle))
# joyval = hcat(throttle, steering, biVrep)
#
# #
# @show prevPs, ps, size(joyval,1)
#
# joyval
# itpJoy = interpToOutx4(joyval)
#
# JlOdoPredictorPoint2(itpJoy,allModels)
#
#
# prevPs = ps
# end
# end
#
# # interpolate joyvals to right size
# # joyval
#
# # the naive model (should be the camera)
# DXmvn = MvNormal(zeros(3),diagm([0.4;0.1;0.4].^2))
#
# fmp2 = FluxModelsPose2Pose2(x->JlOdoPredictorPoint2(x,allModels), joyval25, DXmvn,0.6)
#
#
# # convert to [25 x 4] input


text:
##
# MIME stuff
##

export as_mime

struct AsMIME{M <: MIME}
    mime::M
    x::Any
end

Base.show(io::IO, m::M, am::AsMIME{M}) where M <: MIME = Base.show(io, m, am.x)

function as_mime(m::MIME)
    x -> as_mime(m, x)
end

function as_mime(m::MIME, x)
    AsMIME(m, x)
end

export as_text, as_html, as_svg, as_png
const as_text = as_mime(MIME"text/plain"())
const as_html = as_mime(MIME"text/html"())
const as_svg = as_mime(MIME"image/svg+xml"())
const as_png = as_mime(MIME"image/png"())
# as_jpg = as_mime(MIME"image/jpg"())
# as_bmp = as_mime(MIME"image/bmp"())
# as_gif = as_mime(MIME"image/gif"())

##
# terminal classics
##

export Dump, Show, Print

"""
    Dump(x; maxdepth=8)

Every part of the representation of a value. The depth of the output is truncated at maxdepth.

This is a variant of [`Base.dump`](@ref) that returns the representation directly, instead of printing it to stdout.

See also: [`Print`](@ref) and [`with_terminal`](@ref).
"""
struct Dump
    x::Any
    maxdepth::Integer
end

Dump(x; maxdepth=8) = Dump(x, maxdepth)

function Base.show(io::IO, ::MIME"text/plain", d::Dump)
    dump(io, d.x; maxdepth=d.maxdepth)
end

"""
    Show(mime::MIME, data)

An object that can be rendered as the `mime` MIME type, by writing `data` to the IO stream. For use in environments with rich output support. Read more about [`Base.show`](@ref).

# Examples

```julia
Show(MIME"text/html"(), "I can be <b>rendered</b> as <em>HTML</em>!")
```

```julia
Show(MIME"image/png"(), read("dog.png"))
```

`Base.showable` and `Base.show` are defined for a `Show`.

```julia
s = Show(MIME"text/latex"(), "\\\\frac{hello}{world}")

showable(MIME"text/latex"(), s) == true

repr(MIME"text/latex"(), s) == "\\\\frac{hello}{world}"
```

"""
struct Show{M <: MIME}
    mime::M
    data
end

Base.show(io::IO, ::M, x::Show{M}) where M <: MIME = write(io, x.data)

"""
    Print(xs...)

The text that would be printed when calling `print(xs...)`. Use `string(xs...)` if you want to use the result as a `String`.

See also: [`Dump`](@ref) and [`with_terminal`](@ref).
"""
Print(xs...) = Text(string(xs...))

##
# IO context
##

export WithIOContext

"""
    WithIOContext(x, properties::Pair...)

A wrapper around `x` with extra IOContext properties set, just for the display of `x`.

# Examples

```julia
WithIOContext(rand(100,100), :compact => false)
```

```julia
large_df = DataFrame(rand(100,100))
WithIOContext(large_df, :displaysize => (9999,9999))
```
"""
struct WithIOContext{T}
    x::T
    context::IOContext{Base.DevNull}
    WithIOContext(x::T, props::Pair...; moreprops...) where T = new{T}(x, IOContext(devnull, props..., moreprops...))
end

Base.:(==)(w1::WithIOContext, w2::WithIOContext) = w1.x == w2.x && w1.context == w2.context

function Base.show(io::IO, m::MIME, w::WithIOContext)
	Base.show(IOContext(io, w.context), m, w.x)
end
# we need to add a more specific method for text/plain because Julia base has a fallback method here
function Base.show(io::IO, m::MIME"text/plain", w::WithIOContext)
	Base.show(IOContext(io, w.context), m, w.x)
end

Base.showable(m::MIME, w::WithIOContext) = Base.showable(m, w.x)


text:
using DropboxSDK

filename(entry) =
    entry.path_display === nothing ? entry.name : entry.path_display

@testset "Get authorization" begin
    global auth = get_authorization()
    @test auth isa Authorization
end

@testset "Calculate content hash" begin
    content = read(joinpath("..", "data", "ssc2006-02a1_Lrg.jpg"))
    @assert length(content) == 6056980
    content_hash = calc_content_hash(content)
    @test (content_hash ==
           "8a6ebea6983dc68be1575676d3a8ec0d664cfee69b2dbcdf44087cf5d455fe12")

    len = length(content)
    di = 1234567
    contents = []
    for i in 1:di:len
        push!(contents, @view content[i : min(len, i + di - 1)])
    end
    @assert sum(length.(contents)) == length(content)
    # @assert vcat(contents...) == content

    data_channel = Channel{AbstractVector{UInt8}}(0)
    content_hash_task =
        start_task(() -> calc_content_hash(data_channel),
                   ("Calculate content hash", :calc_content_hash))
    for content in contents
        put!(data_channel, content)
    end
    close(data_channel)
    content_hash = fetch(content_hash_task)
    @test (content_hash ==
           "8a6ebea6983dc68be1575676d3a8ec0d664cfee69b2dbcdf44087cf5d455fe12")
end

@testset "Get current account" begin
    account = users_get_current_account(auth)
    first = account.name.given_name
    last = account.name.surname
    display = account.name.display_name
    @test first == "Erik"
    @test last == "Schnetter"
    @test display == "Erik Schnetter (PI)"
end

@testset "Get space usage" begin
    usage = users_get_space_usage(auth)
    used = usage.used
    @test used isa Integer
    @test used >= 0
end

@testset "List folder" begin
    entries = files_list_folder(auth, "", recursive=false)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == 0
end

@testset "Create folder" begin
    files_create_folder(auth, "/$folder")
    entries = files_list_folder(auth, "", recursive=false)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == 1
    @test count(entry -> entry.path_display == "/$folder", entries) == 1
end

@testset "Upload file" begin
    metadata =
        files_upload(auth, "/$folder/file", Vector{UInt8}("Hello, World!\n"))
    @test metadata isa FileMetadata
    @test metadata.size == length("Hello, World!\n")
    entries = files_list_folder(auth, "/$folder", recursive=true)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == 2
    @test count(entry -> entry.path_display == "/$folder", entries) == 1
    @test count(entry -> entry.path_display == "/$folder/file", entries) == 1
end

@testset "Download file" begin
    metadata, content = files_download(auth, "/$folder/file")
    @test metadata.path_display == "/$folder/file"
    @test String(content) == "Hello, World!\n"
end

@testset "Get file metadata" begin
    metadata = files_get_metadata(auth, "/$folder/file")
    @test metadata.path_display == "/$folder/file"
    @test (metadata.content_hash ==
           calc_content_hash(Vector{UInt8}("Hello, World!\n")))
end

@testset "Upload empty file" begin
    metadata = files_upload(auth, "/$folder/file0", Vector{UInt8}(""))
    @test metadata isa FileMetadata
    @test metadata.size == 0
    entries = files_list_folder(auth, "/$folder", recursive=true)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == 3
    @test count(entry -> entry.path_display == "/$folder", entries) == 1
    @test count(entry -> entry.path_display == "/$folder/file0", entries) == 1
end

@testset "Download empty file" begin
    metadata, content = files_download(auth, "/$folder/file0")
    @test metadata.path_display == "/$folder/file0"
    @test String(content) == ""
end

@testset "Upload file in chunks" begin
    upload_channel = Channel{Vector{UInt8}}(0)
    upload_task =
        start_task(() -> files_upload(auth, "/$folder/file1", upload_channel),
                   ("Upload file in chunks", :files_upload, "/$folder/file1"))
    put!(upload_channel, Vector{UInt8}("Hello, "))
    put!(upload_channel, Vector{UInt8}("World!\n"))
    close(upload_channel)
    metadata = fetch(upload_task)
    @test metadata isa FileMetadata
    @test metadata.size == length("Hello, World!\n")
    entries = files_list_folder(auth, "/$folder", recursive=true)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == 4
    @test count(entry -> entry.path_display == "/$folder", entries) == 1
    @test count(entry -> entry.path_display == "/$folder/file1", entries) == 1
end

@testset "Download file" begin
    metadata, content = files_download(auth, "/$folder/file1")
    @test metadata.path_display == "/$folder/file1"
    @test String(content) == "Hello, World!\n"
end

@testset "Upload empty file in chunks" begin
    upload_channel = Channel{Vector{UInt8}}(0)
    upload_task =
        start_task(() -> files_upload(auth, "/$folder/file2", upload_channel),
                   ("Upload empty file in chunks", :files_upload,
                    "/$folder/file2"))
    close(upload_channel)
    metadata = fetch(upload_task)
    @test metadata isa FileMetadata
    @test metadata.size == 0
    entries = files_list_folder(auth, "/$folder", recursive=true)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == 5
    @test count(entry -> entry.path_display == "/$folder", entries) == 1
    @test count(entry -> entry.path_display == "/$folder/file2", entries) == 1
end

@testset "Download empty file" begin
    metadata, content = files_download(auth, "/$folder/file2")
    @test metadata.path_display == "/$folder/file2"
    @test String(content) == ""
end

const numfiles = 4
@testset "Upload several files" begin
    upload_spec_channel = Channel{UploadSpec}(0)
    upload_task = start_task(() -> files_upload(auth, upload_spec_channel),
                             ("Upload several files", :files_upload))
    for i in 0:numfiles-1
        data_channel = Channel{Vector{UInt8}}(0)
        put!(upload_spec_channel,
             UploadSpec("/$folder/files$i", data_channel))
        for j in 1:i
            put!(data_channel, Vector{UInt8}("Hello, World!\n"))
        end
        close(data_channel)
    end
    close(upload_spec_channel)
    metadatas = fetch(upload_task)
    @test length(metadatas) == numfiles
    @test all(metadata isa FileMetadata for metadata in metadatas)
    @test all(metadata.size == (i-1) * length("Hello, World!\n")
              for (i,metadata) in enumerate(metadatas))
    entries = files_list_folder(auth, "/$folder", recursive=true)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == numfiles + 5
    @test count(entry -> entry.path_display == "/$folder", entries) == 1
    for i in 0:numfiles-1
        @test count(entry -> entry.path_display == "/$folder/files$i",
                    entries) == 1
    end
end

@testset "Download files" begin
    for i in 0:numfiles-1
        metadata, content = files_download(auth, "/$folder/files$i")
        @test metadata.path_display == "/$folder/files$i"
        @test String(content) == repeat("Hello, World!\n", i)
    end
end

@testset "Upload zero files" begin
    upload_spec_channel = Channel{UploadSpec}(0)
    upload_task = start_task(() -> files_upload(auth, upload_spec_channel),
                             ("Upload zero files", :files_upload))
    close(upload_spec_channel)
    metadatas = fetch(upload_task)
    @test isempty(metadatas)
    entries = files_list_folder(auth, "/$folder", recursive=true)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == numfiles + 5
end

@testset "Delete folder" begin
    files_delete(auth, "/$folder")
    entries = files_list_folder(auth, "", recursive=false)
    @test count(entry -> startswith(entry.path_display, "/$folder"),
                entries) == 0
end


text:
import LessOLS: normal_sampler, ols, show, desc, yhat, r2, equation

function almost_equals(a, b, threshold)
    diff = map(abs, (a-b))
    return all(diff .< threshold)
end

function one_lm_normal(n=100)
    gen = normal_sampler(a=0, b=50, β_0=5, β=[1, 2.5], sd_e=0.1)
    sam = gen(n)
    return ols(sam, intercept=true)
end

@testset "sampler_normal: estimated coef are close to true values" begin
    lm = one_lm_normal()
    @test almost_equals(lm.beta, [5, 1, 2.5], 0.1)
end

@testset "methods yhat(), r2() and desc() are callable" begin
    lm = one_lm_normal()
    @test isa(yhat(lm), Vector)
    @test 0<r2(lm)<1
    @test occursin("Linear model", desc(lm))
end

@testset "equation formatter returns spaced string" begin
    text = equation(true, [-1/3, -2/7, 1/29])
    @test text == "Y = - 0.333 - 0.286*X1 + 0.034*X2"
end

@testset "values of beta, ols without intercept " begin
    ns = normal_sampler(a=0, b=20, β_0=5, β=[1, 2.5], sd_e=0.5)
    sample1 = ns(100000)
    lm1 = ols(sample1, intercept=false)
    @test almost_equals(lm1.beta, [1.2143, 2.713], 0.01)
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(layer(Geom.line, x=1:10, y=1:10, Theme(default_color=colorant"red")),
     layer(Geom.point, Geom.label, x=[5], y=[5], label=["5"]))


text:
function matrix_projection(χ::AbstractClassFunction)
    U = matrix_projection(values(χ), conjugacy_classes(χ))
    deg = degree(χ)
    ordG = sum(length, conjugacy_classes(χ))
    return deg // ordG, U
end

function matrix_projection(
    vals::AbstractVector{T},
    ccls::AbstractVector{<:AbstractOrbit{<:Perm}},
    d::Integer = degree(first(first(ccls))),
) where {T}

    result = zeros(T, d, d)

    for (val, cc) in zip(vals, ccls)
        for g in cc
            for i = 1:d
                result[i, i^g] += val
            end
        end
    end

    return result
end

function _conjugate_pairs(chars::Vector{<:AbstractClassFunction})
    vals = values.(chars)
    tovisit = trues(length(vals))
    conjugate_pairs = Dict{Int, Int}()
    for (i, v) in enumerate(vals)
        tovisit[i] || continue
        k = findfirst(==(conj(v)), vals)
        conjugate_pairs[i] = k
        tovisit[i] = false
        tovisit[k] = false
    end
    conjugate_pairs
end

function real_vchars(chars::AbstractVector{<:AbstractClassFunction})
    pairs = _conjugate_pairs(chars)
    res = [VirtualCharacter(chars[i]) for (i,j) in pairs if i == j] # real ones
    for (i,j) in pairs
        i == j && continue
        χ, χ_bar = chars[i], chars[j]
        push!(res, χ + χ_bar)
        push!(res, -E(4)*(χ - χ_bar))
    end

    return res
end


text:
using QuadratureRules
using Test

include("test_utils.jl")
include("test_quadrature_rule.jl")
include("test_gauss_chebyshev.jl")
include("test_gauss_legendre.jl")
include("test_lobatto_chebyshev.jl")
include("test_lobatto_legendre.jl")
include("test_clenshaw_curtis.jl")
include("test_tabulated_quadratures.jl")


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from Euler angle and axis to DCM.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/angleaxis_to_dcm.jl
# ===========================================

# Functions: angleaxis_to_dcm
# ---------------------------

@testset "Euler angle and axis => DCM (Float64)" begin
    T = Float64

    # Sample a random Euler angle and axis.
    v = @SVector randn(T, 3)
    v = v / norm(v)
    a = _rand_ang(T)
    av = EulerAngleAxis(a, v)
    D = angleaxis_to_dcm(av)
    @test eltype(D) === T

    # Test the conversion using the theory.
    sθ, cθ = sincos(a)
    aux = 1 - cθ
    @test D[1, 1] ≈ cθ + v[1] * v[1] * aux
    @test D[1, 2] ≈ v[1] * v[2] * aux + v[3] * sθ
    @test D[1, 3] ≈ v[1] * v[3] * aux - v[2] * sθ
    @test D[2, 1] ≈ v[1] * v[2] * aux - v[3] * sθ
    @test D[2, 2] ≈ cθ + v[2] * v[2] * aux
    @test D[2, 3] ≈ v[2] * v[3] * aux + v[1] * sθ
    @test D[3, 1] ≈ v[1] * v[3] * aux + v[2] * sθ
    @test D[3, 2] ≈ v[2] * v[3] * aux - v[1] * sθ
    @test D[3, 3] ≈ cθ + v[3] * v[3] * aux
end

@testset "Euler angle and axis => DCM (Float32)" begin
    T = Float32

    # Sample a random Euler angle and axis.
    v = @SVector randn(T, 3)
    v = v / norm(v)
    a = _rand_ang(T)
    av = EulerAngleAxis(a, v)
    D = angleaxis_to_dcm(av)
    @test eltype(D) === T

    # Test the conversion using the theory.
    sθ, cθ = sincos(a)
    aux = 1 - cθ
    @test D[1, 1] ≈ cθ + v[1] * v[1] * aux
    @test D[1, 2] ≈ v[1] * v[2] * aux + v[3] * sθ
    @test D[1, 3] ≈ v[1] * v[3] * aux - v[2] * sθ
    @test D[2, 1] ≈ v[1] * v[2] * aux - v[3] * sθ
    @test D[2, 2] ≈ cθ + v[2] * v[2] * aux
    @test D[2, 3] ≈ v[2] * v[3] * aux + v[1] * sθ
    @test D[3, 1] ≈ v[1] * v[3] * aux + v[2] * sθ
    @test D[3, 2] ≈ v[2] * v[3] * aux - v[1] * sθ
    @test D[3, 3] ≈ cθ + v[3] * v[3] * aux
end

@testset "Euler angle and axis => DCM (Errors)" begin
    @test_throws ArgumentError angleaxis_to_dcm(0, [1, 2])
    @test_throws ArgumentError angleaxis_to_dcm(0, [1, 2, 3, 4])
end


text:
# ----------------------------------------------------
# Handlers
# Messages can come from either the "Simulation" or the Client.
# If it comes from the sim then we route to the socket for transmission
# to the client.
# If it comes from the client then we route to either the server
# or the sim.
# ----------------------------------------------------
function handle_msg(soc::SockServer, data::Dict{String,Any})
    if data["From"] == "Client" && data["To"] == "Server"
        handle_client_to_server(soc, data)
        return
    end

    if data["From"] == "Server" && data["To"] == "Client"
        if isopen(soc.socket)
            println(soc.socket, JSON.json(data))
        else
            println("Server-to-Client must have shutdown. Can't send.")
        end
        return
    end

    if data["From"] == "Simulation" && data["To"] == "Client"
        if isopen(soc.socket)
            println(soc.socket, JSON.json(data))
        else
            println("Simulation-to-Client must have shutdown. Can't send.")
        end
        return
    end
end

function handle_client_to_server(soc::SockServer, data::Dict{String,Any})
    if data["Type"] == "Cmd"
        if data["Data"] == "Shutdown Server"
            # Get protocol
            data = JSON.parsefile("../data/com_protocol_basic.json")

            # Populate
            data["From"] = "Server"
            data["To"] = "Client"
            data["Type"] = "Response"
            data["Data"] = "Shutdown in progress"

            # respond back to client
            if isopen(soc.socket)
                println(soc.socket, JSON.json(data))
            else
                println("Client-to-Server must have shutdown. Can't send.")
            end
        elseif data["Data"] == "Simulate"
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            # --------------- Main entry into simulation -------------
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

            # Simulation.run(soc.chan, data["Data1"])

            # This version is used during development.
            Simulation.run_debug(soc.chan, data["Data1"])
        elseif data["Data"] == "Re-Simulate"
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            # --------------- Secondary entry into simulation -------------
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

            # Simulation.re_run(soc.chan, data["Data1"])

            # This version is used during development.
            Simulation.re_run_debug(soc.chan, data["Data1"])
        else
            println("Unknown Cmd: ", data)
        end
    elseif data["Type"] == "Response"
        # Client has acknowledged the shutdown
        soc.running = false
    end
end


text:
"""

    EarlyStopper(c...; verbosity=0)

Instantiate an object for tracking whether one or more stopping
criterion `c` apply, given a sequence of losses.

For a list of possible criterion, do
`subtypes(EarlyStopping.StoppingCriterion)`.

### Sample usage

    stopper = EarlyStopper(Patience(1), NotANumber())
    done!(stopper, 0.123) # false
    done!(stopper, 0.234) # true

    julia> message(stopper)
    "Early stop triggered by Patience(1) stopping criterion. "

### Training losses

For criteria tracking both an "out-of-sample" loss and a "training"
loss (eg, stopping criterion of type `PQ`), specify `training=true` if
the update is for training, as in

    done!(stopper, 0.123; training=true)

Zero or more training updates may precede each out-of-sample update.

The state of the stopper can be reset or restored to a prior state using `reset!`

"""
mutable struct EarlyStopper{S}
    criterion::S
    verbosity::Int
    state
    EarlyStopper(criterion::S; verbosity=0) where S =
        new{S}(criterion, verbosity, nothing)
end
EarlyStopper(criteria...; kwargs...) = EarlyStopper(sum(criteria); kwargs...)

# Dispatch message, done, update and update_training to wrapped criterion
for f in [:message, :done]
    @eval $f(stopper::EarlyStopper) = $f(stopper.criterion, stopper.state)
end
for f in [:update, :update_training]
    @eval $f(stopper::EarlyStopper, loss) =
        $f(stopper.criterion, loss, stopper.state)
end

"""
    done!(stopper::EarlyStopper, loss; training = false)
"""
function done!(stopper::EarlyStopper, loss; training=false)
    if training
        stopper.state = update_training(stopper, loss)
    else
        stopper.state = update(stopper, loss)
    end
    if stopper.verbosity > 0
        suffix = training ? "training " : ""
        loss_str = suffix*"loss"
        @info "$loss_str: $loss\t state: $(stopper.state)"
    end
    return done(stopper)
end

"""
    reset!(stopper::EarlyStopper)
    reset!(stopper::EarlyStopper, state)

Reset a stopper to it's uninitialized state or to a particular state
"""
reset!(stopper::EarlyStopper) = stopper.state = nothing
reset!(stopper::EarlyStopper, state) = stopper.state = state


text:
export Geometry

struct Geometry

    n1 :: Int64
    n2 :: Int64

    x1min  :: Float64
    x2min  :: Float64

    delta1 :: Float64
    delta2 :: Float64

    x1grid :: Array{Float64,1}
    x2grid :: Array{Float64,1}

    bc :: Symbol

"""
initialize geometry from origin cordinates,
space steps and number of points
"""
    function Geometry(n1, n2, x1min, x2min, delta1, delta2)

        x1grid = collect(range(x1min, length=n1, step=delta1))
        x2grid = collect(range(x2min, length=n2, step=delta2))

        new(n1, n2, x1min, x2min, delta1, delta2, x1grid, x2grid, :none)

    end

"""
initialize geometry from origin and end coordinates and
number of points
"""
    function Geometry(n1, n2, x1, x2, bc::Symbol)

	x1min, x1max = x1
	x2min, x2max = x2

        if ( (bc == :perxy) || (bc == :perx) )
           x1grid = collect(range(x1min, stop=x1max, length=n1+1)[1:end-1])
        else
           x1grid = collect(range(x1min, stop=x1max, length=n1))
        end

        if ( (bc == :perxy) || (bc == :pery) )
           x2grid = collect(range(x2min, stop=x2max, length=n2+1)[1:end-1])
        else
           x2grid = collect(range(x2min, stop=x2max, length=n2))
        end

        delta1 = x1grid[2] - x1grid[1]
        delta2 = x2grid[2] - x2grid[1]

        new(n1, n2, x1min, x2min, delta1, delta2, x1grid, x2grid, bc)

    end

end


text:
export nodeCompositionPlot

function nodeCompositionPlot( mpSim::MPsim, timeGrid::Vector{T},
    nodes::AbstractString...; plotType::Symbol=:normal, showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="",
    timeFactor::Real=12.0 ) where T <: Real

    savePlot |= filename != ""
    plotType = haskey( plotTypes, plotType ) ? plotType : :normal
    plotFunction = plotTypes[plotType][2]

    # Don't do anything if the plot mustn't be shown or saved.
    if !( showPlot || savePlot )
        return
    end  # if !( showPlot || savePlot )

    # Generate the report.
    compReport = nodeCompositionReport( mpSim, timeGrid, nodes... )

    if isempty( compReport )
        return
    end  # if isempty( compReport )

    generateCompPlot( compReport, nodes, timeFactor, showPlot, savePlot,
        filename, plotFunction )

end  # nodeCompositionPlot( mpSim, timeGrid, nodes, plotType, showPlot,
     #   savePlot, filename, timeFactor )

nodeCompositionPlot( mpSim::MPsim, timeRes::Real,
    nodes::AbstractString...; plotType::Symbol=:normal, showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="", timeFactor::Real=12.0 ) =
    nodeCompositionPlot( mpSim, MP.generateTimeGrid( mpSim, timeRes ),
    nodes..., plotType=plotType, showPlot=showPlot, savePlot=savePlot,
    filename=filename, timeFactor=timeFactor )

include( "private/nodecomposition.jl" )


text:
module Tmp

    # add exported definitions below
    export get_baseline_points

    # get user input with a custom msg
    function get_baseline_points()
        println("Enter desired points for baseline correction, separated by spaces:")
        string = readline()
        splitstring = split(string, " ")
        return [parse(Int, i) for i in splitstring]
    end

end

max_stars_repo_name: mtmcrobie/ultra-proc

max_issues_repo_name: mtmcrobie/ultra-proc

max_forks_repo_name: mtmcrobie/ultra-proc


text:
@generated function isulp{N,ES}(x::Vnum{N, ES})
  inc = increment(Vnum{N, ES})
  :( (@u(x) & $inc) != 0 )
end

@generated function Base.isempty{N,ES}(x::Valid{N,ES})
  inc = increment(Vnum{N,ES})
  quote
    ((x.lower == zero(Vnum{N,ES})) || (x.upper == zero(Vnum{N,ES}))) && (@u(x.lower) == @u(x.upper) + $inc)
  end
end

Base.isnan{N,ES}(x::Valid{N,ES}) = isempty(x)

@generated function isallreals{N,ES}(x::Valid{N,ES})
  inc = increment(Vnum{N,ES})
  quote
    (x.lower != zero(Vnum{N,ES})) && (x.upper != zero(Vnum{N,ES})) && (@u(x.lower) == @u(x.upper) + $inc)
  end
end

istile{N,ES}(x::Valid{N,ES}) = (x.lower == x.upper)

doc"""
    roundsinf(x::Valid)
    checks if the value "rounds infinity".  Table:

    ===================================
    [(+, 0])      |           true
    [(+, +])      |  lower > upper
    [(+, Inf)     |          false
    [(+, Inf]     |           true
    [(+, -])      |           true
    [Inf, 0])     |           true
    (Inf, 0])     |          false
    [Inf, +])     |           true
    (Inf, +])     |          false
    [Inf, -])     |           true
    (Inf, -])     |          false
    [Inf, Inf]    |           true
    (Inf, Inf]    |           true (ℝp)
    [Inf, Inf)    |           true (ℝp)
    (Inf, Inf)    |          false
    [(-, 0])      |          false
    [(-, +])      |          false
    [(-, Inf)     |          false
    [(-, Inf]     |           true
    [(-, -])      |  lower > upper
    [0, 0]        |          false
    [0, 0)        |          false (∅, may be mapped to noisy NaN in the future)
    (0, 0]        |          false (∅)
    (0, 0)        |           true
    [(0, +])      |          false
    [(0, Inf)     |          false
    [(0, Inf]     |           true
    [(0, -])      |           true
"""
function roundsinf{N,ES}(x::Valid{N,ES})
    (!isempty(x)) && ((@s(x.upper) < @s(x.lower)) || isinf(x.lower))
end

#amazing symmetrical structure to "roundsinf"
function containszero{N,ES}(x::Valid{N,ES})
    (!isempty(x)) && ((@u(x.upper) < @u(x.lower)) || x.lower == zero(Vnum{N,ES}))
end

function ispositive{N,ES}(x::Valid{N,ES})
  zero(@Int) < (@s x.lower) <= (@s x.upper)
end

function isnegative{N,ES}(x::Valid{N,ES})
  (@u Vnum{N,ES}(Inf))< (@u x.lower) <= (@u x.upper)
end

"""
  nonnegative(::Valid) is true if no values in x are negative.
"""
nonnegative{N,ES}(x::Valid{N,ES}) = (x.lower <= zero(Vnum{N,ES})) && ((!isfinite(x.upper) || x.lower <= x.upper <= maxpos(Vnum{N,ES})))

"""
  nonpositive(::Valid) is true if no values in x are positive.
"""
nonpositive{N,ES}(x::Valid{N,ES}) = (x.upper <= zero(Vnum{N,ES})) && ((!isfinite(x.lower) || minneg(Vnum{N,ES}) <= x.lower <= x.upper))

"""
  rounds_positive(::Valid) is true if x contains both zero, infinity, and passes through the positive numbers.
"""
rounds_positive{N,ES}(x::Valid{N,ES}) = isfinite(x.lower) && (x.lower <= zero(Vnum{N,ES})) && (x.upper < x.lower)

"""
  rounds_negative(::Valid) is true if x contains both zero, infinity, and passes through the negative numbers.
"""
rounds_negative{N,ES}(x::Valid{N,ES}) = isfinite(x.upper) && (x.upper >= zero(Vnum{N,ES})) && (x.lower > x.upper)

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
export iauFk45z
"""
Convert a B1950.0 FK4 star position to J2000.0 FK5, assuming zero
proper motion in the FK5 system.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

This function converts a star's catalog data from the old FK4
(Bessel-Newcomb) system to the later IAU 1976 FK5 (Fricke) system,
in such a way that the FK5 proper motion is zero.  Because such a
star has, in general, a non-zero proper motion in the FK4 system,
the routine requires the epoch at which the position in the FK4
system was determined.

Given:
   - r1950,d1950    double   B1950.0 FK4 RA,Dec at epoch (rad)
   - bepoch         double   Besselian epoch (e.g. 1979.3D0)

Returned:
   - r2000,d2000    double   J2000.0 FK5 RA,Dec (rad)

Notes:

   1. The epoch bepoch is strictly speaking Besselian, but if a
   Julian epoch is supplied the result will be affected only to a
   negligible extent.

   2. The method is from Appendix 2 of Aoki et al. (1983), but using
   the constants of Seidelmann (1992).  See the routine iauFk425
   for a general introduction to the FK4 to FK5 conversion.

   3. Conversion from equinox B1950.0 FK4 to equinox J2000.0 FK5 only
   is provided for.  Conversions for different starting and/or
   ending epochs would require additional treatment for precession,
   proper motion and E-terms.

   4. In the FK4 catalog the proper motions of stars within 10 degrees
   of the poles do not embody differential E-terms effects and
   should, strictly speaking, be handled in a different manner from
   stars outside these regions.  However, given the general lack of
   homogeneity of the star data available for routine astrometry,
   the difficulties of handling positions that may have been
   determined from astrometric fields spanning the polar and non-
   polar regions, the likelihood that the differential E-terms
   effect was not taken into account when allowing for proper motion
   in past astrometry, and the undesirability of a discontinuity in
   the algorithm, the decision has been made in this SOFA algorithm
   to include the effects of differential E-terms on the proper
   motions for all stars, whether polar or not.  At epoch 2000.0,
   and measuring "on the sky" rather than in terms of RA change, the
   errors resulting from this simplification are less than
   1 milliarcsecond in position and 1 milliarcsecond per century in
   proper motion.

References:

   Aoki, S. et al., 1983, "Conversion matrix of epoch B1950.0
   FK4-based positions of stars to epoch J2000.0 positions in
   accordance with the new IAU resolutions".  Astron.Astrophys.
   128, 263-267.

   Seidelmann, P.K. (ed), 1992, "Explanatory Supplement to the
   Astronomical Almanac", ISBN 0-935702-68-7.

Called:
   - iauAnp       normalize angle into range 0 to 2pi
   - iauC2s       p-vector to spherical
   - iauEpb2jd    Besselian epoch to Julian date
   - iauEpj       Julian date to Julian epoch
   - iauPdp       scalar product of two p-vectors
   - iauPmp       p-vector minus p-vector
   - iauPpsp      p-vector plus scaled p-vector
   - iauPvu       update a pv-vector
   - iauS2c       spherical to p-vector

This revision:   2018 December 5

SOFA release 2019-07-22

Copyright (C) 2019 IAU SOFA Board.  See notes at end.
"""

# void iauFk45z(double r1950, double d1950, double bepoch,
#               double *r2000, double *d2000)

function iauFk45z(r1950::Real, d1950::Real, bepoch::Real)
   ref_r2000 = Ref{Float64}(0.0)
   ref_d2000 = Ref{Float64}(0.0)

   ccall((:iauFk45z, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble,
         Ref{Cdouble}, Ref{Cdouble}, ),
         convert(Float64, r1950), convert(Float64, d1950),
         convert(Float64, bepoch),
         ref_r2000, ref_d2000)

   return ref_r2000[], ref_d2000[]
end


text:
#!/usr/bin/env julia

# Packages must be loaded in this order!
using MPI
using HDF5
using PencilArrays
using PencilArrays.PencilIO

import JSON3

if !PencilIO.hdf5_has_parallel()
    @warn "HDF5 has no parallel support. Skipping HDF5 tests."
    exit(0)
end

using Random
using Test

include("include/MPITools.jl")
using .MPITools

function test_write_mpiio(filename, u::PencilArray)
    comm = get_comm(u)
    root = 0
    rank = MPI.Comm_rank(comm)

    X = (u, u .+ 1, u .+ 2, u .+ 3, u .+ 4)

    kws = Iterators.product((false, true), (false, true))

    @test_nowarn open(MPIIODriver(), filename, comm, write=true, create=true) do ff
        pos = 0
        for (i, (collective, chunks)) in enumerate(kws)
            name = "field_$i"
            ff[name, collective=collective, chunks=chunks] = X[i]
            pos += sizeof_global(X[i])
            @test position(ff) == pos
        end
    end

    # Append some data.
    open(MPIIODriver(), filename, comm, write=true, append=true) do ff
        ff["field_5", chunks=false] = X[5]
        ff["collection"] = X
    end

    @test isfile("$filename.json")
    meta = open(JSON3.read, "$filename.json", "r").datasets

    # Test file contents in serial mode.
    # First, gather data from all processes.
    # Note that we may need to permute indices of data, since data on disk is
    # written in memory (not logical) order.
    perm = Tuple(permutation(u))
    Xg = map(X) do x
        xg = gather(x, root)  # note: data is in logical order
        xg === nothing && return xg
        perm === nothing && return xg
        PermutedDimsArray(xg, perm)
    end

    @test (Xg[1] === nothing) == (rank != root)
    if rank == root
        open(filename, "r") do ff
            y = similar(Xg[1])
            for (i, (collective, chunks)) in enumerate(kws)
                mpiio_read_serial!(ff, y, meta, "field_$i")
                if !chunks  # if chunks = true, data is reordered into blocks
                    @test y == Xg[i]
                end
            end
            # Verify appended data
            mpiio_read_serial!(ff, y, meta, "field_5")
            @test y == Xg[5]
            let y = similar.(Xg)
                mpiio_read_serial!(ff, y, meta, "collection")
                @test all(y .== Xg)
            end
        end
    end

    # Read stuff
    y = similar(X[1])
    @test_nowarn open(MPIIODriver(), filename, comm, read=true) do ff
        @test_throws ErrorException read!(ff, y, "field not in file")
        let y = similar.(X)
            read!(ff, y, "collection")
            @test all(y .== X)
        end
        for (i, (collective, _)) in enumerate(kws)
            name = "field_$i"
            read!(ff, y, name, collective=collective)
            @test y == X[i]
        end
        read!(ff, y, "field_5")
        @test y == X[5]
    end

    nothing
end

read_array!(ff, x) = read!(ff, x)
read_array!(ff, t::Tuple) = map(x -> read_array!(ff, x), t)

function mpiio_read_serial!(ff, x, meta, name)
    offset = meta[Symbol(name)].offset_bytes :: Int
    seek(ff, offset)
    read_array!(ff, x)
end

function test_write_hdf5(filename, u::PencilArray)
    comm = get_comm(u)
    rank = MPI.Comm_rank(comm)
    v = u .+ 1
    w = u .+ 2

    # Open file in serial mode first.
    if rank == 0
        h5open(filename, "w") do ff
            # "HDF5 file was not opened with the MPIO driver"
            @test_throws ErrorException ff["scalar"] = u
        end
    end

    MPI.Barrier(comm)

    @test_nowarn open(PHDF5Driver(), filename, comm, write=true) do ff
        @test isopen(ff)
        @test_nowarn ff["scalar", collective=true, chunks=false] = u
        @test_nowarn ff["vector_tuple", collective=false, chunks=true] = (u, v, w)
        @test_nowarn ff["vector_array", collective=true, chunks=true] = [u, v, w]
    end

    @test_nowarn open(PHDF5Driver(), filename, comm, append=true) do ff
        @test isopen(ff)
        @test_nowarn ff["scalar_again"] = u
    end

    @test_nowarn open(PHDF5Driver(), filename, comm, read=true) do ff
        @test isopen(ff)
        uvw = (u, v, w)
        uvw_r = similar.(uvw)
        ur, vr, wr = uvw_r

        read!(ff, ur, "scalar")
        @test u == ur
        read!(ff, vr, "scalar_again")
        @test vr == ur

        let perm = Tuple(permutation(ur))
            @test haskey(attributes(ff["scalar"]), "permutation")
            expected = perm === nothing ? false : collect(perm)
            @test read(ff["scalar"]["permutation"]) == expected
        end

        read!(ff, uvw_r, "vector_tuple")
        @test all(uvw .== uvw_r)

        fill!.(uvw_r, 0)
        read!(ff, collect(uvw_r), "vector_array", collective=false)
        @test all(uvw .== uvw_r)
    end

    nothing
end

function main()
    MPI.Init()

    Nxyz = (16, 21, 41)
    comm = MPI.COMM_WORLD
    Nproc = MPI.Comm_size(comm)
    myrank = MPI.Comm_rank(comm)

    silence_stdout(comm)

    @show HDF5.libhdf5

    rng = MersenneTwister(42)
    perms = (NoPermutation(), Permutation(2, 3, 1))

    @testset "$perm" for perm in perms
        pen = Pencil(Nxyz, (1, 3), comm; permute = perm)
        u = PencilArray{Float64}(undef, pen)
        randn!(rng, u)
        u .+= 10 * myrank

        @testset "MPI-IO" begin
            filename = MPI.bcast(tempname(), 0, comm)
            test_write_mpiio(filename, u)
        end

        @testset "HDF5" begin
            filename = MPI.bcast(tempname(), 0, comm)
            test_write_hdf5(filename, u)
        end
    end

    HDF5.h5_close()
    MPI.Finalize()
end

main()


text:
@testset "f increase behaviour" begin
    f(x) = 2*x[1]^2
    g!(G, x) = copyto!(G, 4*x[1])

    # Returned "minimizers" from one iteration of Gradient Descent
    minimizers = [0.3,  # allow_f_increases = true, alpha = 0.1
                  -1.5, # allow_f_increases = true, alpha = 1.0
                  0.3,  # allow_f_increases = false, alpha = 0.1
                  0.5]  # allow_f_increases = false, alpha = 1.0
    k = 0
    for allow in [true, false]
        for alpha in [0.1, 1.0]
            k += 1
            method = GradientDescent(
                alphaguess = LineSearches.InitialStatic(alpha=alpha),
                linesearch = LineSearches.Static(),
            )
            opts = Optim.Options(iterations=1,allow_f_increases=allow)
            res = optimize(f, g!, [0.5], method, opts)

            @test minimizers[k] == Optim.minimizer(res)[1]
        end
    end
end


text:
# NOTE: there's a `kmeans!` function that updates centers, maybe a candidate
# for the `update` machinery. Same for `kmedoids!`
# NOTE: if the prediction is done on the original array, just the assignment
# should be returned, unclear what's the best way of doing this.

module Clustering_

export KMeans
export KMedoids

import MLJBase

import ..Clustering # strange sytax for lazy-loading

using Distances
using LinearAlgebra: norm

const C = Clustering

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

const KMFitResultType = AbstractMatrix{<:Real} # the medoids

mutable struct KMeans{M<:SemiMetric} <: MLJBase.Unsupervised
    k::Int
    metric::M
end

mutable struct KMedoids{M<:SemiMetric} <: MLJBase.Unsupervised
    k::Int
    metric::M
end

function MLJBase.clean!(model::Union{KMeans, KMedoids})
    warning = ""
    if model.k < 2
        warning *= "Need k >= 2. Resetting k=2.\n"
        model.k = 2
    end
    return warning
end

####
#### KMEANS: constructor, fit, transform and predict
####

function KMeans(; k=3, metric=SqEuclidean())
    model = KMeans(k, metric)
    message = MLJBase.clean!(model)
    isempty(message) || @warn message
    return model
end

function MLJBase.fit(model::KMeans
                   , verbosity::Int
                   , X)

    Xarray = MLJBase.matrix(X)

    # NOTE see https://github.com/JuliaStats/Clustering.jl/issues/136
    # this has been updated but only on #master, in the future replace permutedims
    # with transpose (lazy)
    result    = C.kmeans(permutedims(Xarray), model.k; distance=model.metric)
    fitresult = result.centers # centers (p x k)
    cache     = nothing
    report    = Dict(:assignments => result.assignments) # size n

    return fitresult, cache, report
end

function MLJBase.transform(model::KMeans
                         , fitresult::KMFitResultType
                         , X)

    Xarray = MLJBase.matrix(X)
    (n, p), k = size(Xarray), model.k
    # pairwise distance from samples to centers
    X̃ = pairwise(model.metric, transpose(Xarray), fitresult)
    return MLJBase.table(X̃, prototype=X)
end

####
#### KMEDOIDS: constructor, fit and predict
#### NOTE there is no transform in the sense of kmeans
####

function KMedoids(; k=3, metric=SqEuclidean())
    model = KMedoids(k, metric)
    message = MLJBase.clean!(model)
    isempty(message) || @warn message
    return model
end

function MLJBase.fit(model::KMedoids
                   , verbosity::Int
                   , X)

    Xarray = MLJBase.matrix(X)
    # cost matrix: all the pairwise distances
    Carray = pairwise(model.metric, transpose(Xarray)) # n x n

    result    = C.kmedoids(Carray, model.k)
    fitresult = permutedims(view(Xarray, result.medoids, :)) # medoids (p x k)
    cache     = nothing
    report    = Dict(:assignments => result.assignments) # size n

    return fitresult, cache, report
end

function MLJBase.transform(model::KMedoids
                         , fitresult::KMFitResultType
                         , X)

    Xarray = MLJBase.matrix(X)
    (n, p), k = size(Xarray), model.k
    # pairwise distance from samples to medoids
    X̃ = pairwise(model.metric, transpose(Xarray), fitresult)
    return MLJBase.table(X̃, prototype=X)
end

####
#### Predict methods
####

function MLJBase.predict(model::Union{KMeans,KMedoids}, fitresult::KMFitResultType, Xnew)

    Xarray = MLJBase.matrix(Xnew)
    (n, p), k = size(Xarray), model.k

    pred = zeros(Int, n)
    @inbounds for i ∈ 1:n
        minv = Inf
        for j ∈ 1:k
            curv    = evaluate(model.metric,
                                  view(Xarray, i, :), view(fitresult, :, j))
            P       = curv < minv
            pred[i] =    j * P + pred[i] * !P # if P is true --> j
            minv    = curv * P +    minv * !P # if P is true --> curvalue
        end
    end
    return pred
end

####
#### METADATA
####

MLJBase.load_path(::Type{<:KMeans}) = "MLJModels.Clustering_.KMeans" # lazy-loaded from MLJ
MLJBase.package_url(::Type{<:KMeans}) = "https://github.com/JuliaStats/Clustering.jl"
MLJBase.package_name(::Type{<:KMeans}) = "Clustering"
MLJBase.package_uuid(::Type{<:KMeans}) = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
MLJBase.is_pure_julia(::Type{<:KMeans}) = :yes
MLJBase.input_kinds(::Type{<:KMeans}) = [:continuous,]
MLJBase.output_kind(::Type{<:KMeans}) = :continuous
MLJBase.output_quantity(::Type{<:KMeans}) = :multivariate

MLJBase.load_path(::Type{<:KMedoids}) = "MLJModels.Clustering_.KMedoids" # lazy-loaded from MLJ
MLJBase.package_url(::Type{<:KMedoids}) = MLJBase.package_url(KMeans)
MLJBase.package_name(::Type{<:KMedoids}) = MLJBase.package_name(KMeans)
MLJBase.package_uuid(::Type{<:KMedoids}) = MLJBase.package_uuid(KMeans)
MLJBase.is_pure_julia(::Type{<:KMedoids}) = :yes
MLJBase.input_kinds(::Type{<:KMedoids}) = [:continuous,]
MLJBase.output_kind(::Type{<:KMedoids}) = :continuous
MLJBase.output_quantity(::Type{<:KMedoids}) = :multivariate

end # module

max_stars_repo_name: ablaom/MLJModels

max_issues_repo_name: ablaom/MLJModels

max_forks_repo_name: ablaom/MLJModels


text:
using Pkg
Pkg.activate(".")

using BenchmarkTools

x = rand(20000, 20000)

function ssum(x)
    r, c = size(x)
    y = zeros(c)
    for i in 1:c
        for j in 1:r
            y[i] += x[j, i]
        end
    end
    y
end

function tsum(x)
    r, c = size(x)
    y = zeros(c)
    Threads.@threads for i in 1:c
        for j in 1:r
            y[i] += x[j, i]
        end
    end
    y
end

println("threads: ", Threads.nthreads())
@time ssum(x)
@time ssum(x)
@time tsum(x)
@time tsum(x)