Commit 34dd89a0 authored by Antoine Lucas's avatar Antoine Lucas
Browse files

suppression of duplicated code(used for no thread / multiple thread)

parent d77f8fea
#include "distance.h"
#include "distance_template.h"
#include "distance_T.h"
......@@ -19,15 +19,9 @@
*/
void R_distance(double *x, int *nr, int *nc, double *d, int *diag, int *method,int *nbprocess, int * ierr)
{
#ifndef __MINGW_H
if(*nbprocess != 1)
amaptemplate::distancepar<double>(x,nr,nc, d,diag,method,
nbprocess, ierr);
else
#endif
amaptemplate::T_distance<double>(x,nr, nc, d, diag,
method, ierr);
distance_T::distance<double>(x,nr,nc, d,diag,method,
nbprocess, ierr);
}
......
......@@ -35,7 +35,7 @@
#include <config.h>
#endif
#include "distance_template.h"
#include "distance_T.h"
#include "distance.h"
#include <float.h>
......@@ -48,8 +48,9 @@
#define MAX( A , B ) ( ( A ) > ( B ) ? ( A ) : ( B ) )
#define MIN( A , B ) ( ( A ) < ( B ) ? ( A ) : ( B ) )
#define __MINGW_H 1
namespace amaptemplate
namespace distance_T
{
/** \brief Distance euclidean (i.e. sqrt(sum of square) )
......@@ -311,7 +312,7 @@ namespace amaptemplate
return( diffrang );
/*
* verifcation in R:
* verification in R:
* Dist(x,method='spearman') ; n =dim(x)[2]
* l=c(x[3,],x[4,]); sum((rank(l[1:n])-rank(l[(n+1):(2*n)]))^2)
* cor.test(x[3,],x[4,],method="spearm")$statistic
......@@ -334,13 +335,12 @@ namespace amaptemplate
* \param nbprocess: number of threads to create
* \param ierr error return; 1 good; 0 missing values
*/
template <class T > void distancepar(double *x, int *nr, int *nc, T *d, int *diag, int *method,int *nbprocess, int * ierr)
template <class T > void distance(double *x, int *nr, int *nc, T *d, int *diag, int *method,int *nbprocess, int * ierr)
{
#ifndef __MINGW_H
int dc;
short int * jobs;
int i;
pthread_t * th;
void ** arguments;
dc = (*diag) ? 0 : 1; /* diag=1: we do the diagonal */
......@@ -357,30 +357,37 @@ namespace amaptemplate
* *ierr
*/
/* numero_thread=0;*/
arguments = (void **) malloc ( 9 * sizeof( void *));
jobs = (short int *) malloc ( *nr * sizeof(short int));
th = (pthread_t *) malloc ( *nbprocess * sizeof(pthread_t));
for(i=0; i< *nr; i++){jobs[i]=0;}
/* arguments[0] = &numero_thread ;*/
arguments[0] = jobs;
arguments[1] = nr ;
arguments[2] = nc ;
arguments[3] = &dc ;
arguments[4] = x ;
arguments[5] = d ;
arguments[6] = method ;
arguments[7] = nbprocess;
arguments[8] = ierr;
arguments = (void **) malloc (9 * (*nbprocess) * sizeof( void *));
jobs = (short int *) malloc ( *nbprocess * sizeof(short int));
for(i=0; i< *nbprocess; i++){jobs[i]=i;}
//printf("nb processs %d\n",*nbprocess);
for(i=0; i< *nbprocess; ++i)
{
arguments[0 + 9*i] = jobs+i;
arguments[1 + 9*i] = nr ;
arguments[2 + 9*i] = nc ;
arguments[3 + 9*i] = &dc ;
arguments[4 + 9*i] = x ;
arguments[5 + 9*i] = d ;
arguments[6 + 9*i] = method ;
arguments[7 + 9*i] = nbprocess;
arguments[8 + 9*i] = ierr;
}
*ierr = 1; /* res = 1 => no missing values
res = 0 => missings values */
#ifndef __MINGW_H
pthread_t * th = (pthread_t *) malloc ( *nbprocess * sizeof(pthread_t));
for (i=0; i < *nbprocess ; i++)
{
pthread_create(th+i,0,thread_dist<T>,(void *)arguments);
pthread_create(th+i,0,thread_dist<T>,(void *)(arguments+(i*9)));
}
/* Attends la fin */
......@@ -388,27 +395,32 @@ namespace amaptemplate
{
pthread_join(*(th+i),NULL);
}
free( arguments );
free( jobs);
free( th);
#else
T_distance<T>(x, nr, nc, d,diag,method,ierr);
// p_thread not yet used on windows.
int nombre=1;
arguments[7 ] = &nombre;
thread_dist<T>((void *)arguments);
#endif
free( arguments );
free( jobs);
}
#ifndef __MINGW_H
/** thread_dist function that compute distance.
*
*/
template <class T> void* thread_dist(void* arguments)
{
long int nbprocess,nr,nc,i,j,dc,ij;
short int *jobs;
int nbprocess,nr,nc,i,j,dc,ij;
void ** arguments2;
T * d;
double * x;
......@@ -425,7 +437,7 @@ namespace amaptemplate
arguments2 = (void **) arguments;
jobs = (short int *) arguments2[0];
short int no = * (short int *) arguments2[0];
nr = * (int *) arguments2[1];
nc = * (int *) arguments2[2];
dc = * (int *) arguments2[3];
......@@ -435,13 +447,7 @@ namespace amaptemplate
nbprocess = * (int *) arguments2[7];
ierr = (int *) arguments2[8];
/*
no = * (int *) arguments2[0];*/
/* Increment du thread */
/* tmp = (int *) arguments2[0];
*tmp = no+1;
*/
switch(*method) {
case EUCLIDEAN:
distfun = R_euclidean<T>;
......@@ -490,18 +496,17 @@ namespace amaptemplate
/* debut des boucles 0
fin: nr+1 */
/*
debut = (long int) floor( ((nr+1.)*nbprocess - sqrt( (nr+1.)*(nr+1.) * nbprocess * nbprocess - (nr+1.)*(nr+1.) * nbprocess * no ) )/nbprocess);
fin = (long int) floor(((nr+1.)*nbprocess - sqrt( (nr+1.)*(nr+1.) * nbprocess * nbprocess - (nr+1.)*(nr+1.) * nbprocess * (no+1.) ) )/nbprocess);
printf("Thread %d debut %d fin %d\n",no,debut,fin);
int debut = (int) floor( ((nr+1.)*nbprocess - sqrt( (nr+1.)*(nr+1.) * nbprocess * nbprocess - (nr+1.)*(nr+1.) * nbprocess * no ) )/nbprocess);
int fin = (int) floor(((nr+1.)*nbprocess - sqrt( (nr+1.)*(nr+1.) * nbprocess * nbprocess - (nr+1.)*(nr+1.) * nbprocess * (no+1.) ) )/nbprocess);
printf("Thread %d debut %d fin %d\n",no,debut,fin);
*/
for(j = 0 ; j <= nr ; j++)
if(jobs[j]==0)
// for(j = 0 ; j <= nr ; j++)
for(j = debut ; j < fin ; j++)
{
jobs[j]=1;
ij = (2 * (nr-dc) - j +1) * (j) /2 ;
for(i = j+dc ; i < nr ; i++)
{
......@@ -522,82 +527,10 @@ namespace amaptemplate
}
template <class T> void T_distance(double *x, int *nr, int *nc,
T *d, int *diag,
int *method, int * ierr)
{
int dc, i, j, ij;
/* for spearman dist */
void ** opt;
double * data_tri;
int * order_tri;
int * rank_tri;
T (*distfun)(double*, int, int, int, int, int *, void **) = NULL;
switch(*method) {
case EUCLIDEAN:
distfun = R_euclidean<T>;
break;
case MAXIMUM:
distfun = R_maximum<T>;
break;
case MANHATTAN:
distfun = R_manhattan<T>;
break;
case CANBERRA:
distfun = R_canberra<T>;
break;
case BINARY:
distfun = R_dist_binary<T>;
break;
case PEARSON:
distfun = R_pearson<T>;
break;
case CORRELATION:
distfun = R_correlation<T>;
break;
case SPEARMAN:
distfun = R_spearman<T>;
opt = (void ** ) malloc ( 3 * sizeof(void*));
data_tri = (double * ) malloc (2* (*nc) * sizeof(double));
order_tri = (int * ) malloc (2 * (*nc) * sizeof(int));
rank_tri = (int * ) malloc (2 * (*nc) * sizeof(int));
if( (data_tri == NULL) || (order_tri == NULL) || (rank_tri == NULL))
error("distance(): unable to alloc memory");
opt[0] = (void *) data_tri;
opt[1] = (void *) order_tri;
opt[2] = (void *) rank_tri;
break;
default:
error("distance(): invalid distance");
}
*ierr = 1; /* res = 1 => no missing values
res = 0 => missings values */
dc = (*diag) ? 0 : 1; /* diag=1: we do the diagonal */
ij = 0;
for(j = 0 ; j <= *nr ; j++)
for(i = j+dc ; i < *nr ; i++)
d[ij++] = distfun(x, *nr, *nc, i, j,ierr,opt);
if((*method) == SPEARMAN)
{
free(opt);
free(data_tri);
free(order_tri);
free(rank_tri);
}
/* If spearman: free(vecteur pour tri) */
}
}
}
#endif
......@@ -5,9 +5,16 @@
/* == 1,2,..., defined by order in the r function dist */
enum { EUCLIDEAN=1, MAXIMUM, MANHATTAN, CANBERRA, BINARY ,PEARSON, CORRELATION, SPEARMAN};
namespace amaptemplate
namespace distance_T
{
template <class T> struct arguments
{
double i;
T j;
};
/** \brief Distance euclidean (i.e. sqrt(sum of square) )
*/
template <class T> T R_euclidean(double * x, int nr, int nc, int i1, int i2,int * flag, void ** opt);
......@@ -44,23 +51,17 @@ namespace amaptemplate
template <class T> T R_spearman(double * x, int nr, int nc, int i1, int i2,int * flag, void ** opt);
template <class T > void distancepar(double *x, int *nr, int *nc, T *d, int *diag, int *method,int *nbprocess, int * ierr);
template <class T > void distance(double *x, int *nr, int *nc, T *d, int *diag, int *method,int *nbprocess, int * ierr);
template <class T> void* thread_dist(void* arguments);
template <class T> void T_distance(double *x, int *nr, int *nc,
T *d, int *diag,
int *method, int * ierr);
}
#endif
#ifndef _AMAP_DISTANCE_TEMPLATE_CPP
#define _AMAP_DISTANCE_TEMPLATE_CPP 1
#include "distance_template.cpp"
#include "distance_T.cpp"
#endif
......@@ -47,7 +47,7 @@
#include "mva.h"
#include "hclust.h"
#include "distance.h"
#include "hclust_template.h"
#include "hclust_T.h"
#define max( A , B ) ( ( A ) > ( B ) ? ( A ) : ( B ) )
#define min( A , B ) ( ( A ) < ( B ) ? ( A ) : ( B ) )
......@@ -73,8 +73,8 @@
*/
void hclust(int *n,int *len, int *iopt ,int *ia , int *ib,int *iorder,double *crit,double *membr,double *diss,int *result)
{
amaptemplate::hclust<double>(n,len,iopt ,ia ,ib,iorder,
crit,membr,diss, result);
hclust_T::hclust<double>(n,len,iopt ,ia ,ib,iorder,
crit,membr,diss, result);
} /* end function hclust */
......@@ -105,29 +105,18 @@ void hclust(int *n,int *len, int *iopt ,int *ia , int *ib,int *iorder,double *cr
void hcluster(double *x, int *nr, int *nc, int *diag, int *method, int *iopt ,int *ia , int *ib,int *iorder,double *crit,double *membr,int *nbprocess,int * precision, int * result)
{
#ifndef __MINGW_H
if(*nbprocess > 1)
{
if(*precision == 1)
amaptemplate::hclusterpar<float>(x,nr,nc,diag,method, iopt ,
ia ,ib,iorder,crit,membr,
nbprocess, result);
else
amaptemplate::hclusterpar<double>(x,nr,nc,diag,method, iopt ,
ia ,ib,iorder,crit,membr,
nbprocess, result);
return;
}
#endif
if(*precision == 1)
amaptemplate::hcluster<float>(x,nr,nc,diag,method, iopt ,
hclust_T::hcluster<float>(x,nr,nc,diag,method, iopt ,
ia ,ib,iorder,crit,membr,
result);
nbprocess, result);
else
amaptemplate::hcluster<double>(x,nr,nc,diag,method, iopt ,
hclust_T::hcluster<double>(x,nr,nc,diag,method, iopt ,
ia ,ib,iorder,crit,membr,
result);
nbprocess, result);
return;
}
......
......@@ -2,15 +2,34 @@
#define _AMAP_HCLUST_TEMPLATE_CPP 1
#include "hclust_template.h"
#include "distance_template.h"
#include "hclust_T.h"
#include "distance_T.h"
#include "hclust.h"
namespace amaptemplate
namespace hclust_T
{
template <class T> void hclusterpar(double *x, int *nr, int *nc, int *diag, int *method, int *iopt ,int *ia , int *ib,int *iorder,double *crit,double *membr,int *nbprocess, int * result)
/** hierarchical clustering
* \brief allocate distance matrix execute function R_distance, launch
* hclust on this distance matrix, and free memory of distance matrix.
* \param x: data nr x nc
* \param nr,nc number of row and columns
* \param membr: member, vector 1:nr
* \param method integer -> distance method
* \param diag integer: if we compute diagonal in distance matrix (usually no)
* \param iopt integer -> link used
* \param ia, ib result (merge)
* \param iorder result (order)
* \param crit result (height)
* \param nbprocess number of processes (thread) used
* \param result flag 0 => correct
* 1 => Pb
* 2 => Cannot allocate memory
* 3 => Pb with distance matrix
*
*/
template <class T> void hcluster(double *x, int *nr, int *nc, int *diag, int *method, int *iopt ,int *ia , int *ib,int *iorder,double *crit,double *membr,int *nbprocess, int * result)
{
int len;
int flag;
......@@ -18,7 +37,6 @@ namespace amaptemplate
*result = 1;
#ifndef __MINGW_H
len = (*nr * (*nr-1) )/2;
d = (T*) malloc (len * sizeof(T));
......@@ -31,7 +49,7 @@ namespace amaptemplate
/*
* Calculate d: distance matrix
*/
amaptemplate::distancepar<T>(x,nr,nc,d,diag,method,nbprocess,&flag);
distance_T::distance<T>(x,nr,nc,d,diag,method,nbprocess,&flag);
if(flag == 0)
{
printf("AMAP: Unable to compute Hierarchical Clustering: missing values in distance matrix\n");
......@@ -42,18 +60,11 @@ namespace amaptemplate
/*
* Hierarchical clustering
*/
amaptemplate::hclust<T>(nr,&len,iopt ,ia ,ib,iorder,crit,membr,d,result);
hclust_T::hclust<T>(nr,&len,iopt ,ia ,ib,iorder,crit,membr,d,result);
free(d);
*result = 0;
#else
/* No threads on Windows yet */
amaptemplate::hcluster<T>(x, nr, nc, diag, method, iopt ,ia ,ib,iorder, crit, membr, result);
#endif
}
......@@ -304,64 +315,4 @@ namespace amaptemplate
}
/** hierarchical clustering
* \brief allocate distance matrix execute function R_distance, launch
* hclust on this distance matrix, and free memory of distance matrix.
* \param x: data nr x nc
* \param nr,nc number of row and columns
* \param membr: member, vector 1:nr
* \param method integer -> distance method
* \param diag integer: if we compute diagonal in distance matrix (usually no)
* \param iopt integer -> link used
* \param ia, ib result (merge)
* \param iorder result (order)
* \param crit result (height)
* \param result flag 0 => correct
* 1 => Pb
* 2 => Cannot allocate memory
* 3 => Pb with distance matrix
*/
template <class T> void hcluster(double *x, int *nr, int *nc,
int *diag, int *method, int *iopt ,
int *ia , int *ib,int *iorder,
double *crit,double *membr,int *result)
{
int len;
int flag;
T *d;
*result = 1;
len = (*nr * (*nr-1) )/2;
d = (T*) malloc (len * sizeof(double));
if(d == NULL )
{
printf("AMAP: Not enought system memory to allocate matrix distance\n");
*result = 2;
return;
}
/*
* Calculate d: distance matrix
*/
amaptemplate::T_distance<T>(x,nr,nc,d,diag,method,&flag);
if(flag == 0)
{
printf("AMAP: Unable to compute Hierarchical Clustering: missing values in distance matrix\n");
*result = 3;
return;
}
/*
* Hierarchical clustering
*/
amaptemplate::hclust<T> (nr,&len,iopt ,ia ,ib,iorder,crit,membr,d,result);
free(d);
*result = 0;
}
}
......@@ -3,17 +3,19 @@
#ifndef _AMAP_HCLUST_TEMPLATE
#define _AMAP_HCLUST_TEMPLATE 1
namespace amaptemplate
namespace hclust_T
{
template <class T> void hclusterpar(double *x, int *nr, int *nc, int *diag, int *method, int *iopt ,int *ia , int *ib,int *iorder,double *crit,double *membr,int *nbprocess, int * result);
template <class T> void hclust(int *n,int *len, int *iopt ,int *ia , int *ib,int *iorder,double *crit,double *membr,T *diss,int *result);
template <class T> void hcluster(double *x, int *nr, int *nc,
int *diag, int *method, int *iopt ,
int *ia , int *ib,int *iorder,
double *crit,double *membr,int *result);
double *crit,double *membr,
int *nbprocess, int * result);
template <class T> void hclust(int *n,int *len, int *iopt ,int *ia ,
int *ib,int *iorder,double *crit,
double *membr,T *diss,int *result);
}
......@@ -22,5 +24,5 @@ namespace amaptemplate
#ifndef _AMAP_HCLUST_TEMPLATE_CPP
#define _AMAP_HCLUST_TEMPLATE_CPP 1
#include "hclust_template.cpp"
#include "hclust_T.cpp"
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment