Skip to content

Commit

Permalink
Introduce an optional non linear mapping before normalization (#8)
Browse files Browse the repository at this point in the history
  • Loading branch information
@jlelong
jlelong authored Dec 20, 2023
1 parent 33c7fbf commit 6efb848
Show file tree
Hide file tree
Showing 4 changed files with 292 additions and 180 deletions.
78 changes: 73 additions & 5 deletions examples/basis_test.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include "pnl/pnl_basis.h"
#include "pnl/pnl_mathtools.h"
#include "pnl/pnl_random.h"
#include "pnl/pnl_cdf.h"

#define DATA_DIR "Data_basis"
#include "tests_utils.h"
Expand Down Expand Up @@ -200,12 +201,22 @@ static double function2d_3(double *x)
return (x[0] * x[1]);
}

static double function2d_4(double *x)
{
return exp(1 + 2 * x[0] + x[1]);
}

static double myfunction(const PnlVect *x, void *p)
{
double a = ((double *) p)[0];
return log(1 + a * SQR(pnl_vect_norm_two(x)));
}

static double map(double x, int _dim, void *_p)
{
return exp(x);
}

enum {ABS_ERR = 1, MEAN_ERR = 2};
static double regression_multid_aux(double(*f)(double *), PnlBasis *basis, int err_type)
{
Expand Down Expand Up @@ -246,8 +257,7 @@ static double regression_multid_aux(double(*f)(double *), PnlBasis *basis, int e
mean = err = 0.;
for (i = 0; i < n; i++)
{
double tmp = (*f)(pnl_mat_lget(t, i, 0)) -
pnl_basis_eval(basis, alpha, pnl_mat_lget(t, i, 0));
double tmp = (*f)(pnl_mat_lget(t, i, 0)) - pnl_basis_eval(basis, alpha, pnl_mat_lget(t, i, 0));
mean += tmp * tmp;
if (fabs(tmp) > err) err = fabs(tmp);
}
Expand Down Expand Up @@ -320,6 +330,18 @@ static void regression_multid()
pnl_test_eq_abs(err, 0., 1E-5, "pnl_basis_eval (extra functions)", "log (1+x[0]*x[0] + x[1]*x[1]) + (x[0]*x[0] + x[1]*x[1]) on [-2,2]^2");
pnl_basis_free(&basis);
}

{
double err;
int basis_name = PNL_BASIS_HERMITE;
int nb_variates = 2; /* functions with values in R^2 */
int degree = 3; /* total product degree */
PnlBasis *basis = pnl_basis_create_from_prod_degree(basis_name, degree, nb_variates);
pnl_basis_set_map(basis, map, NULL, NULL, NULL, 0);
err = regression_multid_aux(function2d_4, basis, ABS_ERR);
pnl_test_eq_abs(err, 0., 1E-5, "pnl_basis_eval (map)", "exp(1 + 2 * x[0] + x[1]) on [-2,2]^2");
pnl_basis_free(&basis);
}
}

static double fonction_a_retrouver(double t, double x)
Expand Down Expand Up @@ -411,7 +433,7 @@ static void pnl_basis_eval_test()
pnl_rng_sseed(rng, 0);

/*
* Random points where the function will be evaluted
* Random points where the function will be evaluated
*/
pnl_vect_rng_uni(x, n, -5, 4, rng);
pnl_vect_rng_uni(t, n, 0, 1, rng);
Expand Down Expand Up @@ -524,11 +546,30 @@ static void pnl_basis_eval_test()
pnl_mat_free(&X);
}

static double piecewise_constant1d(double x, double a, double b, int n)
{
if (x < a || x > b) return 0.;
int i = (x - a) / (b - a) * n;
return (double) i;
}


static double piecewise_constant2d(double x1, double x2, double a1, double b1, double a2, double b2, int n)
{
return piecewise_constant1d(pnl_cdfnor(x1), a1, b1, n) + piecewise_constant1d(pnl_cdfnor(x2), a2, b2, n);
}

static double map_local(double x, int _dim, void *_p)
{
return pnl_cdfnor(x);
}

