cluster.xs

聚类分析的源码集

	}
	return ( newRV_noinc( (SV*) matrix_av ) );
}

/* -------------------------------------------------
 *
 */
static SV *
row_c2perl_dbl(pTHX_ double * row, int ncols) {

	int j;
	AV * row_av = newAV();
	for(j=0; j<ncols; ++j) {
		av_push(row_av, newSVnv(row[j]));
		/* printf("%d: %7.3f\n", j, row[j]); */
	}
	return ( newRV_noinc( (SV*) row_av ) );
}

/* -------------------------------------------------
 *
 */
static SV*
row_c2perl_int(pTHX_ int * row, int ncols) {

	int j;
	AV * row_av = newAV();
	for(j=0; j<ncols; ++j) {
		av_push(row_av, newSVnv(row[j]));
	}
	return ( newRV_noinc( (SV*) row_av ) );
}

/* -------------------------------------------------
 *
 */
static SV*
matrix_c2perl_int(pTHX_ int ** matrix, int nrows, int ncols) {

	int i;
	AV * matrix_av = newAV();
	SV * row_ref;
	for(i=0; i<nrows; ++i) {
		row_ref = row_c2perl_int(aTHX_ matrix[i], ncols);
		av_push(matrix_av, row_ref);
	}
	return ( newRV_noinc( (SV*) matrix_av ) );
}

/* -------------------------------------------------
 *
 */
static SV*
matrix_c2perl_dbl(pTHX_ double ** matrix, int nrows, int ncols) {

	int i;
	AV * matrix_av = newAV();
	SV * row_ref;
	for(i=0; i<nrows; ++i) {
		row_ref = row_c2perl_dbl(aTHX_ matrix[i], ncols);
		av_push(matrix_av, row_ref);
	}
	return ( newRV_noinc( (SV*) matrix_av ) );
}

/* -------------------------------------------------
 * Convert the 'data' and 'mask' matrices and the 'weight' array
 * from C to Perl.  Also check for errors, and ignore the
 * mask or the weight array if there are any errors. 
 * Print warnings so the user will know what happened. 
 */
static int
malloc_matrices(pTHX_
	SV *  weight_ref, double  ** weight, int * nweights_ptr, 
	SV *  data_ref,   double *** matrix,
	SV *  mask_ref,   int    *** mask,
	int   nrows,      int        ncols
) {

	int error_count;
	int dummy;
	const int nweights = *nweights_ptr; /* The correct number of weights */

	if(SvTYPE(SvRV(mask_ref)) == SVt_PVAV) { 
		error_count = malloc_matrix_perl2c_int(aTHX_ mask_ref, mask, &dummy, &dummy);
		if(error_count > 0) {
			free_matrix_int(*mask, nrows);
			*mask = malloc_matrix_int(aTHX_ nrows,ncols,1);
		}
	} else {
			*mask = malloc_matrix_int(aTHX_ nrows,ncols,1);
	}

	/* We don't check data_ref because we expect the caller to check it 
	 */
	error_count = malloc_matrix_perl2c_dbl(aTHX_ data_ref, matrix, &dummy, &dummy, *mask);
	if (error_count > 0 && warnings_enabled(aTHX)) 
		Perl_warn(aTHX_ "%d errors when parsing input matrix.\n", error_count);      

	if(SvTYPE(SvRV(weight_ref)) == SVt_PVAV) { 
		error_count = malloc_row_perl2c_dbl(aTHX_ weight_ref, weight, nweights_ptr);
		if(error_count > 0 || *nweights_ptr != nweights) {
			Perl_warn(aTHX_ "Weight array has %d items, should have %d. "
				"%d errors detected.\n", *nweights_ptr, nweights, error_count);      
			free(*weight);
			*weight = malloc_row_dbl(aTHX_ nweights,1.0);
			*nweights_ptr = nweights;
		}
	} else {
			*weight = malloc_row_dbl(aTHX_ nweights,1.0);
			*nweights_ptr = nweights;
	}

	return 0;
}

/******************************************************************************/
/**                                                                          **/
/** XS code begins here                                                      **/
/**                                                                          **/
/******************************************************************************/
/******************************************************************************/

MODULE = Algorithm::Cluster	PACKAGE = Algorithm::Cluster
PROTOTYPES: ENABLE


SV *
_hello()
   CODE:
   printf("Hello, world!\n");
	RETVAL = newSVpv("Hello world!!\n", 0);

	OUTPUT:
	RETVAL

