Commit b1c29e65 authored by Antoine Lucas's avatar Antoine Lucas
Browse files

version 0.8.9 - remove memory leaks

parent 157da0cd
Tue Dec 03 2013
* remove [again] printf in C++ code
* fix memory leaks
* add non-regression tests.
Wed Sep 11 2013
* remove printf in C/c++ code
......
Package: amap
Version: 0.8-8
Date: 2013-12-02
Version: 0.8-9
Date: 2013-12-03
Title: Another Multidimensional Analysis Package
Author: Antoine Lucas
Maintainer: Antoine Lucas <antoinelucas@gmail.com>
......
......@@ -4,7 +4,7 @@
* \brief Robust principal component analysis
*
* \date Created : 06/11/02
* \date Last Modified : Time-stamp: <2013-09-11 20:15:45 antoine>
* \date Last Modified : Time-stamp: <2013-12-03 22:05:39 antoine>
*
* This Message must remain attached to this source code at all times.
* This source code may not be redistributed, nor may any executable
......@@ -24,7 +24,8 @@
#ifndef M_PIl
#define M_PIl 3.1415926535897932384626433832795029L /* pi */
#endif
#include "smartPtr.h"
#include "acprob.h"
/* Compiliation: */
/* R CMD SHLIB acprob.c */
......@@ -136,22 +137,12 @@ void mult(double *x,int *p,double *res)
*/
void W(double *x,double *h,double *d,int *n,int *p,char **kernel,double *res, int * result)
{
double *delta;
double *temp;
SmartPtr<double> delta (*p);
SmartPtr<double> temp (*p * *p);
double N=0,K=0,som=0;
int i,j,k,l;
*result = 1;
delta = (double*) malloc (*p * sizeof(double));
temp = (double*) malloc (*p * *p * sizeof(double));
if(delta == NULL || temp == NULL )
{
Rprintf("AMAP: Not enought system memory \n");
*result = 2;
return;
}
for (l=0; l < (*p * *p) ; l++)
res[l]=0;
......@@ -162,14 +153,14 @@ void W(double *x,double *h,double *d,int *n,int *p,char **kernel,double *res, in
/* delta = Xi-Xj (ligne i et j de la matrice X) */
for (k=0; k < *p ; k++)
delta[k]=x[i+(k * *n)]- x[j+(k * *n)];
N = norm(delta,p,d)/ *h;
N = norm(delta.get(),p,d)/ *h;
/* tmp2 = K ( |delta/h|^2 ) */
noyau(&N,kernel,&K);
som += K;
/* temp = delta * delta' (matrice) */
mult(delta,p,temp);
mult(delta.get(),p,temp.get());
for (l=0; l < (*p * *p) ; l++)
res[l] += temp[l] * K ;
}
......@@ -178,8 +169,7 @@ void W(double *x,double *h,double *d,int *n,int *p,char **kernel,double *res, in
for (l=0; l < (*p * *p) ; l++)
res[l] = res[l] / som ;
free(delta);
free(temp);
*result = 0;
}
......@@ -202,18 +192,11 @@ void W(double *x,double *h,double *d,int *n,int *p,char **kernel,double *res, in
void VarRob(double *x,double *h,double *d,int *n,int *p,char **kernel,double *res, int * result)
{
int i,j;
double *temp, *Xi;
SmartPtr<double> temp (*p * *p);
SmartPtr<double> Xi (*p);
double N=0,K=0,som=0;
*result = 1;
temp = (double*) malloc (*p * *p * sizeof(double));
Xi = (double*) malloc (*p * sizeof(double));
if(temp == NULL || Xi == NULL )
{
printf("AMAP: Not enought system memory \n");
*result = 2;
return;
}
som = 0;
......@@ -222,10 +205,10 @@ void VarRob(double *x,double *h,double *d,int *n,int *p,char **kernel,double *re
for (j=0; j < *p ; j++)
Xi[j]=x[i+(j * *n)];
N = norm(Xi,p,d)/ *h;
N = norm(Xi.get(),p,d)/ *h;
noyau(&N,kernel,&K);
mult(Xi,p,temp);
mult(Xi.get(),p,temp.get());
for (j=0; j < (*p * *p) ; j++)
res[j] += temp[j] * K ;
som += K;
......@@ -234,8 +217,7 @@ void VarRob(double *x,double *h,double *d,int *n,int *p,char **kernel,double *re
for (j=0; j < (*p * *p) ; j++)
res[j] = res[j] / som ;
free(Xi);
free(temp);
*result = 0;
}
......
#ifndef ACPROB_H
#define ACPROB_H 1
extern "C" {
/*! noyau: base function for kernel computation
* \brief compute a kernel. called by W or VarRob
* \param u: input data (scalar)
* \param k: char, g, q, t, e, c, u for:
* gaussien = (2*pi)^(-1/2) * exp(-u^2/2),
* quartic = 15/16 * (1-u^2)^2 * (abs(u)<1),
* triweight = 35/32 * (1-u^2)^3 * (abs(u)<1),
* epanechikov = 3/4 * (1-u^2) * (abs(u)<1),
* cosinus = pi/4 * cos (u*pi/2) * (abs(u)<1),
* uniform = 1/2 * (abs(u)<1),
* \param res: output data: one value of type double
*/
void noyau(double *u, char **k,double *res) ;
/** W compute a "local" variance This function is called directly by R
* \brief this functions compute a "local" variance, using a specific kernel.
* This function depends on function mult, norm and noyau
* \param x: matrix; input data n x p
* \param h: window size (scalar)
* \param d: scalar product matrix p x p
* \param n: length x
* \param p: number of columns of x
* \param kernel: kernel utilised
* \param res: matrix returned (allocated in R)
*
* \note
* matrix x[i,j] (n x p) is substitute
* by vector vecteur x[i + j x n] (i=0..n-1 ; j=0..p-1)
*
* \param result flag 0 => correct
* 1 => Pb
* 2 => Cannot allocate memory
*/
void W(double *x,double *h,double *d,int *n,int *p,char **kernel,double *res, int * result);
/** VarRob: compute robust variance. Function called birecly by R
* \brief Robust variance... gives a low ponderation to isolated
* values. This ponderation is determined with kernel and window size.
* This function depends on function mult, norm and noyau
* \param x: data matrix n x p
* \param h: kernel window size (scalar)
* \param d: matrix of scalar product p x p
* \param n: lenght of x
* \param p: length of x
* \param kernel: kernel used
* \param res result (matrix)
* \param result flag 0 => correct
* 1 => Pb
* 2 => Cannot allocate memory
*/
void VarRob(double *x,double *h,double *d,int *n,int *p,char **kernel,double *res, int * result);
};
#endif
......@@ -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: <2012-01-08 12:16:41 antoine>
* \date Last modified: Time-stamp: <2013-12-03 21:38:11 antoine>
*
* \author R core members, and lately: Antoine Lucas
*
......@@ -516,12 +516,12 @@ template<class T> T distance_T<T>::R_spearman(double * x, double * y , int nr_x
int * flag, T_tri & opt)
{
int j;
double * data_tri_x = opt.data_tri_x;
int * order_tri_x = opt.order_tri_x;
int * rank_tri_x = opt.rank_tri_x;
double * data_tri_y = opt.data_tri_y;
int * order_tri_y = opt.order_tri_y;
int * rank_tri_y = opt.rank_tri_y;
double * data_tri_x = opt.data_tri_x.get();
int * order_tri_x = opt.order_tri_x.get();
int * rank_tri_x = opt.rank_tri_x.get();
double * data_tri_y = opt.data_tri_y.get();
int * order_tri_y = opt.order_tri_y.get();
int * rank_tri_y = opt.rank_tri_y.get();
int n;
T diffrang=0;
......@@ -658,12 +658,12 @@ template<class T> T distance_T<T>::R_kendall(double * x, double * y , int nr_x,
int * flag, T_tri & opt)
{
int j,k;
double * data_tri_x = opt.data_tri_x;
int * order_tri_x = opt.order_tri_x;
int * rank_tri_x = opt.rank_tri_x;
double * data_tri_y = opt.data_tri_y;
int * order_tri_y = opt.order_tri_y;
int * rank_tri_y = opt.rank_tri_y;
double * data_tri_x = opt.data_tri_x.get();
int * order_tri_x = opt.order_tri_x.get();
int * rank_tri_x = opt.rank_tri_x.get();
double * data_tri_y = opt.data_tri_y.get();
int * order_tri_y = opt.order_tri_y.get();
int * rank_tri_y = opt.rank_tri_y.get();
int n;
T dist,P=0;
bool ordre_x,ordre_y;
......@@ -763,7 +763,7 @@ template<class T> void distance_T<T>::distance(double *x, int *nr,
arguments[i].method = method;
arguments[i].nbprocess= *nbprocess;
arguments[i].ierr=ierr;
arguments[i].i2=i2;
arguments[i].i2=i2;
}
*ierr = 1; /* res = 1 => no missing values
res = 0 => missings values */
......@@ -889,15 +889,7 @@ template <class T> void* distance_T<T>::thread_dist(void* arguments_void)
if( (*method == SPEARMAN) || (*method == KENDALL))
{
opt.data_tri_x = (double * ) malloc ( (nc) * sizeof(double));
opt.order_tri_x = (int * ) malloc ( (nc) * sizeof(int));
opt.rank_tri_x = (int * ) malloc ( (nc) * sizeof(int));
opt.data_tri_y = (double * ) malloc ( (nc) * sizeof(double));
opt.order_tri_y = (int * ) malloc ( (nc) * sizeof(int));
opt.rank_tri_y = (int * ) malloc ( (nc) * sizeof(int));
if( (opt.data_tri_x == NULL) || (opt.order_tri_x == NULL) || (opt.rank_tri_x == NULL) ||
(opt.data_tri_y == NULL) || (opt.order_tri_y == NULL) || (opt.rank_tri_y == NULL))
error("distance(): unable to alloc memory");
opt.reset(nc);
}
/*
......@@ -943,17 +935,6 @@ template <class T> void* distance_T<T>::thread_dist(void* arguments_void)
}
}
if( (*method == SPEARMAN) || (*method == KENDALL))
{
free(opt.data_tri_x);
free(opt.rank_tri_x);
free(opt.order_tri_x);
free(opt.data_tri_y);
free(opt.rank_tri_y);
free(opt.order_tri_y);
}
return (void*)0;
}
......
......@@ -2,7 +2,7 @@
#ifndef _AMAP_DISTANCE_TEMPLATE
#define _AMAP_DISTANCE_TEMPLATE 1
#include "smartPtr.h"
template<class T> class distance_T
{
......@@ -12,14 +12,32 @@ template<class T> class distance_T
enum { EUCLIDEAN=1, MAXIMUM, MANHATTAN, CANBERRA, BINARY ,PEARSON, CORRELATION, SPEARMAN, KENDALL, ABSPEARSON, ABSCORRELATION};
struct T_tri
class T_tri
{
double * data_tri_x;
int * order_tri_x;
int * rank_tri_x;
double * data_tri_y;
int * order_tri_y;
int * rank_tri_y;
public:
SmartPtr<double> data_tri_x;
SmartPtr<int> order_tri_x;
SmartPtr<int> rank_tri_x;
SmartPtr<double> data_tri_y;
SmartPtr<int> order_tri_y;
SmartPtr<int> rank_tri_y;
T_tri() :
data_tri_x(0),
order_tri_x(0),
rank_tri_x(0),
data_tri_y(0),
order_tri_y(0),
rank_tri_y(0) {};
void reset(int size) {
data_tri_x.reset(size);
order_tri_x.reset(size);
rank_tri_x.reset(size);
data_tri_y.reset(size);
order_tri_y.reset(size);
rank_tri_y.reset(size);
}
};
......
......@@ -7,8 +7,10 @@
#include "distance.h"
#include "hclust_T.h"
#include <stdio.h>
#include "smartPtr.h"
#include "R.h"
#include <new>
namespace hclust_T
{
......@@ -34,40 +36,34 @@ namespace hclust_T
*/
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 len = (*nr * (*nr-1) )/2;
int flag;
T *d;
SmartPtr<T> d (len);
*result = 1;
len = (*nr * (*nr-1) )/2;
d = (T*) malloc (len * sizeof(T));
if(d == NULL )
{
Rprintf("AMAP: Not enought system memory to allocate matrix distance\n");
*result = 2;
return;
}
/*
* Calculate d: distance matrix
*/
distance_T<T>::distance(x,nr,nc,d,diag,method,nbprocess,&flag,-1);
distance_T<T>::distance(x,nr,nc,d.get(),diag,method,nbprocess,&flag,-1);
if(flag == 0)
{
Rprintf("AMAP: Unable to compute Hierarchical Clustering: missing values in distance matrix\n");
*result = 3;
return;
}
/*
* Hierarchical clustering
*/
hclust_T::hclust<T>(nbprocess,x,*nr,*nc,method,nr,&len,iopt ,ia ,ib,iorder,crit,membr,d,result);
free(d);
hclust_T::hclust<T>(nbprocess,x,*nr,*nc,method,nr,&len,iopt ,ia ,ib,iorder,crit,membr,d.get(),result);
*result = 0;
}
......@@ -78,37 +74,28 @@ namespace hclust_T
{
int im,jm,jj,i,j,ncl,ind,i2,j2,k,ind1,ind2,ind3;
double inf,dmin,x,xx;
int *nn;
int *items = NULL;
double *disnn;
short int *flag;
int *iia;
int *iib;
SmartPtr<int> nn (*n);
SmartPtr<int> items (0);
SmartPtr<double> disnn(*n);
SmartPtr<short int> flag(*n);
int count,h,idx1,idxitem1,idx2;
*result = 1;
nn = (int*) malloc (*n * sizeof(int));
if(*iopt==CENTROID2)
{
items = (int*) malloc (*n*nc * sizeof(int));
items.reset(*n * nc);
}
disnn = (double*) malloc (*n * sizeof(double));
flag = (short int*) malloc (*n * sizeof(short int));
if(nn == NULL || disnn == NULL || flag == NULL )
{
Rprintf("AMAP: Not enought system memory \n");
*result = 2;
return;
}
/* Initialisation */
for ( i=0; i<*n ; i++)
{
flag[i]=1;
idxitem1=i;
if(items != NULL)
if(items.get() != NULL)
{
for ( j=0; j<nc ; j++)
{
......@@ -375,25 +362,11 @@ namespace hclust_T
} /* end of while */
free(nn);
free(disnn);
free(flag);
if(items != NULL )
free(items);
iia = (int*) malloc (*n * sizeof(int));
iib = (int*) malloc (*n * sizeof(int));
if(iia == NULL || iib == NULL )
{
Rprintf("AMAP: Not enought system memory \n");
*result = 2;
return;
}
SmartPtr<int> iia (*n );
SmartPtr<int> iib (*n );
hierclust::hcass2(n,ia,ib,iorder,iia,iib);
hierclust::hcass2(n,ia,ib,iorder,iia.get(),iib.get());
/*
......@@ -406,8 +379,6 @@ namespace hclust_T
ib[i]= iib[i];
}
free(iia);
free(iib);
*result = 0;
}
......
......@@ -4,7 +4,7 @@
* \brief K-means clustering
*
* \date Created : before 2005
* \date Last Modified : Time-stamp: <2011-11-03 22:19:14 antoine>
* \date Last Modified : Time-stamp: <2013-12-03 21:42:43 antoine>
*
* \author R core team. Modified by A. Lucas for distance choice.
*
......@@ -68,16 +68,7 @@ void kmeans_Lloyd2(double *x, int *pn, int *pp, double *cen, int *pk, int *cl,
if( (*method == distance_T<double>::SPEARMAN) || (*method == distance_T<double>::KENDALL))
{
opt.data_tri_x = (double * ) malloc ( p * sizeof(double));
opt.order_tri_x = (int * ) malloc ( p * sizeof(int));
opt.rank_tri_x = (int * ) malloc ( p * sizeof(int));
opt.data_tri_y = (double * ) malloc ( p * sizeof(double));
opt.order_tri_y = (int * ) malloc ( p * sizeof(int));
opt.rank_tri_y = (int * ) malloc ( p * sizeof(int));
if( (opt.data_tri_x == NULL) || (opt.order_tri_x == NULL) || (opt.rank_tri_x == NULL) ||
(opt.data_tri_y == NULL) || (opt.order_tri_y == NULL) || (opt.rank_tri_y == NULL))
error("distance(): unable to alloc memory");
opt.reset(p);
}
for(i = 0; i < n; i++) cl[i] = -1;
......@@ -126,13 +117,4 @@ void kmeans_Lloyd2(double *x, int *pn, int *pp, double *cen, int *pk, int *cl,
}
if( (*method == distance_T<double>::SPEARMAN) || (*method == distance_T<double>::KENDALL))
{
free(opt.data_tri_x);
free(opt.rank_tri_x);
free(opt.order_tri_x);
free(opt.data_tri_y);
free(opt.rank_tri_y);
free(opt.order_tri_y);
}
}
library(amap)
set.seed(1234)
data(USArrests)
METHODS <- c("euclidean", "maximum", "manhattan", "canberra",
"binary","pearson","correlation","spearman","kendall",
"abspearson","abscorrelation")
METHODSLINKS <- c("ward", "single", "complete", "average", "mcquitty",
"median", "centroid","centroid2")
for (mymethod in METHODS) {
d = Dist(USArrests, method = mymethod)
print(d)
k = Kmeans(USArrests, centers = 4, method = mymethod)
print(k)
for (mylink in METHODSLINKS)
{
cat(mylink)
cat(mymethod)
hc <- hcluster(USArrests,link = mylink, method = mymethod, nbproc=4)
print(hc)
}
}
hc <- hcluster(USArrests, nbproc=1)
print(hc)
KERNELS = c("gaussien", "quartic", "triweight", "epanechikov" ,
"cosinus", "uniform")
for(myKernel in KERNELS) {
myacp = acprob(USArrests, kernel = myKernel);
print(myacp)
}
This diff is collapsed.
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