Sampling from (non-)Normal distribution

[luca-fiorito-11]

@AitorBengoechea , @EnricaBelfiore , @GrimFe
When we choose a distribution for sampling that can go in the negative range we should apply some sort of truncation (even if the pdf is normal).
Examples if perturbation coefficients are centered in 1:
• Any coefficient <0 or >2 will be reset respectively to 0 or 2
• Any coefficient <0 or >2 will be reset to 1
• Any coefficient <0 or >2 will be resampled

I guess this should be handled in the CategoryCov.sampling algorithm.

Please provide feedback

[EnricaBelfiore]

As discussed in a previous meeting with @luca-fiorito-11 and @GrimFe, the first option looks promising as it precerves the symmetry and it gives lower changes of the standard deviation with respect to the other options.

I can share a notebook that today @GrimFe and I have developed to analyze how this could be implemented in SANDY. It will be on teams for both of you @luca-fiorito-11 and @AitorBengoechea and I'm now discussing the results:

• Cutting distribution queues at the stage in which the sample population has mean zero and absolute std 1 gives worse results as long as nuclear data relative standard deviations lower than 1 are considered;
• Once you have applied the nuclear data relative std to the sample population, the results will be not affected by the order in which the operations are performed.

I will now focus on thinking on a structured way to implement this in sandy (in matricial form) and to consider Decay Data, for which the absolute standard deviation is given.

[luca-fiorito-11]

Uniform distribution implemented in #145

[AitorBengoechea]

Which method has the smallest relative difference to CatCov?

[EnricaBelfiore]

Edit: Sorry, by accident the previous post in the discussion was deleted. This allowed us to upload a newer version of the notebook.
@AitorBengoechea hopefully this answers your quesition.

@GrimFe and I developed a notebook with a proposal for a maybe even easier sampling procedure, which allows the implementation of both normal and lognormal distribution.
Uniform distribution sampling still requires work.

The proposed implementation allows to:

• avoid matrix decomposition and eigenvalue calculation
• process quite big matrices (@AitorBengoechea would you mind cheking with matrices generated by NJOY?)
• set the negative eigenvalues to zero (just by passing a worked covariance - implemented in the latter part of notebook)
• possibly reduce the number of methods needed in CategoryCov
• introduce the population mean in a different part of SANDY code, as the population mean is imposed to 1 and relative covariances are considered.

You can find the notebook here: https://github.com/GrimFe/SandyProposals/blob/main/random_multivariate.ipynb

As you can see in the notebook, the proposed method was tested for a number of nuclides, among which Eu155. Results for the latter are visible in the notebook since they appear to be rather visually explanatory.
The last section of the notebook shows a comparison with what is now implemented in CategoryCov.sampling(). The results seem to be rather comparable. Some fluctuations corresponding to zeros in the evaluated covariance matrix are found (i.e: the one in the ENDF6 file).

Covarinces of the sampled populations generated with both methods (the proposed one and CategoryCov.sampling() were compared with the evaluated covariance matrix.
We tested 12 nuclides and almost 50% of the times CategoryCov.sampling() gives better representation. This depends on how the comparsion of the covariance matrices is performed (either comparing the largest eigenvalue or the maximum difference in the matrix, namely) and on how the population is distributed (i.e.: it depends on the considered number of samples of course).

[luca-fiorito-11]

This looks really promising!
Indeed if it works we could remove a lot of CategoryCov methods and improve maintainability.

Independently whether it can sample form uniform distributions, this method has the potential to become the default one.
I remind you that sampling from a multivariate normal is the only requirement. Lognormal and uniform are nice-to-have options, if they accommodate well enough within the default method.

A few more things:

• can you check the source/docs to see what decomposition they use, if any?
• can you check how they handle non-SPD matrices?
• can you compare our sampling unit tests when adopting this new method?
• can you check (on a notebook) how this method integrates with our sampling pipeline (real case)?
• can you highlight what existing methods we can remove?

If all these checks are positive, we will integrate it as default

[EnricaBelfiore]

Numpy suggests to use the Random Generator objects for sampling which @GrimFe and I have implemented in our Notebook.

• As we can see in the API reference (https://numpy.org/doc/stable/reference/random/generated/numpy.random.Generator.multivariate_normal.html),
  the input covariance matrix (C) is used to compute a factor matrix L such that L*L.T = C, exactly what we want to
  do in SANDY.
  This new proposed method can use different decompositions: ‘svd’, ‘eigh’, ‘cholesky’. The default one is the
  singular value decomposition, which is also the one implemented in this moment in CategoryCov, followed by QR decomposition to obtain L.
• They seem to handle non-SPD matrices by setting the negative eigenvalues equal to their absolute values (see more detailed discussion below).
• In the following notebook we can see the comparison with the sampling unit tests, it gives coherent results: https://github.com/GrimFe/SandyProposals/blob/main/sampling_distribution_covbyNJOY.ipynb

Concerning the pipeline and the list of methods to be removed from CategoryCov, we need to decide how to approach the negative eigenvalues problem.

• keep CategoryCov as it is, which does not allow for easy implementation of lognormal distribution sampling and needs all the methods now implemented;
• use the random generator multivariate normal method with nuclear data covariance, so without setting the eigenvalues to zero. This allows an easy implementation of the lognormal distribution sampling but can pose problems with the negative eigenvalues: namely non SPD handling is instable and we did not manage to predict how the eigenvalues are modified for big matrices. As mantioned above, tests on small matrices seem to suggest taking the absolute values of negative eigenvalues to be an option.
• use the random generator multivariate normal method with nuclear data covariance recomposed through eigen recomposition setting the negative eigenvalues to zero. This seems to be the most promisig option as it allows an easy implementation of the lognormal distribution sampling and avoid negative eigenvalues instabilities, but reduce/restore size are still needed for good eigen decomposition. The warning that we can see in the notebook given by non- SPD matrix is due to the numerical fluctuation around zero (order of magnitude: 10^-16) that should pose no problem, no matter how they are handled by numpy.

The latter option allows to remove the following method from CategoryCov: sampling, sample_distribution, get_L.

[luca-fiorito-11]

@EnricaBelfiore , @GrimFe , can you prepare a project for this work as in https://github.com/luca-fiorito-11/sandy/projects/6 ?