00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include <string.h>
00008 #include "mrilib.h"
00009 #include <stdlib.h>
00010 #include <ctype.h>
00011 
00012 /*-------------------------- global data --------------------------*/
00013 
00014 /** inputs **/
00015 
00016 static int PC_dmean      = 0 ; /* default is not to remove means */
00017 static int PC_vmean      = 0 ;
00018 static int PC_vnorm      = 0 ; /* 07 July 1999 */
00019 static int PC_normalize  = 0 ; /* and not to normalize */
00020 static int PC_lprin_save = 0 ; /* # of principal components to save */
00021 static int PC_be_quiet   = 1 ; /* quiet is the default */
00022 static int PC_do_float   = 0 ; /* shorts are the default */
00023 
00024 static char ** PC_dsname = NULL ; /* dataset names */
00025 static int     PC_dsnum  = 0    ; /* number of them */
00026 static int     PC_brnum  = 0    ; /* number of bricks */
00027 static int     PC_1ddum  = 0    ; /* number of dummies for 1D files */
00028 
00029 #define PC_lprin_calc PC_brnum
00030 
00031 static THD_3dim_dataset ** PC_dset = NULL ; /* pointers to datasets */
00032 
00033 static char PC_prefix[THD_MAX_PREFIX] = "pc" ;
00034 
00035 static float ** PC_brickdata ;   /* pointer to data bricks */
00036 
00037 static byte * PC_mask      = NULL ;   /* 15 Sep 1999 */
00038 static int    PC_mask_nvox = 0 ;
00039 static int    PC_mask_hits = 0 ;
00040 
00041 /*--------------------------- useful macros ------------------------*/
00042 
00043    /** i'th element of j'th input brick **/
00044 
00045 #define XX(i,j) PC_brickdata[(j)][(i)]
00046 
00047    /** i,j element of covariance matrix **/
00048 
00049 #define AA(i,j) (aa[(i)+(j)*adim])
00050 
00051    /** i'th element of j'th eigenvector **/
00052 
00053 #define VV(i,j) (zout[(i)+(j)*adim])
00054 
00055 /*--------------------------- prototypes ---------------------------*/
00056 
00057 void PC_read_opts( int , char ** ) ;
00058 void PC_syntax(char *) ;
00059 
00060 #undef USE_LAPACK
00061 #ifdef USE_LAPACK
00062 
00063    /** Eigenvalue routine from LAPACK **/
00064 
00065    extern int dsyevx_( char * , char * , char * , int * ,
00066                        double * , int * , double * , double * , int * ,
00067                        int * , double * , int * , double * ,
00068                        double * , int * , double * , int * ,
00069                        int * , int * , int * ) ;
00070 #else
00071 
00072    /** Use EISPACK instead */
00073 
00074 #  include "cs.h"
00075 #endif
00076 
00077 /*--------------------------------------------------------------------
00078    read the arguments, and load the global variables
00079 ----------------------------------------------------------------------*/
00080 
00081 void PC_read_opts( int argc , char * argv[] )
00082 {
00083    int nopt = 1 ;
00084    float val ;
00085    int  kk, nx, ny, nz, nxyz, mm,nn ;
00086    THD_3dim_dataset * dset ;
00087 
00088    while( nopt < argc && argv[nopt][0] == '-' ){
00089 
00090       /**** -verbose ****/
00091 
00092       if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00093          PC_be_quiet = 0 ;
00094          nopt++ ; continue ;
00095       }
00096 
00097       /**** -float ****/
00098 
00099       if( strncmp(argv[nopt],"-float",6) == 0 ){
00100          PC_do_float = 1 ;
00101          nopt++ ; continue ;
00102       }
00103 
00104       /**** -prefix prefix ****/
00105 
00106       if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00107          nopt++ ;
00108          if( nopt >= argc ) PC_syntax("need argument after -prefix!") ;
00109          MCW_strncpy( PC_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00110          continue ;
00111       }
00112 
00113       /**** -dmean ****/
00114 
00115       if( strncmp(argv[nopt],"-dmean",6) == 0 ){
00116          PC_dmean = 1 ;
00117          nopt++ ; continue ;
00118       }
00119 
00120       /**** -vmean ****/
00121 
00122       if( strncmp(argv[nopt],"-vmean",6) == 0 ){
00123          PC_vmean = 1 ;
00124          nopt++ ; continue ;
00125       }
00126 
00127       /**** -vnorm ****/
00128 
00129       if( strncmp(argv[nopt],"-vnorm",6) == 0 ){
00130          PC_vnorm = 1 ;
00131          nopt++ ; continue ;
00132       }
00133 
00134       /**** -normalize ****/
00135 
00136       if( strncmp(argv[nopt],"-normalize",6) == 0 ){
00137          PC_normalize = 1 ;
00138          nopt++ ; continue ;
00139       }
00140 
00141       /**** -pcsave # ****/
00142 
00143       if( strncmp(argv[nopt],"-pcsave",6) == 0 ){
00144          nopt++ ;
00145          if( nopt >= argc ) PC_syntax("need argument after -pcsave!") ;
00146          PC_lprin_save = strtol( argv[nopt] , NULL , 10 ) ;
00147          if( PC_lprin_save <  0 ) PC_syntax("value after -pcsave is illegal!") ;
00148          nopt++ ; continue ;
00149       }
00150 
00151       /**** -1ddum # ****/
00152 
00153       if( strncmp(argv[nopt],"-1ddum",6) == 0 ){
00154          nopt++ ;
00155          if( nopt >= argc ) PC_syntax("need argument after -1ddum!") ;
00156          PC_1ddum = strtol( argv[nopt] , NULL , 10 ) ;
00157          if( PC_1ddum < 0 ) PC_syntax("value after -1ddum is illegal!") ;
00158          nopt++ ; continue ;
00159       }
00160 
00161       /**** -mask mset [15 Sep 1999] ****/
00162 
00163       if( strncmp(argv[nopt],"-mask",5) == 0 ){
00164          THD_3dim_dataset * mset ; int ii,mc ;
00165          nopt++ ;
00166          if( nopt >= argc ) PC_syntax("need argument after -mask!") ;
00167          mset = THD_open_dataset( argv[nopt] ) ;
00168          if( mset == NULL ) PC_syntax("can't open -mask dataset!") ;
00169          PC_mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00170          PC_mask_nvox = DSET_NVOX(mset) ;
00171          DSET_delete(mset) ;
00172          if( PC_mask == NULL ) PC_syntax("can't use -mask dataset!") ;
00173          for( ii=mc=0 ; ii < PC_mask_nvox ; ii++ ) if( PC_mask[ii] ) mc++ ;
00174          if( mc == 0 ) PC_syntax("mask is all zeros!") ;
00175          printf("--- %d voxels in mask\n",mc) ;
00176          PC_mask_hits = mc ;
00177          nopt++ ; continue ;
00178       }
00179 
00180       /**** unknown switch ****/
00181 
00182       fprintf(stderr,"\n*** unrecognized option %s\n",argv[nopt]) ;
00183       exit(1) ;
00184 
00185    }  /* end of loop over options */
00186 
00187    /*--- a simple consistency check ---*/
00188 
00189    if( PC_vmean && PC_dmean )
00190       PC_syntax("can't have both -dmean and -vmean!") ;
00191 
00192    /*--- rest of inputs are dataset names ---*/
00193 
00194    PC_dsnum  = argc - nopt ;
00195    if( PC_dsnum < 1 ) PC_syntax("no input dataset names?") ;
00196 
00197    PC_dsname = (char **) malloc( sizeof(char *) * PC_dsnum ) ;
00198    for( kk=0 ; kk < PC_dsnum ; kk++ ) PC_dsname[kk] = argv[kk+nopt] ;
00199 
00200    PC_dset = (THD_3dim_dataset **) malloc( sizeof(THD_3dim_dataset *) * PC_dsnum ) ;
00201    for( kk=0 ; kk < PC_dsnum ; kk++ ){
00202       PC_dset[kk] = dset = THD_open_dataset( PC_dsname[kk] ) ;  /* allow for selector */
00203       if( !ISVALID_3DIM_DATASET(dset) ){
00204          fprintf(stderr,"\n*** can't open dataset file %s\n",PC_dsname[kk]) ;
00205          exit(1) ;
00206       }
00207       PC_brnum += DSET_NVALS(dset) ;
00208    }
00209    if( PC_brnum < 2 ) PC_syntax("need at least 2 input bricks!") ;
00210 
00211    /*--- another consistency check ---*/
00212 
00213    if( PC_lprin_save <= 0 )
00214       PC_lprin_save = PC_brnum ;
00215    else if( PC_lprin_save > PC_brnum )
00216       PC_syntax("can't save more components than input bricks!") ;
00217 
00218    /*--- load bricks for all input datasets ---*/
00219 
00220    nx = DSET_NX(PC_dset[0]) ;
00221    ny = DSET_NY(PC_dset[0]) ;
00222    nz = DSET_NZ(PC_dset[0]) ;
00223    nxyz = nx * ny * nz ;      /* Total number of voxels per brick */
00224 
00225    /* 15 Sep 1999: check if mask is right size */
00226 
00227    if( PC_mask_nvox > 0 && PC_mask_nvox != nxyz )
00228       PC_syntax("mask and input dataset bricks don't match in size!") ;
00229 
00230    PC_brickdata = (float **) malloc( sizeof(float *) * PC_brnum ) ;
00231 
00232    nn = 0 ; /* current brick index */
00233 
00234    if( !PC_be_quiet ){ printf("--- read dataset bricks"); fflush(stdout); }
00235 
00236    for( kk=0 ; kk < PC_dsnum ; kk++ ){
00237 
00238       if( DSET_NVOX(PC_dset[kk]) != nxyz ) {
00239          fprintf(stderr,
00240                  "\n*** dataset in file %s nonconformant with first dataset!\n"
00241                  "    nx1=%d ny1=%d nz1=%d nx=%d ny=%d nz=%d\n",
00242                  PC_dsname[kk], nx, ny, nz ,
00243                  DSET_NX(PC_dset[kk]), DSET_NY(PC_dset[kk]), DSET_NZ(PC_dset[kk]) ) ;
00244          exit(1) ;
00245       }
00246 
00247       DSET_load( PC_dset[kk] ) ;
00248       if( !DSET_LOADED( PC_dset[kk] ) ){
00249          fprintf(stderr,"\n*** Can't load dataset %s BRIK from disk!\n",PC_dsname[kk]) ;
00250          exit(1) ;
00251       }
00252       if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00253 
00254       /* copy brick data into float storage */
00255 
00256       for( mm=0 ; mm < DSET_NVALS(PC_dset[kk]) ; mm++,nn++ ){
00257 
00258          PC_brickdata[nn] = (float *) malloc( sizeof(float) * nxyz ) ;
00259 
00260          if( PC_brickdata[nn] == NULL )
00261             PC_syntax("*** can't malloc intermediate storage") ;
00262 
00263          EDIT_coerce_type( nxyz , DSET_BRICK_TYPE(PC_dset[kk],mm) ,
00264                                   DSET_ARRAY(PC_dset[kk],mm) ,
00265                            MRI_float , PC_brickdata[nn] ) ;
00266 
00267          DSET_unload_one( PC_dset[kk] , mm ) ;
00268 
00269          if( PC_mask != NULL ){             /* 15 Sep 1999 */
00270             int kk ;
00271             for( kk=0 ; kk < nxyz ; kk++ )
00272                if( !PC_mask[kk] ) PC_brickdata[nn][kk] = 0.0 ;
00273          }
00274 
00275          if( !PC_be_quiet ){ printf("."); fflush(stdout); }
00276       }
00277 
00278       if( kk == 0 ){
00279          DSET_unload( PC_dset[kk] ) ;  /* don't need dataset's data anymore */
00280       } else {
00281          DSET_delete( PC_dset[kk] ) ;  /* don't need this at all anymore */
00282          PC_dset[kk] = NULL ;
00283       }
00284 
00285    }
00286    if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00287 
00288    free( PC_dsname ) ; PC_dsname = NULL ;
00289    return ;
00290 }
00291 
00292 /*---------------------------------------------------------------------------*/
00293 
00294 void PC_syntax(char * msg)
00295 {
00296    if( msg != NULL ){ fprintf(stderr,"\n*** %s\n",msg) ; exit(1) ; }
00297 
00298    printf(
00299     "Principal Component Analysis of 3D Datasets\n"
00300     "Usage: 3dpc [options] dataset dataset ...\n"
00301     "\n"
00302     "Each input dataset may have a sub-brick selector list.\n"
00303     "Otherwise, all sub-bricks from a dataset will be used.\n"
00304     "\n"
00305     "OPTIONS:\n"
00306     "  -dmean        = remove the mean from each input brick (across space)\n"
00307     "  -vmean        = remove the mean from each input voxel (across bricks)\n"
00308     "                    [N.B.: -dmean and -vmean are mutually exclusive]\n"
00309     "                    [default: don't remove either mean]\n"
00310     "  -vnorm        = L2 normalize each input voxel time series\n"
00311     "                    [occurs after the de-mean operations above,]\n"
00312     "                    [and before the brick normalization below. ]\n"
00313     "  -normalize    = L2 normalize each input brick (after mean subtraction)\n"
00314     "                    [default: don't normalize]\n"
00315     "  -pcsave sss   = 'sss' is the number of components to save in the output;\n"
00316     "                    it can't be more than the number of input bricks\n"
00317     "                    [default = all of them = number of input bricks]\n"
00318     "  -prefix pname = Name for output dataset (will be a bucket type);\n"
00319     "                    also, the eigen-timeseries will be in 'pname'.1D\n"
00320     "                    (all of them) and in 'pnameNN.1D' for eigenvalue\n"
00321     "                    #NN individually (NN=00 .. 'sss'-1, corresponding\n"
00322     "                    to the brick index in the output dataset)\n"
00323     "                    [default value of pname = 'pc']\n"
00324     "  -1ddum ddd    = Add 'ddd' dummy lines to the top of each *.1D file.\n"
00325     "                    These lines will have the value 999999, and can\n"
00326     "                    be used to align the files appropriately.\n"
00327     "                    [default value of ddd = 0]\n"
00328     "  -verbose      = Print progress reports during the computations\n"
00329     "  -float        = Save eigen-bricks as floats\n"
00330     "                    [default = shorts, scaled so that |max|=10000]\n"
00331     "  -mask mset    = Use the 0 sub-brick of dataset 'mset' as a mask\n"
00332     "                    to indicate which voxels to analyze (a sub-brick\n"
00333     "                    selector is allowed) [default = use all voxels]\n"
00334    ) ;
00335 
00336    printf("\n" MASTER_SHORTHELP_STRING ) ;
00337 
00338    exit(0) ;
00339 }
00340 
00341 /*------------------------------------------------------------------*/
00342 
00343 int main( int argc , char * argv[] )
00344 {
00345    int nx,ny,nz , nxyz , ii,jj,ll , nn,mm,mmmm , npos,nneg ;
00346    float fmax , ftem ;
00347    THD_3dim_dataset * dset , * new_dset ;
00348    double * dsdev = NULL ;
00349    float  * fout , * perc ;
00350    short  * bout  ;
00351    register float  sum ;
00352    register double dsum ;
00353    register int    kk ;
00354    int ifirst,ilast,idel ;
00355    int ierror;          /* number of errors in editing data */
00356    MRI_IMAGE * vecim ;
00357    char vname[THD_MAX_NAME] ;
00358 
00359    /** data for eigenvalue routine (some only used with LAPACK) **/
00360 
00361    double * aa , * wout , * zout , * work ;
00362    double abstol , atrace ;
00363    int adim , il , iu , mout , lwork , info ;
00364    int * iwork , * ifail ;
00365 
00366    /*-- read command line arguments --*/
00367 
00368    if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) PC_syntax(NULL) ;
00369 
00370    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
00371 
00372    mainENTRY("3dpc main"); machdep(); PRINT_VERSION("3dpc") ;
00373 
00374    { int new_argc ; char ** new_argv ;
00375      addto_args( argc , argv , &new_argc , &new_argv ) ;
00376      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00377    }
00378 
00379    AFNI_logger("3dpc",argc,argv) ;
00380 
00381    PC_read_opts( argc , argv ) ;
00382 
00383    /*-- get dimensions --*/
00384 
00385    dset = PC_dset[0]    ;
00386    nx   = DSET_NX(dset) ;
00387    ny   = DSET_NY(dset) ;
00388    nz   = DSET_NZ(dset) ;
00389    nxyz = nx * ny * nz  ;
00390 
00391    nn = nxyz ;           /* vector length */
00392    mm = PC_brnum ;       /* number of vectors */
00393 
00394    /*-- space for eigenvalue computations --*/
00395 
00396    adim   = mm ;
00397    aa     = (double *) malloc( sizeof(double) * adim * adim ) ; /* matrix */
00398    wout   = (double *) malloc( sizeof(double) * adim ) ;        /* evals  */
00399 
00400    if( aa == NULL || wout == NULL )
00401       PC_syntax("can't malloc space for covariance matrix") ;
00402 
00403 #ifdef USE_LAPACK
00404    il     = adim + 1 - PC_lprin_calc ;  /* lowest index */
00405    iu     = adim ;                      /* upper index */
00406    abstol = 0.0 ;
00407    zout   = (double *) malloc( sizeof(double) * adim * PC_lprin_calc ) ;
00408    lwork  = 32 * adim ;
00409    work   = (double *) malloc( sizeof(double) * lwork ) ;
00410    iwork  = (int *)    malloc( sizeof(int) * 6 * adim ) ;
00411    ifail  = (int *)    malloc( sizeof(int) * adim ) ;
00412 
00413    if( zout == NULL || work == NULL || iwork == NULL || ifail == NULL )
00414       PC_syntax("can't malloc eigen workspace") ;
00415 #endif
00416 
00417    /*-- remove means, if ordered --*/
00418 
00419    if( PC_vmean ){
00420       float * vxmean = (float *) malloc( sizeof(float) * nn ) ;  /* voxel means */
00421 
00422       if( !PC_be_quiet ){ printf("--- remove timeseries means"); fflush(stdout); }
00423 
00424       for( kk=0 ; kk < nn ; kk++ ) vxmean[kk] = 0.0 ;
00425 
00426       for( jj=0 ; jj < mm ; jj++ ){
00427          for( kk=0 ; kk < nn ; kk++ ) vxmean[kk] += XX(kk,jj) ;
00428          if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00429       }
00430 
00431       sum = 1.0 / mm ;
00432       for( kk=0 ; kk < nn ; kk++ ) vxmean[kk] *= sum ;
00433 
00434       for( jj=0 ; jj < mm ; jj++ ){
00435          for( kk=0 ; kk < nn ; kk++ ) XX(kk,jj) -= vxmean[kk] ;
00436          if( !PC_be_quiet ){ printf("-"); fflush(stdout); }
00437       }
00438 
00439       free(vxmean) ;
00440       if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00441 
00442    } else if( PC_dmean ){
00443       if( !PC_be_quiet ){ printf("--- remove brick means"); fflush(stdout); }
00444 
00445       if( PC_mask == NULL ){
00446          for( jj=0 ; jj < mm ; jj++ ){
00447             sum = 0.0 ;
00448             for( kk=0 ; kk < nn ; kk++ ) sum += XX(kk,jj) ;
00449             if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00450             sum /= nn ;
00451             for( kk=0 ; kk < nn ; kk++ ) XX(kk,jj) -= sum ;
00452             if( !PC_be_quiet ){ printf("-"); fflush(stdout); }
00453          }
00454       } else {                           /* 15 Sep 1999 */
00455          for( jj=0 ; jj < mm ; jj++ ){
00456             sum = 0.0 ;
00457             for( kk=0 ; kk < nn ; kk++ ) if( PC_mask[kk] ) sum += XX(kk,jj) ;
00458             if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00459             sum /= PC_mask_hits ;
00460             for( kk=0 ; kk < nn ; kk++ ) if( PC_mask[kk] ) XX(kk,jj) -= sum ;
00461             if( !PC_be_quiet ){ printf("-"); fflush(stdout); }
00462          }
00463       }
00464       if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00465    }
00466 
00467    /*-- 07 July 1999: vnorm --*/
00468 
00469    if( PC_vnorm ){
00470       float * vxnorm = (float *) malloc( sizeof(float) * nn ) ;  /* voxel norms */
00471 
00472       if( !PC_be_quiet ){ printf("--- normalize timeseries"); fflush(stdout); }
00473 
00474       for( kk=0 ; kk < nn ; kk++ ) vxnorm[kk] = 0.0 ;
00475 
00476       for( jj=0 ; jj < mm ; jj++ ){
00477          for( kk=0 ; kk < nn ; kk++ ) vxnorm[kk] += XX(kk,jj) * XX(kk,jj) ;
00478          if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00479       }
00480 
00481       for( kk=0 ; kk < nn ; kk++ )
00482          if( vxnorm[kk] > 0.0 ) vxnorm[kk] = 1.0 / sqrt(vxnorm[kk]) ;
00483 
00484       for( jj=0 ; jj < mm ; jj++ ){
00485          for( kk=0 ; kk < nn ; kk++ ) XX(kk,jj) *= vxnorm[kk] ;
00486          if( !PC_be_quiet ){ printf("*"); fflush(stdout); }
00487       }
00488 
00489       free(vxnorm) ;
00490       if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00491    }
00492 
00493    /*-- load covariance matrix
00494         (very short code that takes a long time to run!) --*/
00495 
00496    if( !PC_be_quiet ){ printf("--- compute covariance matrix"); fflush(stdout); }
00497 
00498    idel = 1 ;                           /* ii goes forward */
00499    for( jj=0 ; jj < mm ; jj++ ){
00500 
00501       ifirst = (idel==1) ?    0 : jj ;  /* back and forth in ii to   */
00502       ilast  = (idel==1) ? jj+1 : -1 ;  /* maximize use of cache/RAM */
00503 
00504       for( ii=ifirst ; ii != ilast ; ii += idel ){
00505          dsum = 0.0 ;
00506          if( PC_mask == NULL ){
00507             for( kk=0 ; kk < nn ; kk++ ) dsum += XX(kk,ii) * XX(kk,jj)   ;
00508          } else {
00509             for( kk=0 ; kk < nn ; kk++ )
00510                        if( PC_mask[kk] ) dsum += XX(kk,ii) * XX(kk,jj)   ;
00511          }
00512          AA(ii,jj) = AA(jj,ii) = dsum ;
00513       }
00514 
00515       if( !PC_be_quiet ){ printf("+"); fflush(stdout); }
00516 
00517       idel = -idel ;                    /* reverse direction of ii */
00518    }
00519    if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00520 
00521    /*-- check diagonal for OK-ness --**/
00522 
00523    atrace = 0.0 ;
00524    ii     = 0 ;
00525    for( jj=0 ; jj < mm ; jj++ ){
00526       if( AA(jj,jj) <= 0.0 ){
00527          fprintf(stderr,"*** covariance diagonal (%d,%d) = %g\n",
00528                  jj+1,jj+1,AA(jj,jj)) ;
00529          ii++ ;
00530       }
00531       atrace += AA(jj,jj) ;
00532    }
00533    if( ii > 0 ){ printf("*** program exiting right here and now!\n"); exit(1); }
00534 
00535    if( !PC_be_quiet ){ printf("--- covariance trace = %g\n",atrace); fflush(stdout); }
00536 
00537    /*-- normalize, if desired --*/
00538 
00539    if( PC_normalize ){
00540       if( !PC_be_quiet ){ printf("--- normalizing covariance"); fflush(stdout); }
00541 
00542       dsdev = (double *) malloc( sizeof(double) * mm ) ; /* brick stdev */
00543 
00544       for( jj=0 ; jj < mm ; jj++ ) dsdev[jj] = sqrt( AA(jj,jj) ) ;
00545 
00546       for( jj=0 ; jj < mm ; jj++ )
00547          for( ii=0 ; ii < mm ; ii++ ) AA(ii,jj) /= (dsdev[ii]*dsdev[jj]) ;
00548 
00549       atrace = mm ;
00550 
00551       if( !PC_be_quiet ){ printf("\n"); fflush(stdout); }
00552    }
00553 
00554    if( !PC_be_quiet ){ printf("--- compute eigensolution\n"); fflush(stdout); }
00555 
00556 #ifdef USE_LAPACK
00557    (void) dsyevx_( "V"     , /* eigenvalues and vectors */
00558                    "I"     , /* a subrange of eigenvalues */
00559                    "U"     , /* use upper triangle of A */
00560                    &adim   , /* dimension of A */
00561                    aa      , /* the matrix A */
00562                    &adim   , /* leading dimension of A */
00563                    NULL    , /* not used */
00564                    NULL    , /* not used */
00565                    &il     , /* lowest eigen-index */
00566                    &iu     , /* highest eigen-index */
00567                    &abstol , /* tolerance */
00568                    &mout   , /* output # of eigenvalues */
00569                    wout    , /* output eigenvalues */
00570                    zout    , /* output eigenvectors */
00571                    &adim   , /* leading dimension of zout */
00572                    work    , /* double work array */
00573                    &lwork  , /* size of work array */
00574                    iwork   , /* another work array */
00575                    ifail   , /* output failure list */
00576                    &info     /* results code */
00577                  ) ;
00578 
00579    free(aa) ; free(work) ; free(iwork) ; free(ifail) ;
00580 
00581    if( info != 0 ){
00582       fprintf(stderr,"*** DSYEVX returns error code info=%d\n",info);
00583       if( info < 0 ) exit(1) ;
00584    }
00585 #else
00586    symeig_double( mm , aa , wout ) ;  /* eigenvectors go over aa  */
00587    zout = aa ;                        /* eigenvalues go into wout */
00588 #endif
00589 
00590    if( !PC_be_quiet ) printf("\n") ;
00591 
00592    sum = 0.0 ;
00593 
00594    perc = (float *) malloc( sizeof(float) * PC_lprin_calc ) ;
00595 
00596    printf("Num.  --Eigenvalue--  -Var.Fraction-  -Cumul.Fract.-\n") ;
00597    for( jj=0 ; jj < PC_lprin_calc ; jj++ ){
00598       ll = PC_lprin_calc - 1-jj ;      /* reversed order of eigensolution! */
00599       perc[jj] = wout[ll]/atrace ;
00600       sum     += perc[jj] ;
00601       printf("%4d  %14.7g  %14.7g  %14.7g\n",
00602              jj+1 , wout[ll] , perc[jj] , sum ) ;
00603    }
00604    fflush(stdout) ;
00605 
00606    /*--- form and save output dataset ---*/
00607 
00608    dset = PC_dset[0] ;
00609    new_dset = EDIT_empty_copy( dset ) ;
00610 
00611    if( PC_dsnum == 1 ) tross_Copy_History( dset , new_dset ) ;
00612    tross_Make_History( "3dpc" , argc,argv , new_dset ) ;
00613 
00614    EDIT_dset_items( new_dset,
00615                        ADN_prefix    , PC_prefix ,
00616                        ADN_nvals     , PC_lprin_save ,
00617                        ADN_ntt       , 0 ,
00618                        ADN_func_type , ISANAT(dset) ? ANAT_BUCK_TYPE
00619                                                     : FUNC_BUCK_TYPE ,
00620                     ADN_none ) ;
00621 
00622    if( THD_is_file(DSET_HEADNAME(new_dset)) ){
00623       fprintf(stderr,
00624               "\n*** Output dataset %s already exists--will be destroyed!\n",
00625               DSET_HEADNAME(new_dset) ) ;
00626 
00627    } else if( !PC_be_quiet ){
00628       printf("--- output dataset %s" , DSET_BRIKNAME(new_dset) ) ;
00629       fflush(stdout) ;
00630    }
00631 
00632    fout = (float *) malloc( sizeof(float) * nn ) ; /* output buffer */
00633 
00634    for( jj=0 ; jj < PC_lprin_save ; jj++ ){
00635       ll = PC_lprin_calc - 1-jj ;
00636 
00637       /** output = weighted sum of input datasets,
00638                    with weights from the ll'th eigenvector **/
00639 
00640       for( kk=0 ; kk < nn ; kk++ ) fout[kk] = 0.0 ;
00641 
00642       for( ii=0 ; ii < mm ; ii++ ){
00643          sum = VV(ii,ll) ; if( PC_normalize ) sum /= dsdev[ii] ;
00644 
00645          for( kk=0 ; kk < nn ; kk++ ) fout[kk] += XX(kk,ii) * sum ;
00646       }
00647 
00648       fmax = 0.0 ; npos = nneg = 0 ;
00649       for( kk=0 ; kk < nn ; kk++ ){
00650          ftem = fabs(fout[kk]) ; if( fmax < ftem ) fmax = ftem ;
00651               if( fout[kk] > 0 ) npos++ ;
00652          else if( fout[kk] < 0 ) nneg++ ;
00653       }
00654 
00655       if( PC_do_float ){
00656          if( nneg > npos )
00657             for( kk=0 ; kk < nn ; kk++ ) fout[kk] = -fout[kk] ;
00658 
00659          EDIT_substitute_brick( new_dset , jj , MRI_float , fout ) ;
00660 
00661          fout = (float *) malloc( sizeof(float) * nn ) ;  /* new buffer */
00662       } else {
00663 
00664          bout = (short *) malloc( sizeof(short) * nn ) ; /* output buffer */
00665          if( fmax != 0.0 ){
00666             fmax = 10000.49/fmax ; if( nneg > npos ) fmax = -fmax ;
00667             for( kk=0 ; kk < nn ; kk++ ) bout[kk] = fmax * fout[kk] ;
00668          } else {
00669             for( kk=0 ; kk < nn ; kk++ ) bout[kk] = 0.0 ;
00670          }
00671          EDIT_substitute_brick( new_dset , jj , MRI_short , bout ) ;
00672       }
00673 
00674       sprintf(vname,"var=%6.3f%%" , 100.0*perc[jj]+0.499 ) ;
00675       EDIT_BRICK_LABEL( new_dset , jj , vname ) ;
00676 
00677       if( !PC_be_quiet ){ printf(".") ; fflush(stdout); }
00678    }
00679    free(fout) ;
00680 
00681    DSET_write(new_dset) ;
00682    fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
00683    DSET_delete(new_dset) ;
00684    if( !PC_be_quiet ){ printf("!\n") ; fflush(stdout); }
00685 
00686    /*-- write eigenvectors also --*/
00687 
00688    mmmm  = mm + PC_1ddum ;
00689    vecim = mri_new( PC_lprin_save , mmmm , MRI_float ) ;
00690    fout  = MRI_FLOAT_PTR(vecim) ;
00691    for( jj=0 ; jj < PC_lprin_save ; jj++ ){
00692       ll = PC_lprin_calc - 1-jj ;
00693       for( ii=0 ; ii < PC_1ddum ; ii++ )
00694          fout[jj + ii*PC_lprin_save] = 999999.0 ;
00695       for( ii=0 ; ii < mm ; ii++ )
00696          fout[jj + (ii+PC_1ddum)*PC_lprin_save] = VV(ii,ll) ;
00697    }
00698    sprintf(vname,"%s.1D",PC_prefix) ;
00699    mri_write_ascii( vname, vecim ) ;
00700    mri_free(vecim) ;
00701 
00702    for( jj=0 ; jj < PC_lprin_save ; jj++ ){
00703       ll = PC_lprin_calc - 1-jj ;
00704       vecim = mri_new( 1 , mmmm , MRI_float ) ;
00705       fout  = MRI_FLOAT_PTR(vecim) ;
00706       for( ii=0 ; ii < PC_1ddum ; ii++ ) fout[ii] = 999999.0 ;
00707       for( ii=0 ; ii < mm ; ii++ ) fout[ii+PC_1ddum] = VV(ii,ll) ;
00708       sprintf(vname,"%s%02d.1D",PC_prefix,jj) ;
00709       mri_write_ascii( vname, vecim ) ;
00710       mri_free(vecim) ;
00711    }
00712 
00713 #if 0
00714    free(PC_dset) ;
00715    for( ii=0 ; ii < mm ; ii++ ) free(PC_brickdata[ii]) ;
00716    free(PC_brickdata) ;
00717    free(aa) ; free(wout) ; free(perc) ; if( dsdev!=NULL ) free(dsdev) ;
00718 #endif
00719 
00720    exit(0) ;
00721 }