static void local_basis_test()
{
int i;
int dim = 2;
PnlBasis *B = pnl_basis_local_create_regular(30, dim);
int n_intervals = 30;
PnlBasis *B = pnl_basis_local_create_regular(n_intervals, dim);
int n = 100000;
PnlRng *rng = pnl_rng_create(PNL_RNG_MERSENNE);
PnlMat *x = pnl_mat_new();
Expand All @@ -540,7 +581,34 @@ static void local_basis_test()
pnl_basis_fit_ls(B, alpha, x, y);
for (i = 0; i < alpha->size; i++)
{
pnl_test_eq_abs(GET(alpha, i), 1, 1E-6, "local_constant_regression", "");
if (pnl_isequal_abs(GET(alpha, i), 1, 1E-6) == PNL_FALSE)
{
pnl_test_set_fail( "local_constant_regression", GET(alpha, i), 1);
return;
}
}
if (i == alpha->size)
{
pnl_test_set_ok("local_constant_regression");
}
for (i = 0; i < x->m; i++)
{
LET(y, i) = piecewise_constant2d(MGET(x, i, 0),MGET(x, i, 1), -1, 1, -1, 1, n_intervals);
}
pnl_basis_set_map(B, map_local, NULL, NULL, NULL, 0);
pnl_basis_fit_ls(B, alpha, x, y);
for (i = 0; i < x->m; i++)
{
double evalB = pnl_basis_eval(B, alpha, pnl_mat_lget(x, i, 0));
if (pnl_isequal_abs(evalB, GET(y, i), 1E-6) == PNL_FALSE)
{
pnl_test_set_fail( "local_cdfnor_regression", evalB, GET(y, i));
return;
}
}
if (i == alpha->size)
{
pnl_test_set_ok("local_constant_regression");
}
pnl_basis_free(&B);
pnl_vect_free(&y);
Expand Down
53 changes: 37 additions & 16 deletions man/bases.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,17 @@ \subsection{Overview}
/** Extra parameters to pass to basis functions */
void *f_params;
/** Size of params in bytes to be passed to malloc */
size_t f_params_size;
size_t f_params_size;
/** Non linear mapping of the data */
double (*map)(double x, int dim, void *params);
/** First derivative of the non linear mapping */
double (*Dmap)(double x, int dim, void *params);
/** Second dérivate of the linear mapping */
double (*D2map)(double x, int dim, void *params);
/** Extra paramaters for map, Dmap and D2map */
void *map_params;
/** Size of @p map_params in bytes to be passed to malloc */
size_t map_params_size;
};
\end{lstlisting}

Expand Down Expand Up @@ -84,14 +94,14 @@ \subsubsection{Tensor functions}
\item \describefun{\PnlBasis *}{pnl_basis_create}{int index, int
nb_func, int nb_variates}
\sshortdescribe Create a \PnlBasis for the family
defined by \var{index} (see Table~\ref{basis_index} and
defined by \var{index} (see Table~\ref{tab:basis_index} and
\reffun{pnl_basis_type_register}) with \var{nb_variates}
variates. The basis will contain \var{nb_func}.

\item \describefun{\PnlBasis *}{pnl_basis_create_from_degree}{int
index, int degree, int nb_variates}
\sshortdescribe Create a \PnlBasis for the family
defined by \var{index} (see Table~\ref{basis_index} and \reffun{pnl_basis_type_register}) with total degree less
defined by \var{index} (see Table~\ref{tab:basis_index} and \reffun{pnl_basis_type_register}) with total degree less
or equal than \var{degree} and \var{nb_variates} variates. The total degree is
the sum of the partial degrees.\\
For instance, calling \verb!pnl_basis_create_from_degree(index, 2, 4)! is
Expand All @@ -118,7 +128,7 @@ \subsubsection{Tensor functions}
\right) \]
\item \describefun{\PnlBasis *}{pnl_basis_create_from_prod_degree}{int index, int degree, int nb_variates}
\sshortdescribe Create a \PnlBasis for the family
defined by \var{index} (see Table~\ref{basis_index} and \reffun{pnl_basis_type_register}) with total degree less
defined by \var{index} (see Table~\ref{tab:basis_index} and \reffun{pnl_basis_type_register}) with total degree less
or equal than \var{degree} and \var{nb_variates} variates. The total degree is
the product of \var{MAX(1, d_i)} where the \var{d_i} are the partial degrees.

