Test suite to analyze the NcmSplineFunc’s knots distribution.

This module is intended to be a test suite for the NcmSplineFunc object. It performs a brute force approach to check if the required tolerance (NcmSplineFuncTest:rel-error and NcmSplineFuncTest:scale) was achieved throughout the entire desired range [NcmSplineFuncTest:xi, NcmSplineFuncTest:xf] to evaluate a base function $f(x)$. To perform the test, it is created an evenly distributed grid with NcmSplineFuncTest:ngrid knots, including both limiting points.

The criteria applied by NcmSplineFunc to include a knot is given by the condition, \begin{equation}\label{eq:condition} |f(x) - \hat{f}(x)| < \mathrm{rel \_ error} \times \left[ |f(x)| + \mathrm{scale} \right] \,\, . \end{equation} Where, $f(x)$ is the base (true) function to be analyzed and $\hat{f}(x)$ is the estimation given by NcmSplineFunc using the interpolation method NcmSplineCubicNotaknot.

The test can be done on any function provided by the user (#NCM_SPLINE_FUNC_TEST_TYPE_USER), but it is also a stress test tool with some built in base functions (see NcmSplineFuncTestType):

• NCM_SPLINE_FUNC_TEST_TYPE_POLYNOMIAL: it applies a polynomial interpolation with the desired degree upon a given set of points.

  The polynomial degree is given by the number of points provided by the user subtracted by one. It uses GSL‘s polynomial interpolation method encapsulated with NcmSplineGsl.

• NCM_SPLINE_FUNC_TEST_TYPE_COSINE: it applies a summation of cosine functions given by

  \begin{equation} f(x) = \sum_{i=0}^{N-1} A_{i} \cos \left( 2 \pi \, \nu_i \, x \right) \,\, . \end{equation} The user provides the amplitudes $A_i$ and the frequencies $\nu_i$.

• NCM_SPLINE_FUNC_TEST_TYPE_RBF: it is similar to the polynomial, but now it applies a RBF interpolation upon a given set of points.

For both interpolated base functions, polynomial and RBF, all ordinate points $y$ can be drawn from a flat or a normal (gaussian) distribution. The abscissa points $x$ are always drawn from a flat distribution, and its limiting points are fixed at NcmSplineFuncTest:xi and NcmSplineFuncTest:xf.

The randomness goal is to produce a base function with several shapes. The NcmSplineFuncTestTypePDF type provides the switch between both available PDFs.

Setting the parameters # {#set-par}

In order to provide the base function parameters, the user has two options:

1. ncm_spline_func_test_set_params_info(): the user must supply a NcmMatrix with the number of rows as the number of parameters and two columns.

   If it is a flat PDF, each column represents the minimum and maximum parameter value possible to be drawn. If it is a normal (gaussian) PDF, each column represents its mean and standard deviation. The supplied matrix is copied to NcmSplineFuncTest:par-info.

2. ncm_spline_func_test_set_params_info_all(): the user must supply the number of parameters and the same value for each column.

   Therefore, all parameters will have the same statistical properties. It creates the matrix NcmSplineFuncTest:par-info internally.

Fixing a parameter: the user also have the option to fix some or all the parameters.

• NCM_SPLINE_FUNC_TEST_TYPE_PDF_FLAT: set the parameter's minimum and maximum to the same value (both columns with the same value).

• NCM_SPLINE_FUNC_TEST_TYPE_PDF_NORMAL: set the parameter's standard deviation to zero (second column equal to zero).

Built in functions parameters:

• For the two interpolation methods, polynomial and RBF, the provided parameters are the ordinate points $y$.

• For the cosine summation method, the provided parameters are the amplitude and the frequency, alternating in that order. If someone wants to perform a sum of two cosine functions, must provide a matrix with four rows: $A_0$, $\nu_0$, $A_1$ and $\nu_1$, for example. The #NCM_SPLINE_FUNC_TEST_TYPE_COSINE input parameters matrix must have an even number of rows.

One grid examples # {#grid-ex}

To perform one grid statistics, it is needed to place ncm_spline_func_test_set_one_grid_stats() after preparing it (ncm_spline_func_test_prepare()). Two examples are shown below together with their results.

Example: 7th degree polynomial interpolation. # {#7th-ex}

#include <numcosmo/numcosmo.h>

int
main (void)
{
  guint npar = 8, seed = 1, ngrid = 100000;

  gdouble rel_error = 1.e-10, scale = 1.0, mean = 0.0, sigma = 1.0;

  NcmSplineFuncTest *sft = ncm_spline_func_test_new ();

  ncm_spline_func_test_set_seed (sft, seed);

  ncm_spline_func_test_set_ngrid (sft, ngrid);

  ncm_spline_func_test_set_scale (sft, scale);

  ncm_spline_func_test_set_rel_error (sft, rel_error);

  ncm_spline_func_test_set_params_info_all (sft, npar, mean, sigma);

  ncm_spline_func_test_prepare (sft, NCM_SPLINE_FUNCTION_SPLINE, NCM_SPLINE_FUNC_TEST_TYPE_PDF_NORMAL);

  ncm_spline_func_test_set_one_grid_stats (sft);

  ncm_spline_func_test_log_vals_one_grid_stats (sft);

  ncm_spline_func_test_save_grid_functions_to_txt (sft, "functions.txt");

  ncm_spline_func_test_save_knots_to_txt (sft, "knots.txt");

  ncm_spline_func_test_unref (sft);

  return 0;
}

The function ncm_spline_func_test_set_params_info_all() creates the 8 rows matrix, “NcmSplineFuncTest:par-info”, the number of ordinate points $y$ (consequently 8 abscissa $x$), and two columns. The first column is $\mu=0$ and the second $\sigma = 1$ with the same value for all parameters (rows) to be used by the gaussian PDF.

The function ncm_spline_func_test_prepare() sets the NcmSplineFunc method used, NCM_SPLINE_FUNCTION_SPLINE, and also sets the PDF, #NCM_SPLINE_FUNC_TEST_TYPE_PDF_NORMAL.

The functions ncm_spline_func_test_save_grid_functions_to_txt() and ncm_spline_func_test_save_knots_to_txt() create the ascii files “functions.txt”, with all the information about the functions at all grid, and “knots.txt”, which saves NcmSplineFunc‘s knots.

The function ncm_spline_func_test_log_vals_one_grid_stats() prints the statistics evaluated by ncm_spline_func_test_set_one_grid_stats(). In the above example it displays:

Grid statistics

NcmSplineFunc number of knots = 3847

(ncm) diff. = -7.377e-13 +/- 2.582e-11 (abs. max. diff. = 2.917e-10) | outliers = 0.00 % (lin) diff. = -3.826e-13 +/- 3.091e-11 (abs. max. diff. = 1.211e-09) | outliers = 0.09 %

* “NcmSplineFunc number of knots = 3847”: is self explanatory. * “diff.” first value (-7.377e-13 and -3.826e-13): it is the mean signed difference between the base and approximated functions, $f(x)$ and $\hat{f}(x)$. The signed difference is defined as $\Delta f(x)= f(x) - \hat{f}(x)$. Therefore its mean value through all grid knots should be as close to zero as possible, regardless of the required tolerance. * “diff.” second value (2.582e-11 and 3.091e-11): it is the standard deviation of $\Delta f(x)$. This value is related to the required tolerance. If $\Delta f(x)$ distribution is assumed to be normal, which is not, we have $5\sigma \approx 10^{-10}$, the required NcmSplineFuncTest:rel-error in this example. * “abs. max. diff.”: it is the maximum absolute difference between the entire grid knots, $|\Delta f(x)|_{\mathrm{max}}$. In this example, it appears that “ncm” has a higher value than the required tolerance, but in fact it is not true, because scale is not zero. * “outliers”: it is defined as the knots in the linear grid that did not passed the criteria given by Eq. \eqref{eq:condition}. It is given in percentage of the number of knots (NcmSplineFuncTest:ngrid).

Example: cosine summation with 50 terms.

#include <numcosmo/numcosmo.h>

int
main (void)
{
  NcmVector *c           = ncm_vector_new (50);
  NcmMatrix *params      = ncm_matrix_new (50, 2);
  NcmSplineFuncTest *sft = ncm_spline_func_test_new ();

  ncm_vector_set_all (c, -5.0);
  ncm_matrix_set_col (params, 0, c);
  ncm_vector_set_all (c, +5.0);
  ncm_matrix_set_col (params, 1, c);

  ncm_matrix_set (params, 0, 0, 100.); // Sets the first
  ncm_matrix_set (params, 0, 1, 100.); // amplitude fixed: A_0 = 100.
  ncm_matrix_set (params, 1, 0, 0.);   // Sets the first
  ncm_matrix_set (params, 1, 1, 0.);   // frequency fixed: \nu_0 = 0.

  ncm_spline_func_test_set_type (sft, NCM_SPLINE_FUNC_TEST_TYPE_COSINE);

  ncm_spline_func_test_set_params_info (sft, params);

  ncm_spline_func_test_set_ngrid (sft, 1000000);

  ncm_spline_func_test_set_seed (sft, 73649);

  ncm_spline_func_test_prepare (sft, NCM_SPLINE_FUNCTION_4POINTS, NCM_SPLINE_FUNC_TEST_TYPE_PDF_FLAT);

  ncm_spline_func_test_set_one_grid_stats (sft);

  ncm_spline_func_test_log_vals_one_grid_stats (sft);

  ncm_vector_free (c);
  ncm_matrix_free (params);
  ncm_spline_func_test_unref (sft);

  return 0;
}

In this example the user created a matrix with 50 parameters, 25 amplitudes and 25 frequencies. But note that the first amplitude is fixed to $A_0 = 100$ and the first frequency to $\nu_0 = 0$. Therefore, creating a base function with some kind of oscilatory behaviour added to a constant factor, \begin{equation} f(x) = A_0 + \sum_{i=1}^{49} A_i \cos \left( 2 \pi \, \nu_i \, x \right) \,\,\, , \end{equation} with $A_i$ and $\nu_i$ drawn from a flat PDF between the values [-5, 5]. Below is shown the output from ncm_spline_func_test_log_vals_one_grid_stats()

Grid statistics

NcmSplineFunc number of knots = 7318

(ncm) diff. = 2.050e-13 +/- 3.565e-12 (abs. max. diff. = 5.193e-11) | outliers = 0.00 % (lin) diff. = -1.088e-14 +/- 2.517e-12 (abs. max. diff. = 7.323e-11) | outliers = 0.00 %

Those statistics are remarkably better than the required tolerance given by Eq. \eqref{eq:condition}. In this case we have the default values, $\mathrm{rel \ error} = 10^{-8}$ and $\mathrm{scale} = 0$. The result condition is around $|f(x)| \times \mathrm{rel \ error} \approx 10^{-6}$. Note that both maximum absolute error are given by $|\Delta f(x)|_{\mathrm{max}} \approx 10^{-10}$, four orders of magnitude better than expected. This fact is due to the NcmSplineFunc method applied in this example, #NCM_SPLINE_FUNCTION_4POINTS. It is a much more conservative approach compared to #NCM_SPLINE_FUNCTION_SPLINE, applied in the previous example.

Monte Carlo examples # {#mc-ex}

To perform a Monte Carlo statistics it is needed to place ncm_spline_func_test_monte_carlo_and_save_to_txt() or ncm_spline_func_test_monte_carlo() after preparing it (ncm_spline_func_test_prepare()).

Two examples are shown below, together with their results.

Example: 5th degree polynomial interpolation with 1 million realizations.

#include <numcosmo/numcosmo.h>

int
main (void)
{
  guint npar = 6, nsim = 1000000;

  NcmSplineFuncTest *sft = ncm_spline_func_test_new ();

  ncm_spline_func_test_set_params_info_all (sft, npar, -10.0, 10.0);

  ncm_spline_func_test_set_scale (sft, 1.0);

  ncm_spline_func_test_set_ngrid (sft, 10000);

  ncm_spline_func_test_prepare (sft, NCM_SPLINE_FUNCTION_SPLINE, NCM_SPLINE_FUNC_TEST_TYPE_PDF_FLAT);

  ncm_spline_func_test_monte_carlo_and_save_to_txt (sft, nsim, "mc.txt");

  ncm_spline_func_test_log_vals_mc_stats (sft);

  ncm_spline_func_test_unref (sft);

  return 0;
}

The function ncm_spline_func_test_monte_carlo_and_save_to_txt() creats the ascii file “mc.txt”. It saves the same statistics as printed by the function ncm_spline_func_test_log_vals_one_grid_stats () for each realization. Therefore it is going to be quite a big file. In this example it has 165 megabytes. That is the reason why it is not allowed to save all the grids informations and also because the file is filled on the fly.

Below is shown the output of ncm_spline_func_test_log_vals_mc_stats():

Monte Carlo statistics for 1000000 simulations

• NcmSplineFunc number of knots: 750.58 +/- 114.83

• Ncm clean grid (%): 97.59

• Lin clean grid (%): 27.78

• Ncm outliers (%): 0.03 +/- 0.24

• Lin outliers (%): 0.25 +/- 0.56

• Ncm diff. : 1.624e-07 +/- 6.353e-06

• Lin diff. : 7.880e-09 +/- 2.385e-07

• Ncm abs. max. diff. : 3.891e-05 +/- 9.496e-03

• Lin abs. max. diff. : 5.886e-06 +/- 1.099e-03

First, we need to define two different statistics, “one grid” and “Monte Carlo”. The mean value for each will be represented by $\overline{X}$ and $\langle X \rangle$, respectively. The standard deviation $\sigma$ is applied to the Monte Carlo realizations.

• “NcmSplineFunc number of knots:” $\langle \mathrm{spline \ length} \rangle \pm \sigma \left( \mathrm{spline \ length} \right)$.
• “clean grid (%):” it is the proportion of the grids realizations with no outliers at all.
• “outliers (%):” $\Big \langle \mathrm{outlier} \Big \rangle \pm \sigma \left( \mathrm{outlier} \right)$.
• “diff.:” $\Big \langle \overline{\Delta f(x)} \Big \rangle \pm \sigma \left( \overline{\Delta f(x)} \right)$.
• “abs. max. diff.:” $\Big \langle |\Delta f(x)|{\mathrm{max}} \Big \rangle \pm \sigma \left( |\Delta f(x)|{\mathrm{max}} \right)$.

Note that that the “clean grid” shows a much better result for NcmSplineFunc over the linear grid. Also “outliers” shows the same trend but not at the same level. Nevertheless, the “diff.” and “abs. max. diff.” seems to contradict both of them. To check this apparent contradiction, it is needed to look into the “mc.txt” file, which should show that the ~2% “non-clean” NcmSplineFunc $\hat{f}(x)$ should have outliers with high values, probably indicating some issue with those functions.

Example: user supplied function with 100 thousand realizations.

#include <numcosmo/numcosmo.h>

gdouble
f_user (gdouble x, gpointer pin)
{
  NcmSplineFuncTest *sft = NCM_SPLINE_FUNC_TEST (pin);

  NcmVector *params = ncm_spline_func_test_peek_current_params (sft);

  gdouble p1 = ncm_vector_fast_get (params, 0);
  gdouble p2 = ncm_vector_fast_get (params, 1);

  return exp (sin (p1 * x)) * sin (x) * sin (x) + p2;
}

int
main (void)
{
  gsl_function F;

  NcmSplineFuncTest *sft = ncm_spline_func_test_new ();

  ncm_spline_func_test_set_type (sft, NCM_SPLINE_FUNC_TEST_TYPE_USER);

  ncm_spline_func_test_set_rel_error (sft, 1.e-10);
  ncm_spline_func_test_set_scale (sft, 10.0);
  ncm_spline_func_test_set_out_threshold (sft, 10.0);

  ncm_spline_func_test_set_xi (sft, -0.5);
  ncm_spline_func_test_set_xf (sft, +0.5);

  ncm_spline_func_test_set_params_info_all (sft, 2, 0., 10.);

  F.function = &f_user;
  F.params   = sft;

  ncm_spline_func_test_set_user_gsl_function (sft, &F);

  ncm_spline_func_test_prepare (sft, NCM_SPLINE_FUNCTION_SPLINE, NCM_SPLINE_FUNC_TEST_TYPE_PDF_NORMAL);

  ncm_spline_func_test_monte_carlo (sft, 100000);

  ncm_spline_func_test_log_vals_mc_stats (sft);

  ncm_spline_func_test_unref (sft);

  return 0;
}

The user function “f_user” has two parameters, $p_1$ and $p_2$, and is given by \begin{equation} f(x) = \exp \left[ \sin(p_1 \, x) \right] \sin^{2}(x) + p_2 \,\, . \end{equation} Unlike the previous example, the Monte Carlo statistic is not saved in a file. Instead, it is used the ncm_spline_func_test_monte_carlo() function. This procedure has a faster execution time, but the drawback is information lost. The function ncm_spline_func_test_set_out_threshold() saves information only for the functions with outilers above 10\% of the threshold given in Eq. \eqref{eq:condition}.

Below is shown the output of ncm_spline_func_test_log_vals_mc_stats():

Monte Carlo statistics for 100000 simulations

• NcmSplineFunc number of knots: 581.23 +/- 331.82

• Ncm clean grid (%): 79.74

• Lin clean grid (%): 91.69

• Ncm outliers (%): 0.24 +/- 0.68

• Lin outliers (%): 0.01 +/- 0.03

• Ncm diff. : -8.430e-10 +/- 6.675e-09

• Lin diff. : -1.745e-12 +/- 5.589e-11

• Ncm abs. max. diff. : 3.114e-08 +/- 8.864e-06

• Lin abs. max. diff. : 8.927e-10 +/- 6.938e-10

The results shows a better approximation by the linear grid compared to the NcmSplineFunc method. In order to try to understand this behaviour, the user should look into the files created by ncm_spline_func_test_set_out_threshold(), but note that “Ncm abs. max. diff.” has mean near the desired tolerance. In this case we have an absolute error $\left( |f(x)| + 10 \right) \mathrm{rel \_ error} \approx 10^{-9} \rightarrow 10^{-8}$, which is around the “Ncm abs. max. diff.”. Around 10 grids have outliers greater than 10\%.

If, instead, it is applied the method #NCM_SPLINE_FUNCTION_4POINTS, as was found previously, the result is more robust compared to #NCM_SPLINE_FUNCTION_SPLINE. Applying the 4POINTS, we got no grids with outliers, as can be seem below:

Monte Carlo statistics for 100000 simulations

• NcmSplineFunc number of knots: 4173.34 +/- 2363.36

• Ncm clean grid (%): 100.00

• Lin clean grid (%): 100.00

• Ncm outliers (%): 0.00 +/- 0.00

• Lin outliers (%): 0.00 +/- 0.00

• Ncm diff. : 2.707e-14 +/- 3.753e-13

• Lin diff. : -9.547e-16 +/- 1.582e-14

• Ncm abs. max. diff. : 4.948e-12 +/- 7.645e-12

• Lin abs. max. diff. : 9.834e-14 +/- 7.774e-14

It is highly advisable to apply a scale for this method. Its high precision can produce a crash while crossing the abscissa.

Note that both Monte Carlo examples can have slightly different results for any given run, due to the fact that the seed is created internally.

The examples above are stress tests, with some not so commom functions, but even then, the results from NcmSplineFunc are quite remarkable. .