hddc_learn.c

高维数据分类程序！数据挖掘及图像分析利器！matlab版！

#include <assert.h>
#include <math.h>
#include "mex.h"
#include <sys/time.h>


#include <string.h>

/* 

// export PATH=$PATH:/softs/stow/matlab-14/bin
setenv PATH ${PATH}:/softs/stow/matlab-14/bin
mex CFLAGS=-Wall em_hdda_c.c
cp em_hdda_c.mexglx ..
*/


/********************************************************
 * prototypes for BLAS and LAPACK functions
 */

#define doublereal double
#define integer int

extern int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
 integer *lda, doublereal *w, doublereal *work, integer *lwork,
 integer *info);

extern int dsyevx_(char *jobz, char *range, char *uplo, integer *n,
  doublereal *a, integer *lda, doublereal *vl, doublereal *vu, integer *
  il, integer *iu, doublereal *abstol, integer *m, doublereal *w,
  doublereal *z__, integer *ldz, doublereal *work, integer *lwork,
  integer *iwork, integer *ifail, integer *info);

extern int 
dgemm_(char* transA, char* transB, integer* M, integer* N, integer* K,
       doublereal* alpha,
       doublereal* A, integer* lda,
       doublereal* B, integer* ldb,
       doublereal* beta,
       doublereal* C, integer* ldc);


extern int
dsyrk_(char* uplo, char* trans, integer* N, integer* K,
       doublereal* alpha,
       doublereal* A, integer* lda,
       doublereal* beta,
       doublereal* C, integer* ldc);


#undef doublereal
#undef integer

/********************************************************
 * some useful definitions
 */

/*
#define NEWA(type,n) (type*)mxMalloc((printf("mall %d %d %d\n",n,sizeof(type),__LINE__),sizeof(type)*(n)))
*/

#define NEWA(type,n) (type*)mxMalloc(sizeof(type)*(n))

#define FREE(p) mxFree(p)

#define ASSERT(v,comment) mxAssert(v,comment) 

double getMilliSecs() {
  struct timeval tv;
  gettimeofday(&tv,NULL);
  return tv.tv_sec*1e3+tv.tv_usec*1e-3;
}



/********************************************************
 * Subfunctions
 */

/* common flags
 * any combination of the 3 is acceptable */
#define Common_a 1
#define Common_b 2
#define Common_d 4
#define Common_p 8


/* samples k values among 0..n-1 with replacement 
 * result in cls(k)
 */
void randsample(int n,int k,int *cls,double *rand_input) {
  int i;
  if(rand_input) 
    for(i=0;i<k;i++) 
      cls[i]=ceil(rand_input[i]*n)-1;
  else 
    for(i=0;i<k;i++) 
      cls[i]=random()%n;
}




/* access macros for s and t (the same for all the functions) */

#define T(i,j) t[(i)+(j)*n]
#define S(i,j) s[(i)+(j)*n]


/* lwork for compute_params*/
static int compute_params_lwork(int k,int n,int p) {
  int lwork_eig=(32+3)*p;
  
  return p*k
    +p*p
    +p*p
    +lwork_eig
    +p*k
    +p*p*k
    +p
    +2*n
    +2*p*n
    +k
    +k
    +k
    +k;
}





/* computes intrinsic dimensions and estimator parameters 
 * all matrices are stored by columns (