hddc_toolbox/0000755000175000001440000000000010673440706013075 5ustar charlesusershddc_toolbox/source/0000755000175000001440000000000010526622635014375 5ustar charlesusershddc_toolbox/source/hddc_learn_fast.mexglx.gcc-4.x0000755000175000001440000004045410526621723022074 0ustar charlesusersELF 414 (+++;;LP,<<Qtd%-+(',%!$ # " * &)H   ) * +;;;0=8=_93 i g)  PA&n?:_Y,<QX'Lg3f! ZC"5==2uC0f `| __gmon_start____cxa_finalize_Jv_RegisterClassesmxGetStringmexErrMsgTxtsqrtdsyrk_memcpydsyev_dsyevx_dgemm_memsetlogexprandommxMallocmexPrintfmxFreemexFunctionmxGetMmxGetNmxGetPrmxGetScalarmxGetClassIDmxCreateStructMatrixmxCreateDoubleMatrixmxSetFieldByNumbermxCreateNumericMatrixmxGetDatamxCreateCellMatrixmxSetCellgettimeofdaymexVersionlibmx.solibmex.solibmat.solibm.so.6libstdc++.so.6libpthread.so.0libc.so.6hddc_learn_fast.mexglxMEXGLIBC_2.0v7.0GLIBC_2.1.3H6Q ii { 0si ii r%JbiL;  !! ""a"m""" #<#####$D$t$$$$%%%$%.%)))+0=4="*!!""$"#&)n ("?&:&Q&d+++++*+ +Y[' Z!e)!$!!a$$%%&&'')! &# Y# # .$ " T%#%%y&&''(%*%*M&*&*&*4'**&'''U&' '%<))!="=)=,(=",=,U=S hhUS[À0tX[ÐUSS0u8t$} ҋuƃ[]Í&US./tt $҃[]Ë$ÐUWVS E}~|EE}ǾEE} ~2EEvى؃U9M u؍SE9tE9uuE [^_]USTD$2EʉD$E$t $*UʻtIMʸat*ft btdt p Qu؃T[]UWVS|MEݝk]#]}8uEUU8 ‰MUʍuƉx|؉]Ép}<up}< ljppӉ]EÉ]؉EuƉuƉt}(u u(0t}$u"t}$ljt EMMMM]É8E8ߋ@Dzt$\؃;E|ǃ}>M 9E 9]~ }}?݅U‹ ~ظ@ 9u݅0ʋ< ዕpڃ9]؋]~vEE<pNjp<EÍ  ]Ã]9щ؉ڃE;U~~݅0=+]؍ut$$}|$ U؉T$ED$ML$T$D$ t$D$p*${*M~؍E9ul}ED$ED$x$ED$ ED$|L$8\$ut$x|$ t$D${*$*E0E=+Mȉ~+84ʺFfv؃9uރx_)܍8Ddvȃ`v؉֍F))؋8L)9uظM],U)ƉEUEE)ԍL$_E)čD$_U؉T$LD$HL$DED$@||$?MbP?MATLAB R14 native+r{ H )  =8 Ho(oo`ooo;<v 0=;GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux),)H &$)] |/usr/src/packages/BUILD/glibc-2.4/cc-nptl/csu/crti.S/usr/src/packages/BUILD/glibc-2.4/csuGNU AS 2.16.91.0.5|/usr/src/packages/BUILD/glibc-2.4/cc-nptl/csu/crtn.S/usr/src/packages/BUILD/glibc-2.4/csuGNU AS 2.16.91.0.5%%K /usr/src/packages/BUILD/glibc-2.4/cc-nptl/csucrti.S)3!/!=Z!H #!/=  !/!=Z!gg//Z!!!sK /usr/src/packages/BUILD/glibc-2.4/cc-nptl/csucrtn.S)!!!]  !.symtab.strtab.shstrtab.hash.dynsym.dynstr.gnu.version.gnu.version_d.gnu.version_r.rel.dyn.rel.plt.init.text.fini.rodata.eh_frame.ctors.dtors.jcr.dynamic.got.got.plt.data.bss.comment.debug_aranges.debug_info.debug_abbrev.debug_lineP!  )1oZ>o((8Mo``\ He 8 8  nH H i` ` 0t Tz))**++;+;+;+<,=- =-0=0-8=8-8--XP.P/ p/0`6pE <\(`8 H  `   ) * +;;;<==0=8= A  Q\;j;x;8=4=   Q;;+;) <N>  ]j+ r8. z   <`  0=H  8=t)/ ) 8=W  =<= '9@9S cgiy )  A&:_Y,<&'8gHZf! fC~"5=2C0  %`>M initfini.c/usr/src/packages/BUILD/glibc-2.4/cc-nptl/csu/crti.Scall_gmon_startcrtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST__completed.5751p.5749__do_global_dtors_auxframe_dummy__CTOR_END____DTOR_END____FRAME_END____JCR_END____do_global_ctors_aux/usr/src/packages/BUILD/glibc-2.4/cc-nptl/csu/crtn.Shddc_learn_fast.ccompute_paramsmexversion.cversionem_hddcinterpret_commonrandsample_DYNAMICcompute_class__dso_handle_init__bss_startgetMilliSecs_fini_edata__i686.get_pc_thunk.bx_GLOBAL_OFFSET_TABLE__endcompute_probasmxGetScalar@@v7.0dsyev_mxGetClassID@@v7.0mxGetData@@v7.0MEXrandom@@GLIBC_2.0mxGetPr@@v7.0mexVersionmxSetCell@@v7.0mxCreateDoubleMatrix@@v7.0dgemm_mexErrMsgTxt@@v7.0mxMalloc@@v7.0mxGetN@@v7.0exp@@GLIBC_2.0mexPrintf@@v7.0mxFree@@v7.0dsyrk_memcpy@@GLIBC_2.0sqrt@@GLIBC_2.0mxGetString@@v7.0mexFunctiongettimeofday@@GLIBC_2.0__cxa_finalize@@GLIBC_2.1.3mxCreateStructMatrix@@v7.0mxGetM@@v7.0mxCreateCellMatrix@@v7.0memset@@GLIBC_2.0mxCreateNumericMatrix@@v7.0dsyevx__Jv_RegisterClassesmxSetFieldByNumber@@v7.0log@@GLIBC_2.0__gmon_start__hddc_toolbox/source/hddc_learn_fast.c0000600000175000001440000004610310375105645017634 0ustar charlesusers#include #include #include "mex.h" #include #include /* // 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; } #define TICXX lines[n_t]=__LINE__; times[n_t++]+=getMilliSecs() #define TIC /******************************************************** * Subfunctions */ #define MAXD ((p)/2) /* 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 #define Free_a 16 /* 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;ival) { ind=j; val=sc_j; } } for(j=(MAXD-1);j>-1;j--){ sc_j=LL(j)-LL(j+1); if(sc_j>seuil*val) { d_i=j+1; break; } } #undef LL } //if(d_i==0) //printf("In em_hdda_c_fast: d_i is %d\n",d_i); //if(dims) dims[i]=d_i; //TIC; { /* compute corresponding eigenvecs */ int il=p-d_i+1,iu=p; /* eigenvecs nos il to iu (ascending, inclusive) */ double dminusone=-1,dunused; int iunused; int iwork[5*p]; int ifail[n]; int info; dsyevx_("Vectors also","Indexed values","Upper", &p,sigma_i,&p,&dunused,&dunused,&il,&iu,&dminusone, &iunused,eigenvals,eigenvecs,&p,work_eig,&lwork_eig,iwork, ifail,&info); ASSERT(info==0,"dsyevx_ error"); } //TIC; } { /* * %%%%%%%%%%%%%%%%%%%% Compute estimators %%%%%%%%%%%%%%%%%%%% * all models handled in 2 loops */ double sum_n_s1,sum_n_s2,sum_n_pd,sum_n_d; sum_n_s1=sum_n_s2=sum_n_pd=sum_n_d=0; for(i=0;imaxpow) maxpow=powi; clusterprobas[i]=exp(powi); accu=accu+clusterprobas[i]; } if(maxpow>limitpow || accu==0 || accu>=HUGE_VAL){ for(i=0;ival) { val=T(j,l); ind=l; } if(cls[j]!=ind+1) { n_chg++; cls[j]=ind+1; } } return n_chg; } /* the main algorithm * (see previous functions for explanation of parameters) */ void em_hddc(double *data,int k,int n,int p,double *t,double *s, double diff,int maxiter,double seuil,int common, int *cls,double *rand_in, double *init, double *a,double *b, int *dims,int dim, double *prop, double *m,double *w) { int i,j; int it; double q_it_prev=0; //double *s=NEWA(double,n*k); int lwork=compute_params_lwork(k,n,p); double *work=NEWA(double,lwork); //double *ktable=NEWA(double,k); double times[50]; //int lines[50]; //int n_t=0; for(i=0;i<50;i++) times[i]=0; { /* initialization */ if(init){ for(i=0;iHUGE_VAL: %d \n",exp(1000.0)>=HUGE_VAL); printf("Check overflow Inf0 && (fabs(q_it-q_it_prev)?TransposedUpperNo vectorsIndexed valuesVectors alsoNot transposed* User initialization ! * Random initialization ! * Iteration %d / %d abpropmW3 output arguments required.variable name requireddiffmaxiterseuilrand_inputcommoninitunknown variable nameCould not convert common's string data.even nb of input arguments required.Could not convert string data.+&-&j&/&4&6&MATLAB R14 nativekt~ 0%  t8 o ooDooo9`7 8X74'GCC: (GNU) 3.3.3 20040412 (Red Hat Linux 3.3.3-7)GCC: (GNU) 3.3.3 20040412 (Red Hat Linux 3.3.3-7)GCC: (GNU) 3.3.3 20040412 (Red Hat Linux 3.3.3-7)GCC: (GNU) 3.3.3 20040412 (Red Hat Linux 3.3.3-7)GCC: (GNU) 3.3.3 20040412 (Red Hat Linux 3.3.3-7)GCC: (GNU) 3.3.3 20040412 (Red Hat Linux 3.3.3-7).symtab.strtab.shstrtab.hash.dynsym.dynstr.gnu.version.gnu.version_d.gnu.version_r.rel.dyn.rel.plt.init.text.fini.rodata.eh_frame.ctors.dtors.jcr.dynamic.got.got.plt.data.bss.commentL!  )1oX>o  8MoDD\ e   n i  t$ $ z0%0%P%P%H'H'L7L'T7T'\7\'`7`'h8h( t8t(8( 8((2).= 4 D     $  0% P% H'L7T7\7`7h8t888$  L7*T78\7E8I8UH  k  wP7X7H'\7$  K 8  8 R w %  0`79v G8T  Z8f ' s0% y8t881 9 a  $  A%&@G:Z_iYv,<y[ :" 5&=32LC^0z ` call_gmon_startcrtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST__p.0completed.1__do_global_dtors_auxframe_dummy__CTOR_END____DTOR_END____FRAME_END____JCR_END____do_global_ctors_auxhddc_learn_fast.ccompute_params_lworkcompute_paramsmexversion.cversionem_hddcinterpret_commonrandsample_DYNAMICcompute_class__dso_handle_init__bss_startgetMilliSecs_fini_edata_GLOBAL_OFFSET_TABLE__endcompute_probasmxGetScalar@@v7.0dsyev_mxGetClassID@@v7.0mxGetData@@v7.0MEXrandom@@GLIBC_2.0mxGetPr@@v7.0mexVersionmxSetCell@@v7.0mxCreateDoubleMatrix@@v7.0dgemm_mexErrMsgTxt@@v7.0mxMalloc@@v7.0mxGetN@@v7.0exp@@GLIBC_2.0mexPrintf@@v7.0mxFree@@v7.0dsyrk_sqrt@@GLIBC_2.0mxGetString@@v7.0mexFunctiongettimeofday@@GLIBC_2.0__cxa_finalize@@GLIBC_2.1.3mxCreateStructMatrix@@v7.0mxGetM@@v7.0mxCreateCellMatrix@@v7.0memset@@GLIBC_2.0mxCreateNumericMatrix@@v7.0dsyevx__Jv_RegisterClassesmxSetFieldByNumber@@v7.0log@@GLIBC_2.0__gmon_start__hddc_toolbox/source/hddc_learn.c0000644000175000001440000004110010375105560016613 0ustar charlesusers#include #include #include "mex.h" #include #include /* // 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;ival) { ind=j; val=sc_j; } if(j>ind && sc_j=p/2) { d_i=p/2; } #undef LL } if(dims) dims[i]=d_i; TIC; { /* compute corresponding eigenvecs */ int il=p-d_i+1,iu=p; /* eigenvecs nos il to iu (ascending, inclusive) */ double dminusone=-1,dunused; int iunused; int iwork[5*p]; int ifail[n]; int info; dsyevx_("Vectors also","Indexed values","Upper", &p,sigma_i,&p,&dunused,&dunused,&il,&iu,&dminusone, &iunused,eigenvals,eigenvecs,&p,work_eig,&lwork_eig,iwork, ifail,&info); ASSERT(info==0,"dsyevx_ error"); } TIC; } { /* * %%%%%%%%%%%%%%%%%%%% Compute estimators %%%%%%%%%%%%%%%%%%%% * all models handled in 2 loops */ double sum_n_s1,sum_n_s2,sum_n_pd,sum_n_d; sum_n_s1=sum_n_s2=sum_n_pd=sum_n_d=0; for(i=0;ival) { val=T(j,l); ind=l; } if(cls[j]!=ind+1) { n_chg++; cls[j]=ind+1; } } return n_chg; } /* the main algorithm * (see previous functions for explanation of parameters) */ void em_hddc(double *data,int k,int n,int p,double *t,double *s, double diff,int maxiter,double seuil,int common, int *cls,double *rand_in, double *init, double *a,double *b, int *dims,int dim, double *prop, double *m,double *w) { int i,j; int it; double q_it_prev=0; //double *s=NEWA(double,n*k); int lwork=compute_params_lwork(k,n,p); double *work=NEWA(double,lwork); { /* initialization */ if(init){ for(i=0;i0 && (fabs(q_it-q_it_prev)ߋ@Dzt$\؃;E|ǃ}>M 9E 9]~ }}?݅U‹ ~ظ@ 9u݅0ʋ< ዕpڃ9]؋]~vEE<pNjp<EÍ  ]Ã]9щ؉ڃE;U~~݅0=+]؍ut$$}|$ U؉T$ED$ML$T$D$ t$D$p*${*M~؍E9ul}ED$ED$x$ED$ ED$|L$8\$ut$x|$ t$D${*$*E0E=+Mȉ~+84ʺFfv؃9uރx_)܍8Ddvȃ`v؉֍F))؋8L)9uظM],U)ƉEUEE)ԍL$_E)čD$_U؉T$LD$HL$DED$@||$?MbP?MATLAB R14 native+r{ H )  =8 Ho(oo`ooo;<v 0=;GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux)GCC: (GNU) 4.1.0 (SUSE Linux),)H &$)] |/usr/src/packages/BUILD/glibc-2.4/cc-nptl/csu/crti.S/usr/src/packages/BUILD/glibc-2.4/csuGNU AS 2.16.91.0.5|/usr/src/packages/BUILD/glibc-2.4/cc-nptl/csu/crtn.S/usr/src/packages/BUILD/glibc-2.4/csuGNU AS 2.16.91.0.5%%K /usr/src/packages/BUILD/glibc-2.4/cc-nptl/csucrti.S)3!/!=Z!H #!/=  !/!=Z!gg//Z!!!sK /usr/src/packages/BUILD/glibc-2.4/cc-nptl/csucrtn.S)!!!]  !.symtab.strtab.shstrtab.hash.dynsym.dynstr.gnu.version.gnu.version_d.gnu.version_r.rel.dyn.rel.plt.init.text.fini.rodata.eh_frame.ctors.dtors.jcr.dynamic.got.got.plt.data.bss.comment.debug_aranges.debug_info.debug_abbrev.debug_lineP!  )1oZ>o((8Mo``\ He 8 8  nH H i` ` 0t Tz))**++;+;+;+<,=- =-0=0-8=8-8--XP.P/ p/0`6pE <\(`8 H  `   ) * +;;;<==0=8= A  Q\;j;x;8=4=   Q;;+;) <N>  ]j+ r8. z   <`  0=H  8=t)/ ) 8=W  =<= '9@9S cgiy )  A&:_Y,<&'8gHZf! fC~"5=2C0  %`>M initfini.c/usr/src/packages/BUILD/glibc-2.4/cc-nptl/csu/crti.Scall_gmon_startcrtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST__completed.5751p.5749__do_global_dtors_auxframe_dummy__CTOR_END____DTOR_END____FRAME_END____JCR_END____do_global_ctors_aux/usr/src/packages/BUILD/glibc-2.4/cc-nptl/csu/crtn.Shddc_learn_fast.ccompute_paramsmexversion.cversionem_hddcinterpret_commonrandsample_DYNAMICcompute_class__dso_handle_init__bss_startgetMilliSecs_fini_edata__i686.get_pc_thunk.bx_GLOBAL_OFFSET_TABLE__endcompute_probasmxGetScalar@@v7.0dsyev_mxGetClassID@@v7.0mxGetData@@v7.0MEXrandom@@GLIBC_2.0mxGetPr@@v7.0mexVersionmxSetCell@@v7.0mxCreateDoubleMatrix@@v7.0dgemm_mexErrMsgTxt@@v7.0mxMalloc@@v7.0mxGetN@@v7.0exp@@GLIBC_2.0mexPrintf@@v7.0mxFree@@v7.0dsyrk_memcpy@@GLIBC_2.0sqrt@@GLIBC_2.0mxGetString@@v7.0mexFunctiongettimeofday@@GLIBC_2.0__cxa_finalize@@GLIBC_2.1.3mxCreateStructMatrix@@v7.0mxGetM@@v7.0mxCreateCellMatrix@@v7.0memset@@GLIBC_2.0mxCreateNumericMatrix@@v7.0dsyevx__Jv_RegisterClassesmxSetFieldByNumber@@v7.0log@@GLIBC_2.0__gmon_start__hddc_toolbox/hddc_learn.m0000600000175000001440000001455610526631526015340 0ustar charlesusersfunction [prms,T,cls] = hddc_learn(XX,k,varargin); % High Dimensionality Data Clustering (hddc) % % Usage: (1) [prms,T,cls] = hddc_learn(X,k,'model','AiBiQiDi'); % (2) [prms,T,cls] = hddc_learn(X,k,'model','best'); % (3) [prms,T,cls] = hddc_learn(X,k,'model','AiBiQiDi','seuil',0.1); % (4) [prms,T,cls] = hddc_learn(X,k,'model','AiBiQiDi','maxiter',100,); % % Input: % - X: data % - k: number of classes % - model: % % Output: % - cls: label of test data Y.data, % - T: a posteriori probabilities % - prms: parameters of classes % % Authors: C. Bouveyron - 2004-2006 % % Reference: C. Bouveyron, S. Girard and C. Schmid, High-Dimensional Data % Clustering, Technical Report n°1083M, LMC-IMAG, Université J. Fourier % Grenoble 1. % Initialization of parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if isstruct(XX), X = XX.data; if exist('XX.cls'), class = XX.cls; end else X = XX; class = []; end [N,p] = size(X); maxiter = 100; seuil = 0.2; model = 'AiBiQiDi'; common = ''; init = 'random'; dim = []; diff = 1e-5; %%% PARAMETERS MANAGEMENT for i=1:2:length(varargin) if ~isstr(varargin{i}) | ~exist(varargin{i},'var') error(['Invalid parameter ' varargin{i}]); end eval([varargin{i} '= varargin{i+1};']); end %%% if isempty(strmatch(model,strvcat('AijBiQiDi', 'AijBQiDi', 'AiBiQiDi', 'ABiQiDi',... 'AiBQiDi', 'ABQiDi','AijBiQiD', 'AijBQiD', 'AiBiQiD', 'ABiQiD', 'AiBQiD',... 'ABQiD','best'),'exact')) error('--> The parameter ''model'' is not valide: see the help of hddc.'); end % Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if ~ischar(init), init_cls = init; fprintf('* User initialization !\n'); else switch init case 'random' init_cls = randsample(k,N,1); fprintf('* Random initialization !\n'); case 'kmeans' init_cls = kmeans(X,k); fprintf('* Initialization by kmeans !\n'); end end % Calling the main function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if isequal(model,'best') % looking for the model with the smallest BIC value if isempty(dim), % models with free dimensions models = {'AijBiQiDi', 'AijBQiDi', 'AiBiQiDi', 'ABiQiDi', 'AiBQiDi', 'ABQiDi'}; fprintf('--> Bic values of the HDDA models with free dimensions:\n'); for i = 1:length(models) [prms_t{i},T_t{i},cls_t{i}] = learn(X,k,models{i},seuil,maxiter,dim,init_cls,diff); bic(i) = prms_t{i}.bic; fprintf(' - model %s: %g\n',models{i},bic(i)); end else % models with common dimensions models = {'AijBiQiD', 'AijBQiD', 'AiBiQiD', 'ABiQiD', 'AiBQiD', 'ABQiD'}; fprintf('--> Bic values of the HDDA models with common dimensions:\n'); for i = 1:length(models) [prms_t{i},T_t{i},cls_t{i}] = learn(X,k,models{i},seuil,maxiter,dim,init_cls,diff); bic(i) = prms_t{i}.bic; fprintf(' - model %s: %g\n',models{i},bic(i)); end end [val,ind] = min(fliplr(bic)); prms = prms_t{ind}; T = T_t{ind}; cls = cls_t{ind}; fprintf('--> Best model: %s\n',prms_t{ind}.model); else [prms,T,cls] = learn(X,k,model,seuil,maxiter,dim,init_cls,diff); fprintf('--> Used model: %s\n',prms.model); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SUB-FUNCTION % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [prms,T,cls] = learn(X,k,model,seuil,maxiter,dim,init_cls,diff); [N,p] = size(X); % Convert model name in the old format %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% switch model case 'AijBiQiDi' % Model [a_ij b_i Q_i d_i] common = 'f'; case 'AijBiQiD' % Model [a_ij b_i Q_i d] common = 'fd'; case 'AijBQiDi' % Model [a_ij b_i Q_i d_i] common = 'fb'; case 'AijBQiD' % Model [a_ij b_i Q_i d] common = 'fbd'; case 'AiBiQiDi' % Model [a_i b_i Q_i d_i] common = ''; case 'AiBiQiD' % Model [a_i b_i Q_i d] common = 'd'; case 'ABiQiDi' % Model [a b_i Q_i d_i] common = 'a'; case 'ABiQiD' % Model [a b_i Q_i d] common = 'ad'; case 'AiBQiDi' % Model [a_i b Q_i d_i] common = 'b'; case 'AiBQiD' % Model [a_i b Q_i d] common = 'bd'; case 'ABQiDi' % Model [a b Q_i d_i] common = 'ab'; case 'ABQiD' % Model [a b Q_i d] common = 'abd'; otherwise error('Not yet implemented!'); end % Calling the C function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [prms,T,cls,S] = hddc_learn_fast(X,k,'seuil',seuil,'maxiter',maxiter,... 'common',common,'dim',dim,'init',init_cls,'diff',diff); % Compute the BIC value %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% d = double(prms.d); switch model case 'AijBiQiDi' % Model [a_ij b_i Q_i d_i] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + sum(d) + 2*k; case 'AijBiQiD' % Model [a_ij b_i Q_i d] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + sum(d) + k + 1; case 'AijBQiDi' % Model [a_ij b_i Q_i d_i] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + sum(d) + k + 1; case 'AijBQiD' % Model [a_ij b_i Q_i d] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + sum(d) + 2; case 'AiBiQiDi' % Model [a_i b_i Q_i d_i] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + 3*k; case 'AiBiQiD' % Model [a_i b_i Q_i d] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + 2*k + 1; case 'ABiQiDi' % Model [a b_i Q_i d_i] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + 2*k + 1; case 'ABiQiD' % Model [a b_i Q_i d] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + k + 2; case 'AiBQiDi' % Model [a_i b Q_i d_i] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + 2*k + 1; case 'AiBQiD' % Model [a_i b Q_i d] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + k + 2; case 'ABQiDi' % Model [a b Q_i d_i] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + k + 2; case 'ABQiD' % Model [a b Q_i d] q = (k*p+k-1) + sum( d.*(p-(d+1)/2) ) + 3; otherwise error('Not yet implemented!'); end K = S + p * log(2*pi); ll = sum(sum(T .* K)); bic = (ll + q * log(N)) / N; % Return parameters and results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% prms.p = p; prms.model = model; prms.bic = bic; prms.k = k; prms.d = double(prms.d); prms.m = prms.m'; cls = double(cls);hddc_toolbox/hddc_demo.m0000644000175000001440000000163010526635611015157 0ustar charlesusers% High Dimensionality Data Clustering (demonstration) % % Authors: C. Bouveyron - 2004-2006 % % Reference: C. Bouveyron, S. Girard and C. Schmid, High-Dimensional Data % Clustering, Technical Report n°1083M, LMC-IMAG, Université J. Fourier % Grenoble 1. fprintf('\n- Loading of data\n') fprintf('load(''crabs.mat'');\n') load('crabs.mat' ); fprintf('\n- Press the keyboard to continue ...\n\n'), pause() fprintf('- Looking for the model with the smallest BIC value (k=4 groups)\n') fprintf('[prms,T,cls] = hddc_learn(X,4,''model'',''best'',''seuil'',0.2);\n') [prms,T,cls] = hddc_learn(X,4,'model','best','seuil',0.2); fprintf('\n- Press the keyboard to continue ...\n\n'), pause() fprintf('- Learning parameters of a specific model\n') fprintf('[prms,T,cls] = hddc_learn(X,4,''model'',''AiBQiDi'',''seuil'',0.2)\n'); [prms,T,cls] = hddc_learn(X,4,'model','AiBQiDi','seuil',0.2);hddc_toolbox/crabs.mat0000644000175000001440000000416610412200033014653 0ustar charlesusersMATLAB 5.0 MAT-file, Platform: GLNX86, Created on: Fri Mar 24 13:58:09 2006 IMx|YIRcG;lGGv(B f1|M:Yv G\@>@_  p6^ZLA$U+(~=+/EˣDZٷ> N+^xL_P7S_ty)!\t*N1FFENy-r|(Ɨ8^'F?A{B~TVbUV7iޚQ'ZgIM'`[#;89Ο3{AU/D1 ԫ}Ύs;b=Q| A~ҮXgy+s*?]=~Cٷic|t E_|Г?ORMRp˸d-I?:뜷N?^(ݠ供;+[&,O,kR8vZkVg"_K#T2Ύd{^˖zp)2uq_Hz oG;)ƁCϞ^e_qr|^uwb *57sun 3ҹ7L_u+*^a+ YAc#ɱޚ̓: PP ՛muDwߢ/79F#x,6|.wm>w|yw_m;:ixߒ\O%QuwyހC.#M7oOslq|nW|!.{6}vjZfwk)\)֓}K>]z7QK{h 6F{Ym%m_Ex8qp}9<8˟u~x{wIoԉ~٥qE]3n:*o!<>q+A)/|RsjXw}^u=Sޥ5g|2^&}5Y쫼e:S=P(y:'f'!w^>}R۰/D]'Q)WqRe :4$ +c^3z TDo^c߫a~>SA'{"{;׺qֿ>Gbl瑞]x8"C~vO|烿{]PݟէIO5Ys[8F?RWEq|_59ܟZ3W?F?^w;_vS=4*;b\k{[ c$0 I,$!8hddc_toolbox/hddc_classif.m0000755000175000001440000000474110517663514015673 0ustar charlesusersfunction [cls,P] = hdda_classif2_faster(prms,Y,varargin); % High Dimensionality Discriminant Analysis (classification) % % Usage: (1) [cls,P] = hdda_classif(prms,Y); % % Input: % - prms: parameters of classes % - Y: data to classify % % Output: % - cls: label of test data Y.data, % - P: a posteriori probabilities % % Authors: C. Bouveyron - 2004-2006 % % Reference: C. Bouveyron, S. Girard and C. Schmid, "High Dimensional Discriminant Analysis", % Communications in Statistics, Theory and methods, in press, 2007. % Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if isstruct(Y), data = Y.data; if isfield(Y,'cls'), class = Y.cls; else class = []; end else data = Y; class = []; end [nt,pt] = size(data); model = prms.model; k = prms.k; p = prms.p; a = prms.a; b = prms.b; d = prms.d; prop = prms.prop; m = prms.m; Q = prms.Q; % Cost function computing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:k % Projection of test data in the eigenspace Ei Pa = ((data - repmat(m(i,:),nt,1)) * Q{i}) * Q{i}'; Pb = Pa + repmat(m(i,:),nt,1) - data; % Compute cost function K_i(x) switch model case {'AijBiQiDi','AijBiQiD','AijBQiDi','AijBQiD','AjBiQiD','AjBQiD','AjBQD'} ai = a(1:d(i),i)'; K(:,i) = diag(Pa * Q{i} * diag(1./ai) * Q{i}' * Pa') ... + (1/b(i) * sum(Pb.^2,2)) + sum(log(ai)) ... + (p-d(i)) * log(b(i)) - 2 * log(prop(i)) + p * log(2*pi); otherwise K(:,i) = 1/a(i) * sum(Pa.^2,2) + (1/b(i) * sum(Pb.^2,2)) + d(i) * log(a(i)) ... + (p-d(i)) * log(b(i)) - 2 * log(prop(i)) + p * log(2*pi); end end % a posteriori probabilities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% P = exp(-0.5*K) ./ repmat(sum(exp(-0.5*K),2),1,k); if length(find(isnan(P)))~=0 | length(find(isinf(P)))~=0 fprintf('--> overflow switching to O(k*k) mode\n'); clear P; eksp = zeros(nt,k); for i=1:k for l=1:k eksps(:,l) = exp(0.5*(K(:,i) - K(:,l))); end P(:,i) = 1./sum(eksps,2); end end % Classification %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [val,cls] = max(P,[],2); % Classification results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if ~isempty(class), tx = sum(cls == class) / length(cls); fprintf('--> Correct classification rate: %g \n',tx); end