Expand All @@ -133,7 +143,7 @@ \subsubsection{Tensor functions}
\item \describefun{\PnlBasis *}{pnl_basis_create_from_tensor}{int
index, PnlMatInt \ptr T}
\sshortdescribe Create a \PnlBasis for the polynomial family
defined by \var{index} (see Table~\ref{basis_index}) using the basis
defined by \var{index} (see Table~\ref{tab:basis_index}) using the basis
described by the tensor matrix \var{T}. The number of lines of \var{T} is
the number of functions of the basis whereas the numbers of columns of
\var{T} is the number of variates of the functions.
Expand All @@ -160,7 +170,7 @@ \subsubsection{Tensor functions}
\item \describefun{void }{pnl_basis_set_from_tensor}{\PnlBasis \ptr
b, int index, const \PnlMatInt \ptr T}
\sshortdescribe Set an alredy existing basis \var{b} to a polynomial family
defined by \var{index} (see Table~\ref{basis_index}) using the basis
defined by \var{index} (see Table~\ref{tab:basis_index}) using the basis
described by the tensor matrix \var{T}. The number of lines of \var{T} is
the number of functions of the basis whereas the numbers of columns of
\var{T} is the number of variates of the functions. \\
Expand All @@ -185,7 +195,7 @@ \subsubsection{Local basis functions}
\sshortdescribe Return the linear index of the cell containing \var{x}. It is an integer between $0$ and $(\prod_{i=1}^d n_i) - 1$.
\end{itemize}

If the domain you want to consider is not ${[-1,1]}^d$, use the functions \reffun{pnl_basis_set_domain} or \reffun{pnl_basis_set_reduced} to map your product space into ${[-1,1]}^d$.
If the domain you want to consider is not ${[-1,1]}^d$, use the functions \reffun{pnl_basis_set_domain}, \reffun{pnl_basis_set_reduced} and \reffun{pnl_basis_set_map} to map your product space into ${[-1,1]}^d$.

\subsubsection{Standard multivariate functions}

Expand Down Expand Up @@ -254,7 +264,7 @@ \subsection{Functions}
\end{itemize}


Functional regression based on a least square approach often leads to ill conditioned linear systems. One way of improving the stability of the system is to use centered and renormalised polynomials so that the original domain of interest $\cD$ (a subset of $\R^d$) is mapped to $[-1,1]^d$. If the domain $\cD$ is rectangular and writes $[a, b]$ where $a,b \in \R^d$, the mapping is done by
Functional regression based on a least square approach often leads to ill conditioned linear systems. One way of improving the stability of the system is to use centered and renormalized polynomials so that the original domain of interest $\cD$ (a subset of $\R^d$) is mapped to $[-1,1]^d$. If the domain $\cD$ is rectangular and writes $[a, b]$ where $a,b \in \R^d$, the \emph{reduction} mapping is done by
\begin{equation}
\label{basis_reduced}
x \in \cD \longmapsto \left(\frac{x_i - (b_i+a_i)/2}{(b_i - a_i)/2}
Expand All @@ -263,26 +273,37 @@ \subsection{Functions}
\begin{itemize}
\item \describefun{void}{pnl_basis_set_domain}{\PnlBasis \ptr B,
const \PnlVect \ptr a, const \PnlVect \ptr b}
\sshortdescribe This function declares \var{B} as a centered and normalised basis
\sshortdescribe This function declares \var{B} as a centered and normalized basis
as defined by Equation~\ref{basis_reduced}. Calling this function is equivalent to
calling \reffun{pnl_basis_set_reduced} with \var{center=(b+a)/2} and
\var{scale=(b-a)/2}.
\item \describefun{void}{pnl_basis_set_reduced}{\PnlBasis \ptr B,
const \PnlVect \ptr center, const \PnlVect \ptr scale}
\sshortdescribe This function declares \var{B} as a centered and normalised basis
using the mapping
\sshortdescribe This function declares \var{B} as a centered and normalized basis using the mapping
\begin{equation*}
x \in \cD \longmapsto \left(\frac{x_i - \var{center}_i }{\var{scale}_i}
\right)_{i=1,\cdots,d}
\right)_{i=1,\dots,d}
\end{equation*}
\item \describefun{void}{pnl_basis_reset_reduced}{\PnlBasis \ptr B}
\sshortdescribe Reset the reduction settings.
\end{itemize}
Note that this renormlization does not apply to the extra functions by \reffun{pnl_basis_add_function} but only to the functions defined by the tensor \var{T}.
Note that this renormalization does not apply to the extra functions by \reffun{pnl_basis_add_function} but only to the functions defined by the tensor \var{T}.