int
_readprint(input)
	SV *      input;
	PREINIT:
	int       nrows, ncols;
	double ** matrix;  /* two-dimensional matrix of doubles */

	CODE:
	malloc_matrix_perl2c_dbl(aTHX_ input, &matrix, &nrows, &ncols, NULL);

	if(matrix != NULL) {
		print_matrix_dbl(aTHX_ matrix,nrows,ncols);
		free_matrix_dbl(matrix,nrows);
		RETVAL = 1;
	} else {
		RETVAL = 0;
	}

	OUTPUT:
	RETVAL


SV *
_readformat(input)
	SV *      input;
	PREINIT:
	int       nrows, ncols;
	double ** matrix;  /* two-dimensional matrix of doubles */

	CODE:
	malloc_matrix_perl2c_dbl(aTHX_ input, &matrix, &nrows, &ncols, NULL);

	if(matrix != NULL) {
		RETVAL = format_matrix_dbl(aTHX_ matrix,nrows,ncols);
		free_matrix_dbl(matrix,nrows);
	} else {
		RETVAL = newSVpv("",0);
	}

	OUTPUT:
	RETVAL


SV *
_mean(input)
	SV * input;

	PREINIT:
	int array_length;
	double * data;  /* one-dimensinal array of doubles */

	CODE:
	if(SvTYPE(SvRV(input)) != SVt_PVAV) { 
		XSRETURN_UNDEF;
	}

	malloc_row_perl2c_dbl (aTHX_ input, &data, &array_length);

	RETVAL = newSVnv( mean(array_length, data) );

	OUTPUT:
	RETVAL


SV *
_median(input)
	SV * input;

	PREINIT:
	int array_length;
	double * data;  /* one-dimensinal array of doubles */

	CODE:
	if(SvTYPE(SvRV(input)) != SVt_PVAV) { 
		XSRETURN_UNDEF;
	}

	malloc_row_perl2c_dbl (aTHX_ input, &data, &array_length);

	RETVAL = newSVnv( median(array_length, data) );

	free(data);

	OUTPUT:
	RETVAL


void
_treecluster(nrows,ncols,data_ref,mask_ref,weight_ref,applyscale,transpose,dist,method)
    int      nrows;
    int      ncols;
    SV *     data_ref;
    SV *     mask_ref;
    SV *     weight_ref;
    int      applyscale;
    int      transpose;
    char *   dist;
    char *   method;

    PREINIT:
    SV   *    result_ref;
    SV   *    linkdist_ref;
    int       (*result)[2];
    double   * linkdist;
    int       nweights;

    double  * weight;
    double ** matrix;
    int    ** mask;

    PPCODE:
    /* ------------------------
     * Don't check the parameters, because we rely on the Perl
     * caller to check most paramters.
     */

    /* ------------------------
     * Malloc space for result[][2] and linkdist[]. 
     * Don't bother to cast the pointer for 'result', because we can't 
     * cast it to a pointer-to-array anyway. 
     */
    if (transpose==0) {
	nweights = ncols;
	result   = malloc(2 * (nrows-1) * sizeof(int) );
	linkdist = malloc(    (nrows-1) * sizeof(double) );
    } else {
	nweights = nrows;
	result   = malloc(2 * (ncols-1) * sizeof(int) );
	linkdist = malloc(    (ncols-1) * sizeof(double) );
    }

    /* ------------------------
     * Convert data and mask matrices and the weight array
     * from C to Perl.  Also check for errors, and ignore the
     * mask or the weight array if there are any errors. 
     * Set nweights to the correct number of weights.
     */
    malloc_matrices( aTHX_
	weight_ref, &weight, &nweights, 
	data_ref,   &matrix,
	mask_ref,   &mask,  
	nrows,      ncols
    );

    /* ------------------------
     * Run the library function
     */
    treecluster( nrows, ncols, matrix, mask, weight, applyscale, 
		transpose, dist[0], method[0], result, linkdist, 0);

    /* ------------------------
     * Check result to make sure we didn't run into memory problems
     */
    if(result[0][0]==0 && result[0][1]==0) {
        /* treecluster failed due to insufficient memory */
	if(warnings_enabled(aTHX))
            Perl_warn(aTHX_ "treecluster failed due to insufficient memory.\n");
    }
    else {

        /* ------------------------
         * Convert generated C matrices to Perl matrices
         */
        if (transpose==0) {
            result_ref   = matrix_c_array_2perl_int(aTHX_ result,   nrows-1, 2);
            linkdist_ref =           row_c2perl_dbl(aTHX_ linkdist, nrows-1   ); 
        } else {
            result_ref   = matrix_c_array_2perl_int(aTHX_ result,   ncols-1, 2);
            linkdist_ref =           row_c2perl_dbl(aTHX_ linkdist, ncols-1   );
        }

        /* ------------------------
         * Push the new Perl matrices onto the return stack
         */
        XPUSHs(sv_2mortal( result_ref   ));
        XPUSHs(sv_2mortal( linkdist_ref ));
    }

    /* ------------------------
     * Free what we've malloc'ed 
     */
    free_matrix_int(mask,     nrows);
    free_matrix_dbl(matrix,   nrows);

    free(weight);
    free(result);
    free(linkdist);

    /* Finished _treecluster() */


void
_kcluster(nclusters,nrows,ncols,data_ref,mask_ref,weight_ref,transpose,npass,method,dist)
	int      nclusters;
	int      nrows;
	int      ncols;
	SV *     data_ref;
	SV *     mask_ref;
	SV *     weight_ref;
	int      transpose;
	int      npass;
	char *   method;
	char *   dist;

	PREINIT:
	SV  *    centroid_ref;
	SV  *    clusterid_ref;
	int *    clusterid;
	int      nweights;
	double   error;
	int      ifound;

	double  * weight;
	double ** matrix;
	int    ** mask;
	double ** centroid;

	PPCODE:
	/* ------------------------
	 * Don't check the parameters, because we rely on the Perl
	 * caller to check most parameters.
	 */

	/* ------------------------
	 * Malloc space for the return values from the library function
	 */
        if (transpose==0) {
		nweights = ncols;
		clusterid = malloc(nrows * sizeof(int) );
		centroid  = malloc_matrix_dbl(aTHX_ nclusters,ncols,0.0);
	} else {
		nweights = nrows;
		clusterid = malloc(ncols * sizeof(int) );
		centroid  = malloc_matrix_dbl(aTHX_ nrows,nclusters,0.0);
	}

	/* ------------------------
	 * Convert data and mask matrices and the weight array
	 * from C to Perl.  Also check for errors, and ignore the
	 * mask or the weight array if there are any errors. 
	 * Set nweights to the correct number of weights.
	 */
	malloc_matrices( aTHX_
		weight_ref, &weight, &nweights, 
		data_ref,   &matrix,
		mask_ref,   &mask,  
		nrows,      ncols
	);

	/* ------------------------
	 * Run the library function
	 */
	kcluster( 
		nclusters, nrows, ncols, 
		matrix, mask, weight,
		transpose, npass, method[0], dist[0], clusterid, 
		centroid, &error, &ifound
	);

	/* ------------------------
	 * Convert generated C matrices to Perl matrices
	 */
        if (transpose==0) {
		clusterid_ref =    row_c2perl_int(aTHX_ clusterid, nrows           );
		centroid_ref  = matrix_c2perl_dbl(aTHX_ centroid,  nclusters, ncols);
	} else {
		clusterid_ref =    row_c2perl_int(aTHX_ clusterid, ncols           );
		centroid_ref  = matrix_c2perl_dbl(aTHX_ centroid,  nrows, nclusters);
	}

	/* ------------------------
	 * Push the new Perl matrices onto the return stack
	 */
	XPUSHs(sv_2mortal( clusterid_ref   ));
	XPUSHs(sv_2mortal( centroid_ref    ));
	XPUSHs(sv_2mortal( newSVnv(error) ));
	XPUSHs(sv_2mortal( newSViv(ifound) ));

	/* ------------------------
	 * Free what we've malloc'ed 
	 */
        if (transpose==0) {
		free_matrix_dbl(centroid, nclusters);
	} else {
		free_matrix_dbl(centroid, nrows);
	}
	free(clusterid);
	free_matrix_int(mask,     nrows);
	free_matrix_dbl(matrix,   nrows);
	free(weight);

	/* Finished _kcluster() */



double
_clusterdistance(nrows,ncols,data_ref,mask_ref,weight_ref,cluster1_len,cluster2_len,cluster1_ref,cluster2_ref,dist,method,transpose)
	int      nrows;
	int      ncols;
	SV *     data_ref;
	SV *     mask_ref;
	SV *     weight_ref;
	int      cluster1_len;
	int      cluster2_len;
	SV *     cluster1_ref;
	SV *     cluster2_ref;
	char *   dist;
	char *   method;
	int      transpose;

	PREINIT:
	int   error_count;
	int   nweights;
	int   dummy;

	int     * cluster1;
	int     * cluster2;

	double  * weight;
	double ** matrix;
	int    ** mask;

	double distance;

	CODE:
	error_count = 0;

	/* ------------------------
	 * Don't check the parameters, because we rely on the Perl
	 * caller to check most paramters.
	 */

	/* ------------------------
	 * Convert cluster index Perl arrays to C arrays
	 */
	error_count += malloc_row_perl2c_int(aTHX_ cluster1_ref, &cluster1, &dummy);
	error_count += malloc_row_perl2c_int(aTHX_ cluster2_ref, &cluster2, &dummy);

	/* ------------------------
	 * Convert data and mask matrices and the weight array
	 * from C to Perl.  Also check for errors, and ignore the
	 * mask or the weight array if there are any errors. 
	 * Set nweights to the correct number of weights.
	 */
	nweights = (transpose==0) ? ncols : nrows;
	malloc_matrices( aTHX_
		weight_ref, &weight, &nweights, 
		data_ref,   &matrix,
		mask_ref,   &mask,  
		nrows,      ncols
	);


	/* ------------------------
	 * Run the library function
	 */
	distance = clusterdistance( 
		nrows, ncols, 
		matrix, mask, weight,
		cluster1_len, cluster2_len, cluster1, cluster2,
		dist[0], method[0], transpose
	);

	RETVAL = distance;

	/* ------------------------
	 * Free what we've malloc'ed 
	 */
	free_matrix_int(mask,     nrows);
	free_matrix_dbl(matrix,   nrows);
	free(weight);
	free(cluster1);
	free(cluster2);

	/* Finished _clusterdistance() */

	OUTPUT:
	RETVAL



void
_somcluster(nrows,ncols,data_ref,mask_ref,weight_ref,transpose,nxgrid,nygrid,inittau,niter,dist)
	int      nrows;
	int      ncols;
	SV *     data_ref;
	SV *     mask_ref;
	SV *     weight_ref;
	int      transpose;
	int      nxgrid;
	int      nygrid;
	double   inittau;
	int      niter;
	char *   dist;

	PREINIT:
	int      (*clusterid)[2];
	SV *  clusterid_ref;
	double*** celldata;
	SV *  celldata_ref;

	double  * weight;
	double ** matrix;
	int    ** mask;
	int       nweights;

	PPCODE:
	/* ------------------------
	 * Don't check the parameters, because we rely on the Perl
	 * caller to check most paramters.
	 */

	/* ------------------------
	 * Allocate space for clusterid[][2]. 
	 * Don't bother to cast the pointer, because we can't cast
	 * it to a pointer-to-array anway. 
	 */
	if (transpose==0) {
		clusterid  =  malloc(2 * (nrows) * sizeof(int) );
	} else {
		clusterid  =  malloc(2 * (ncols) * sizeof(int) );
	}
	celldata  =  0;
	/* Don't return celldata, for now at least */


	/* ------------------------
	 * Convert data and mask matrices and the weight array
	 * from C to Perl.  Also check for errors, and ignore the
	 * mask or the weight array if there are any errors. 
	 * Set nweights to the correct number of weights.
	 */
	nweights = (transpose==0) ? ncols : nrows;
	malloc_matrices( aTHX_
		weight_ref, &weight, &nweights, 
		data_ref,   &matrix,
		mask_ref,   &mask,  
		nrows,      ncols
	);

	/* ------------------------
	 * Run the library function
	 */
	somcluster( 
		nrows, ncols, 
		matrix, mask, weight,
		transpose, nxgrid, nygrid, inittau, niter,
		dist[0], celldata, clusterid
	);

	/* ------------------------
	 * Convert generated C matrices to Perl matrices
	 */
	clusterid_ref = matrix_c_array_2perl_int(aTHX_ clusterid, nrows, 2); 

	/* ------------------------
	 * Push the new Perl matrices onto the return stack
	 */
	XPUSHs(sv_2mortal( clusterid_ref   ));

	/* ------------------------
	 * Free what we've malloc'ed 
	 */
	free_matrix_int(mask,     nrows);
	free_matrix_dbl(matrix,   nrows);
	free(weight);
	free(clusterid);

	/* Finished _somcluster() */