ajout d'un template pour les matrices triangulaires

src/distance_T.h

@@ -382,7 +382,7 @@ namespace amap {
 
 #ifndef _AMAP_DISTANCE_TEMPLATE_CPP 
 #define _AMAP_DISTANCE_TEMPLATE_CPP 1
-#include "distance_T.cpp"
+#include "distance_T.inl"
 #endif
 
 

src/distance_T.cpp → src/distance_T.inl

@@ -2,7 +2,7 @@
  * \brief all functions requiered for R dist function and C hcluster function.
  *
  *  \date Created: probably in 1995
- *  \date Last modified: Time-stamp: <2014-03-01 13:07:20 antoine>
+ *  \date Last modified: Time-stamp: <2014-03-01 17:21:54 antoine>
  *
  *  \author R core members, and lately: Antoine Lucas 
  *
@@ -36,6 +36,7 @@
 #endif
 
 #include "distance_T.h"
+#include "matriceTriangle.h"
 #include "distance.h"
 
 #include <float.h>
@@ -862,7 +863,7 @@ template <class T> void* distance_T<T>::thread_dist(void* arguments_void)
 
 
   matrice<double> myMatrice (x, nr, nc);
-  
+  matriceTriangle<T> distMatrice(d, nr, false);
   getDistfunction(*method,distfun);
 
   
@@ -896,27 +897,29 @@ template <class T> void* distance_T<T>::thread_dist(void* arguments_void)
     {
       for(j = debut ; j < fin ; j++)
 	{
+	  vecteur<T> distRow = distMatrice.getRow(j);
+	  vecteur<double> rowJ = myMatrice.getRow(j);
+
 	  ij = (2 * (nr-dc) - j +1) * (j) /2 ;
 	  for(i = j+dc ; i < nr ; i++)
 	    {
 	      vecteur<double> rowI = myMatrice.getRow(i);
-	      vecteur<double> rowJ = myMatrice.getRow(j);
-	      d[ij++] = distfun(rowI, rowJ ,ierr,opt);
+	      distRow[i] = distfun(rowI, rowJ ,ierr,opt);
 	    }
 	}
   } 
   else { /* updates the distance only for i2*/
-    //printf("DEBUG AL : %d %d %d\n", i2, nr, nc);
+
+    // a row of distance matrix
+    vecteur<T> distRow = distMatrice.getRow(i2);
+    vecteur<double> rowI = myMatrice.getRow(i2);
+
     for(j = debut ; j < fin ; j++)
       {
 	if (i2!=j) {
-	  int ind1 = j+i2*nr-(i2+1)*(i2+2)/2;
-	  if (i2 > j)
-	    ind1 = i2+j*nr-(j+1)*(j+2)/2;
 
-	  vecteur<double> rowI = myMatrice.getRow(i2);
 	  vecteur<double> rowJ = myMatrice.getRow(j);
-	  d[ind1] = distfun(rowI, rowJ,ierr,opt);
+	  distRow[j] = distfun(rowI, rowJ,ierr,opt);
 	  //printf("updated dist %d %d %f\n",i2,j,d[ind1]);
 	}
 	

src/hclust_T.h

@@ -42,5 +42,5 @@ namespace hclust_T
 
 #ifndef _AMAP_HCLUST_TEMPLATE_CPP
 #define _AMAP_HCLUST_TEMPLATE_CPP 1
-#include "hclust_T.cpp"
+#include "hclust_T.inl"
 #endif

src/matrice.cpp

+#include "matrice.h"
+#include "matriceTriangle.h"
+#include <R.h>
+
+extern "C" {
+  /**
+   * purpose of this function is to display a matrix, for unitary testing.
+   */
+  void checkMatrix(double * values, int * nrow, int * ncol)  {
+    amap::matrice<double> myMatrix(values, *nrow, *ncol);
+
+    Rprintf("\n");
+    for (int i = 0 ; i < *nrow; i++) {
+      amap::vecteur<double> row = myMatrix.getRow(i);
+
+      for (int j  = 0 ; j < row.size(); j++) {
+	Rprintf("%f\t", row[j]);
+      }
+
+      Rprintf("\n");      
+    }
+
+  }
+
+
+  void checkMatrixTriangle(double * values, int * nrow, int * isDiagonal)  {
+    amap::matriceTriangle<double> myMatrix(values, *nrow, *isDiagonal);
+
+    Rprintf("\n");
+    for (int i = 0 ; i < *nrow; i++) {
+      amap::vecteur<double> row = myMatrix.getRow(i);
+
+      for (int j  = 0 ; j < row.size(); j++) {
+	Rprintf("%f\t", row[j]);
+      }
+
+      Rprintf("\n");      
+    }
+
+  }
+
+};

src/matrice.h

@@ -57,7 +57,20 @@ namespace amap {
   };
 
 
-
+  /**
+   * Matrix data.
+   *
+   * a matrix with 3 row and 4 cols
+   * +---+---+---+---+
+   * | 0 | 3 | 6 | 9 |
+   * +---+---+---+---+
+   * | 1 | 4 | 7 | 10|
+   * +---+---+---+---+
+   * | 2 | 5 | 8 | 11|
+   * +---+---+---+---+
+   * 
+   * 
+   */
   template<class T> class matrice : public array<T> {
   private:
     /**
@@ -131,6 +144,22 @@ namespace amap {
       return nrow * ncol;
     };
     
+    /**
+     * accessor.
+     * \return number of rows.
+     */
+    int getNrow() {
+      return nrow;
+    }
+
+    /**
+     * accessor.
+     * \return number of rows.
+     */
+    int getNcol() {
+      return ncol;
+    }
+
 
   };
 

src/matriceTriangle.h

+#ifndef AMAP_TRIANGULAR_MATRIX
+#define AMAP_TRIANGULAR_MATRIX
+
+#include "matrice.h"
+
+namespace amap {
+  /**
+   * Matrix data.
+   *
+   * a triangular matrix with 4 row and 4 cols
+   * +---+---+---+---+
+   * |   |   |   |   |
+   * +---+---+---+---+
+   * | 1 |   |   |   |
+   * +---+---+---+---+
+   * | 2 | 4 |   |   |
+   * +---+---+---+---+
+   * | 3 | 5 | 6 |   |
+   * +---+---+---+---+
+   * 
+   * matrix is  4x4 but dataSet is of a size 6 (when not diagonal).
+   * dataSet is of size 10 if matrix is diagonal.
+   * 
+   */
+
+  template<class T> class matriceTriangle : public matrice<T> {
+
+  private:
+    T blankValue;
+    bool isDiagonal;
+  public:
+    /**
+     * Contructor.
+     * \param values_p the data matrix
+     * \param nrows_p the number of rows
+     */
+    matriceTriangle(T * values_p, int nrow_p, bool isDiagonal_p) 
+      : matrice<T>(values_p, nrow_p, nrow_p),
+	isDiagonal(isDiagonal_p)
+      {
+      };
+  
+    
+    /**
+     * Accessor on data.
+     */
+    virtual T & operator[] (int index) {
+      blankValue = 0;
+
+
+      int i = index % this->getNrow();
+      int j = index / this->getNrow();
+
+      if (i == j && !isDiagonal) {
+	return blankValue;
+      }
+
+      if (i <= j)
+	{
+	  int temp = i;
+	  i = j;
+	  j = temp;
+	}
+      int maxRowSize = this->getNrow();
+      if (!isDiagonal) {
+	i--;
+	maxRowSize--;
+      }
+
+      int ij = i + j * maxRowSize - j * (j + 1) / 2;
+
+
+      return matrice<T>::operator[](ij);
+    };
+
+
+  };
+
+
+};
+
+
+#endif

tests/amap-test.R

@@ -43,3 +43,17 @@ for(myKernel in KERNELS) {
   print(myacp)
 } 
 
+x = matrix (1:12, 3, 4)
+ncol = as.integer(dim(x)[2])
+nrow = as.integer(dim(x)[1])
+rien = .C("checkMatrix",as.double(x), nrow, ncol, PACKAGE="amap")
+
+m = matrix(1:16,4,4)
+d = as.dist(m)
+rien = .C("checkMatrixTriangle",as.double(d),as.integer(4), as.integer(0))
+
+d=c(1,2,3,4,6,7,8,11,12,16)
+rien = .C("checkMatrixTriangle",as.double(d),as.integer(4), as.integer(1))
+
+
+

tests/amap-test.Rout.save

@@ -6425,7 +6425,35 @@ Standard deviations:
 Eigen values:
 [1] 3.2817172 1.0106413 0.4610868 0.3944977
 > 
+> x = matrix (1:12, 3, 4)
+> ncol = as.integer(dim(x)[2])
+> nrow = as.integer(dim(x)[1])
+> rien = .C("checkMatrix",as.double(x), nrow, ncol, PACKAGE="amap")
+
+1.000000	4.000000	7.000000	10.000000	
+2.000000	5.000000	8.000000	11.000000	
+3.000000	6.000000	9.000000	12.000000	
+> 
+> m = matrix(1:16,4,4)
+> d = as.dist(m)
+> rien = .C("checkMatrixTriangle",as.double(d),as.integer(4), as.integer(0))
+
+0.000000	2.000000	3.000000	4.000000	
+2.000000	0.000000	7.000000	8.000000	
+3.000000	7.000000	0.000000	12.000000	
+4.000000	8.000000	12.000000	0.000000	
+> 
+> d=c(1,2,3,4,6,7,8,11,12,16)
+> rien = .C("checkMatrixTriangle",as.double(d),as.integer(4), as.integer(1))
+
+1.000000	2.000000	3.000000	4.000000	
+2.000000	6.000000	7.000000	8.000000	
+3.000000	7.000000	11.000000	12.000000	
+4.000000	8.000000	12.000000	16.000000	
+> 
+> 
+> 
 > 
 > proc.time()
    user  system elapsed 
-  2.250   0.700   2.876 
+  2.160   0.890   3.015