Remove the full tensor basis representation

man/bases.tex

@@ -20,11 +20,9 @@ \subsection{Overview}
   int           nb_variates;
   /** The total number of elements in the basis */
   int           nb_func;
-  /** The tensor matrix */
-  PnlMatInt    *T;
   /** The sparse Tensor matrix */
   PnlSpMatInt  *SpT;
-  /** The number of functions in the tensor #T */
+  /** The number of functions in the tensor @p SpT */
   int           len_T;
   /** Compute the i-th element of the one dimensional basis.  As a convention, (*f)(x, 0) MUST be equal to 1 */
   double      (*f)(double x, int i, int dim, void *params);
@@ -40,7 +38,7 @@ \subsection{Overview}
   double       *scale;
   /** An array of additional functions */
   PnlRnFuncR   *func_list;
-  /** The number of functions in #func_list */
+  /** The number of functions in @p func_list */
   int           len_func_list;
   /** Extra parameters to pass to basis functions */
   void         *f_params;
@@ -50,9 +48,9 @@ \subsection{Overview}
   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 */
+  /** Second derivate of the linear mapping */
   double       (*D2map)(double x, int dim, void *params);
-  /** Extra paramaters for map, Dmap and D2map */
+  /** Extra parameters 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;
@@ -102,30 +100,7 @@ \subsubsection{Tensor functions}
       index, int degree, int nb_variates}
     \sshortdescribe Create a \PnlBasis for the family
     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
-    equivalent to calling \verb!pnl_basis_create_from_tensor(index, T)! where
-    \var{T} is given by
-    \[ \left(
-      \begin{array}{cccc}
-        0 & 0 & 0 & 0\\
-        1 & 0 & 0 & 0\\
-        0 & 1 & 0 & 0\\
-        0 & 0 & 1 & 0\\
-        0 & 0 & 0 & 1\\
-        1 & 1 & 0 & 0\\
-        1 & 0 & 1 & 0\\
-        1 & 0 & 0 & 1\\
-        0 & 1 & 1 & 0\\
-        0 & 1 & 0 & 1\\
-        0 & 0 & 1 & 1\\
-        2 & 0 & 0 & 0\\
-        0 & 2 & 0 & 0\\
-        0 & 0 & 2 & 0\\
-        0 & 0 & 0 & 2\\
-      \end{array}
-    \right) \]
+    or equal than \var{degree} and \var{nb_variates} variates. The total degree is the sum of the partial degrees.
   \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{tab:basis_index} and \reffun{pnl_basis_type_register}) with total degree less
@@ -147,9 +122,7 @@ \subsubsection{Tensor functions}
     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.
-    Note that \var{T} is not copied inside this function but only its address is
-    stored, so {\bf never} free \var{T}. It will be freed when calling
-    \reffun{pnl_basis_free} on the returned object. i\\
+    \\
     Here is an example of a tensor matrix. Assume you are working with three
     variate functions, the basis \verb!{ 1, x, y, z, x^2, xy, yz, z^3}! is
     decomposed in the one dimensional canonical basis using the following tensor
@@ -167,6 +140,12 @@ \subsubsection{Tensor functions}
       \end{array}
     \right) \]
 
+  \item \describefun{\PnlBasis *}{pnl_basis_create_from_sparse_tensor}{int
+      index, PnlSpMatInt \ptr SpT}
+    \sshortdescribe Create a \PnlBasis for the polynomial family
+    defined by \var{index} (see Table~\ref{tab:basis_index}) using the basis
+    described by the sparse tensor matrix \var{SpT}. The number of lines of \var{SpT} is the number of functions of the basis whereas the numbers of columns of \var{SpT} is the number of variates of the 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
@@ -182,7 +161,7 @@ \subsubsection{Tensor functions}
 
 \subsubsection{Local basis functions}
 
-A local basis is a family of indicator functions of a Cartesian partition of ${[-1,1]}^d$. Let $(n_i)_{i \in \{1, \dots, d\}}$ be the number of interval along each direction. An element of the partition write $A_k = \prod_{i=1}^d [-1 + k_i/n_i, -1 + (k_i + 1) / n_i]$ where $k$ is the multi-index defined by $k_i \in \{0, \dots, n_i - 1\}$ for all $i\in \{1,\dots, d\}$. These functions are orthogonal; this property is used by \reffun{pnl_basis_fit_ls}. Note that they are not differentiable.
+A local basis is a family of indicator functions of a Cartesian partition of ${[-1,1]}^d$. Let $(n_i)_{i \in \{1, \dots, d\}}$ be the number of interval along each direction. An element of the partition write \[A_k = \prod_{i=1}^d \left[-1 + \frac{k_i}{n_i}, -1 + \frac{k_i + 1}{n_i}\right]\] where $k$ is the multi-index defined by $k_i \in \{0, \dots, n_i - 1\}$ for all $i\in \{1,\dots, d\}$. These functions are orthogonal for the standard $L^2$ scalar product; this property is used by \reffun{pnl_basis_fit_ls}. Note that they are not differentiable.
 
 We provide the helper function \reffun{pnl_basis_local_get_index} to compute the linear representation of the multi-index $k$.
 

man/pnl.css

@@ -25,10 +25,14 @@ h5 { font-size: 110%;   font-style: bold; }
 
 
 /* Font settings for page elements */
-code, pre, tt,  {
+code, pre, tt {
     font-family:  monospace;
 }
 
+pre.lstlisting, span.verb {
+    font-size: 80%;
+}
+
 /* Indented block */
 div.indent {  margin-left: 50px;}
 

src/include/pnl/pnl_basis.h

@@ -62,10 +62,8 @@ struct _PnlBasis
   /** The total number of elements in the basis */
   int           nb_func;
   /** The tensor matrix */
-  PnlMatInt    *T;
-  /** The sparse Tensor matrix */
   PnlSpMatInt  *SpT;
-  /** The number of functions in the tensor #T */
+  /** The number of functions in the tensor @p SpT */
   int           len_T;
   /** Compute the i-th element of the one dimensional basis.  As a convention, (*f)(x, 0) MUST be equal to 1 */
   double      (*f)(double x, int i, int dim, void *params);
@@ -81,7 +79,7 @@ struct _PnlBasis
   double       *scale;
   /** An array of additional functions */
   PnlRnFuncR   *func_list;
-  /** The number of functions in #func_list */
+  /** The number of functions in @p func_list */
   int           len_func_list;
   /** Extra parameters to pass to basis functions */
   void         *f_params;
@@ -91,9 +89,9 @@ struct _PnlBasis
   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 */
+  /** Second derivate of the linear mapping */
   double       (*D2map)(double x, int dim, void *params);
-  /** Extra paramaters for map, Dmap and D2map */
+  /** Extra parameters 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;
@@ -111,6 +109,7 @@ extern void pnl_basis_clone(PnlBasis *dest, const PnlBasis *src);
 extern PnlBasis* pnl_basis_copy(const PnlBasis *B);
 extern void pnl_basis_set_type(PnlBasis *B, int index);
 extern void  pnl_basis_set_from_tensor(PnlBasis *b, const PnlMatInt *T);
+extern void  pnl_basis_set_from_sparse_tensor(PnlBasis *b, const PnlSpMatInt *T);
 extern PnlBasis* pnl_basis_create_from_tensor( int index, const PnlMatInt *T);
 extern void pnl_basis_del_elt(PnlBasis *B, const PnlVectInt *d);
 extern void pnl_basis_del_elt_i(PnlBasis *B, int i);

src/interpol/basis.c

@@ -227,48 +227,6 @@ static PnlMatInt *compute_tensor_from_prod_degree(int degree, int nb_variates)
   return compute_tensor_from_degree_function(degree, nb_variates, nb_elements, NULL, count_prod_degree);
 }
 
-/**
- * Compute the tensor for a tensor local basis
- *
- * @param space_dim the dimension of the state space
- * @param n_intervals this is an array of size \a space_dim describing the number of intervals for every dimension
- * @return PnlMatInt
- */
-static PnlMatInt* compute_tensor_permutation(int space_dim, int *n_elements)
-{
-  int permutation_length, i, col;
-  int inner_loop_length, outer_loop_length;
-  PnlMatInt *T;
-  permutation_length = 1;
-  for (i = 0; i < space_dim; i++)
-    {
-      permutation_length *= n_elements[i];
-    }
-  T = pnl_mat_int_create(permutation_length, space_dim);
-  outer_loop_length = permutation_length;
-  inner_loop_length = 1;
-  for (col = 0; col < space_dim; col++)
-    {
-      int outer_index;
-      outer_loop_length /= n_elements[col];
-      for (outer_index = 0; outer_index < outer_loop_length; outer_index++)
-        {
-          int middle_loop;
-          for (middle_loop = 0; middle_loop < n_elements[col]; middle_loop++)
-            {
-              int inner_loop;
-              for (inner_loop = 0; inner_loop < inner_loop_length; inner_loop++)
-                {
-                  int row = outer_index * inner_loop_length * n_elements[col] + middle_loop * inner_loop_length + inner_loop;
-                  PNL_MLET(T, row, col) = middle_loop + 1;
-                }
-            }
-        }
-      inner_loop_length *= n_elements[col];
-    }
-  return T;
-}
-
 /**
  *  Canonical polynomials
  *  @param x the address of a real number
@@ -829,7 +787,6 @@ PnlBasis *pnl_basis_new()
   o->id = 0;
   o->label = "";
   o->nb_variates = 0;
-  o->T = NULL;
   o->SpT = NULL;
   o->isreduced = 0;
   o->center = NULL;
@@ -854,25 +811,35 @@ PnlBasis *pnl_basis_new()
 }
 
 /**
- * Create a PnlBasis and stores it into its first argument
+ * Set an existing basis to a given full tensor representation
  *
  * @param b an already allocated basis (as returned by pnl_basis_new for instance)
- * @param T the tensor of the multi-dimensional basis. No copy of T is done, so
- * do not free T. It will be freed transparently by pnl_basis_free
+ * @param T the full tensor of the multi-dimensional basis.
  */
 void pnl_basis_set_from_tensor(PnlBasis *b, const PnlMatInt *T)
 {
-  b->nb_func = T->m;
-  b->nb_variates = T->n;
+  PnlSpMatInt *SpT = pnl_sp_mat_int_create_from_mat(T);
+  pnl_basis_set_from_sparse_tensor(b, SpT);
+}
+
+/**
+ * Set an existing basis to a given sparse tensor representation
+ *
+ * @param b an already allocated basis (as returned by pnl_basis_new for instance)
+ * @param SpT the sparse tensor of the multi-dimensional basis. No copy of SpT is done, so
+ * do not free it. It will be freed transparently by pnl_basis_free
+ */
+void pnl_basis_set_from_sparse_tensor(PnlBasis *b, const PnlSpMatInt *SpT)
+{
+  b->nb_func = SpT->m;
+  b->nb_variates = SpT->n;
   if(b->func_list) free(b->func_list);
   b->len_func_list = 0;
-  b->len_T = T->m;
+  b->len_T = SpT->m;
 
   /* Not sure this is the right place to put it */
-  pnl_mat_int_free(&(b->T));
   pnl_sp_mat_int_free(&(b->SpT));
-  b->T = (PnlMatInt *) T;
-  b->SpT = pnl_sp_mat_int_create_from_mat(T);
+  b->SpT = (PnlSpMatInt *) SpT;
 }
 
 /**
@@ -1057,7 +1024,7 @@ void pnl_basis_add_function(PnlBasis *b, PnlRnFuncR *f)
  */
 void pnl_basis_clone(PnlBasis *dest, const PnlBasis *src)
 {
-  pnl_basis_set_from_tensor(dest, src->T);
+  pnl_basis_set_from_sparse_tensor(dest, src->SpT);
   dest->id = src->id;
   dest->label = src->label;
   dest->f = src->f;
@@ -1091,7 +1058,6 @@ void pnl_basis_clone(PnlBasis *dest, const PnlBasis *src)
 void pnl_basis_del_elt_i(PnlBasis *B, int i)
 {
   PNL_CHECK(i >= B->len_T, "size mismatch",  "basis_del_elt_i");
-  pnl_mat_int_del_row(B->T, i);
   pnl_sp_mat_int_del_row(B->SpT, i);
   --(B->nb_func);
   --(B->len_T);
@@ -1109,8 +1075,17 @@ void pnl_basis_del_elt(PnlBasis *B, const PnlVectInt *d)
   PNL_CHECK(B->nb_variates != d->size, "size mismatch",  "basis_del_elt");
   for (i = 0 ; i < B->nb_func ; i++)
     {
-      PnlVectInt Ti = pnl_vect_int_wrap_mat_row(B->T, i);
-      if (pnl_vect_int_eq(&Ti, d) == PNL_TRUE) break;
+      int j;
+      int test = 1;
+      for (j = 0; j < d->size; j++)
+        {
+          if (pnl_sp_mat_int_get(B->SpT, i, j) != pnl_vect_int_get(d, j))
+            {
+              test = 0;
+              break;
+            }
+        }
+      if (test) break;
     }
   if (i < B->nb_func) pnl_basis_del_elt_i(B, i);
 }
@@ -1124,7 +1099,6 @@ void pnl_basis_del_elt(PnlBasis *B, const PnlVectInt *d)
 void pnl_basis_add_elt(PnlBasis *B, const PnlVectInt *d)
 {
   PNL_CHECK(B->nb_variates != d->size, "size mismatch",  "basis_del_elt");
-  pnl_mat_int_add_row(B->T, B->T->m, d);
   pnl_sp_mat_int_add_row(B->SpT, B->SpT->m, d);
   ++(B->nb_func);
   ++(B->len_T);
@@ -1246,7 +1220,6 @@ void pnl_basis_set_map(
 void pnl_basis_free(PnlBasis **B)
 {
   if (*B == NULL) return;
-  pnl_mat_int_free(&((*B)->T));
   pnl_sp_mat_int_free(&(*B)->SpT);
   if ((*B)->f_params_size > 0)
     {
@@ -1326,7 +1299,7 @@ void pnl_basis_print(const PnlBasis *B)
  */
 double pnl_basis_ik(const PnlBasis *b, const double *x, int i, int k)
 {
-  const int Tik = PNL_MGET(b->T, i, k);
+  const int Tik = pnl_sp_mat_int_get(b->SpT, i, k);
   if (Tik == 0) return 1.;
   return f_reduction_map(b, x[k], Tik, k);
 }
@@ -1380,7 +1353,7 @@ double pnl_basis_i_D(const PnlBasis *b, const double *x, int i, int j)
 
   CHECK_IS_DIFFERENTIABLE(b);
   /* Test if the partial degree is small enough to return 0 */
-  if (PNL_MGET(b->T, i, j) == 0) return 0.;
+  if (pnl_sp_mat_int_get(b->SpT, i, j) == 0) return 0.;
 
   for (l = b->SpT->I[i] ; l < b->SpT->I[i + 1] ; l++)
     {
@@ -1413,8 +1386,8 @@ double pnl_basis_i_D2(const PnlBasis *b, const double *x, int i, int j1, int j2)
   double aux = 1;
   CHECK_IS_DIFFERENTIABLE(b);
   /* Test if the partial degree is small enough to return 0 */
-  if ((j1 == j2) && PNL_MGET(b->T, i, j1) <= 1) return 0.;
-  if (PNL_MGET(b->T, i, j1) == 0 || PNL_MGET(b->T, i, j2) == 0) return 0.;
+  if ((j1 == j2) && pnl_sp_mat_int_get(b->SpT, i, j1) <= 1) return 0.;
+  if (pnl_sp_mat_int_get(b->SpT, i, j1) == 0 || pnl_sp_mat_int_get(b->SpT, i, j2) == 0) return 0.;
 
   if (j1 == j2)
     {
@@ -1655,7 +1628,7 @@ void pnl_basis_eval_derivs(const PnlBasis *b, const PnlVect *coef, const double
               const int k = b->SpT->J[l];
               if (k != j) auxf *= f[k];
             }
-            const int Tij = PNL_MGET(b->T, i, j);
+            const int Tij = pnl_sp_mat_int_get(b->SpT, i, j);
             /* gradient */
             if (Tij >= 1)
               {

src/mpi/mpi_ops.c

@@ -380,8 +380,8 @@ static int size_basis(const PnlObject *Obj, MPI_Comm comm, int *size)
   *size += count;
   if (BASIS_HAS_TENSOR_REP(B))
     {
-      /* B->T */
-      if ((info = pnl_object_mpi_pack_size(PNL_OBJECT(B->T), comm, &count))) return info;
+      /* B->SpT */
+      if ((info = pnl_object_mpi_pack_size(PNL_OBJECT(B->SpT), comm, &count))) return info;
       *size += count;
     }
   else
@@ -744,7 +744,7 @@ static int pack_basis(const PnlObject *Obj, void *buf, int bufsize, int *pos, MP
   if ((info = MPI_Pack(&B->id, 1, MPI_INT, buf, bufsize, pos, comm))) return info;
   if (BASIS_HAS_TENSOR_REP(B))
     {
-      if ((info = pnl_object_mpi_pack(PNL_OBJECT(B->T), buf, bufsize, pos, comm))) return info;
+      if ((info = pnl_object_mpi_pack(PNL_OBJECT(B->SpT), buf, bufsize, pos, comm))) return info;
     }
   else
     {
@@ -941,6 +941,7 @@ static int unpack_sp_matrix(PnlObject *Obj, void *buf, int bufsize, int *pos, MP
   if ((info = MPI_Unpack(buf, bufsize, pos, &n, 1, MPI_INT, comm))) return info;
   if ((info = MPI_Unpack(buf, bufsize, pos, &nz, 1, MPI_INT, comm))) return info;
   pnl_sp_mat_object_resize(M, m, n, nz);
+  M->nz = nz;
   if ((info = MPI_Unpack(buf, bufsize, pos, M->I, m + 1, MPI_INT, comm))) return info;
   if ((info = MPI_Unpack(buf, bufsize, pos, M->J, nz, MPI_INT, comm))) return info;
   switch (PNL_GET_TYPE(M))
@@ -1117,13 +1118,12 @@ static int unpack_basis(PnlObject *Obj, void *buf, int bufsize, int *pos, MPI_Co
   pnl_basis_set_type(B, basis_id);
   if (BASIS_HAS_TENSOR_REP(B))
     {
-      PnlMatInt *T = pnl_mat_int_new();
-      if ((info = pnl_object_mpi_unpack(PNL_OBJECT(T), buf, bufsize, pos, comm))) return info;
-      pnl_basis_set_from_tensor(B, T);
+      PnlSpMatInt *SpT = pnl_sp_mat_int_new();
+      if ((info = pnl_object_mpi_unpack(PNL_OBJECT(SpT), buf, bufsize, pos, comm))) return info;
+      pnl_basis_set_from_sparse_tensor(B, SpT);
     }
   else
     {
-      B->T = NULL;
       B->SpT = NULL;
       if ((info = MPI_Unpack(buf, bufsize, pos, &B->nb_func, 1, MPI_INT, comm))) return info;
       if ((info = MPI_Unpack(buf, bufsize, pos, &B->nb_variates, 1, MPI_INT, comm))) return info;