It is also possible to apply a non linear map $\varphi_i:\R \to \R$ to the $i-$th coordinate of the input variable before the reduction operation if any. Then, the input variables are transformed according to
\begin{equation}
\label{basis_non_linear_reduced}
x \in \cD \longmapsto \left(\frac{\varphi_i(x_i) - (b_i+a_i)/2}{(b_i - a_i)/2}
\right)_{i=1,\dots,d}
\end{equation}


\begin{itemize}
\item \describefun{void}{pnl_basis_set_map}{\PnlBasis \ptr B, double (\ptr map)(double x, int i, void\ptr params), double (\ptr Dmap)(double x, int i, void\ptr params), double (\ptr D2map)(double x, int i, void\ptr params), void \ptr params, size_t size_params}
\sshortdescribe Define the non linear applied to the input variables before the reduction operation. The parameters \var{Dmap} and \var{D2map} can be \var{NULL}. The three functions \var{map}, \var{Dmap} and \var{D2map} take as second argument the index of the coordinate, while their third argument is used to passe extra parameters defined by \var{params}.
\end{itemize}

\begin{itemize}
\item \describefun{int}{pnl_basis_fit_ls}{\PnlBasis \ptr
P, \PnlVect \ptr coef, \PnlMat \ptr x,
\PnlVect \ptr y}
\item \describefun{int}{pnl_basis_fit_ls}{\PnlBasis \ptr P, \PnlVect \ptr coef, \PnlMat \ptr x, \PnlVect \ptr y}
\sshortdescribe Compute the coefficients \var{coef} defined by
\begin{equation*}
\var{coef} = \arg\min_\alpha \sum_{i=1}^n
Expand Down
15 changes: 13 additions & 2 deletions src/include/pnl/pnl_basis.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,8 +85,18 @@ struct _PnlBasis
int len_func_list;
/** Extra parameters to pass to basis functions */
void *f_params;
/** Size of params in bytes to be passed to malloc */
size_t f_params_size;
/** Size of @p f_params in bytes to be passed to malloc */
size_t f_params_size;
/** Non linear mapping of the data */
double (*map)(double x, int dim, void *params);
/** First derivative of the non linear mapping */
double (*Dmap)(double x, int dim, void *params);
/** Second dérivate of the linear mapping */
double (*D2map)(double x, int dim, void *params);
/** Extra paramaters for map, Dmap and D2map */
void *map_params;
/** Size of @p map_params in bytes to be passed to malloc */
size_t map_params_size;
};

extern int pnl_basis_type_register(const char *name, double (*f)(double, int, int, void*), double (*Df)(double, int, int, void*), double (*D2f)(double, int, int, void*), int is_orthogonal);
Expand All @@ -109,6 +119,7 @@ extern int pnl_basis_local_get_index(const PnlBasis *basis, const double *x);
extern void pnl_basis_set_domain(PnlBasis *B, const PnlVect *xmin, const PnlVect *xmax);
extern void pnl_basis_set_reduced(PnlBasis *B, const PnlVect *center, const PnlVect *scale);
extern void pnl_basis_reset_reduced(PnlBasis *B);
extern void pnl_basis_set_map(PnlBasis *B, double (*map)(double, int, void*), double (*Dmap)(double, int, void*), double (*D2map)(double, int, void*), void *params, size_t size_params);
extern void pnl_basis_free(PnlBasis **basis);
extern void pnl_basis_print(const PnlBasis *B);
extern int pnl_basis_fit_ls(const PnlBasis *f, PnlVect *coef, const PnlMat *x, const PnlVect *y);
Expand Down
Loading

0 comments on commit 6efb848

Please sign in to comment.