Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

analyzer.cpp

Go to the documentation of this file.
00001 
00006 /*
00007  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00008  * Copyright (c) 2000-2006 Baylor College of Medicine
00009  *
00010  * This software is issued under a joint BSD/GNU license. You may use the
00011  * source code in this file under either license. However, note that the
00012  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00013  * so you are responsible for compliance with the licenses of these packages
00014  * if you opt to use BSD licensing. The warranty disclaimer below holds
00015  * in either instance.
00016  *
00017  * This complete copyright notice must be included in any revised version of the
00018  * source code. Additional authorship citations may be added, but existing
00019  * author citations must be preserved.
00020  *
00021  * This program is free software; you can redistribute it and/or modify
00022  * it under the terms of the GNU General Public License as published by
00023  * the Free Software Foundation; either version 2 of the License, or
00024  * (at your option) any later version.
00025  *
00026  * This program is distributed in the hope that it will be useful,
00027  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00028  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00029  * GNU General Public License for more details.
00030  *
00031  * You should have received a copy of the GNU General Public License
00032  * along with this program; if not, write to the Free Software
00033  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00034  *
00035  * */
00036 
00037 #include <ctime>
00038 #include <memory>
00039 #include "emdata.h"
00040 #include "analyzer.h"
00041 #include "sparx/analyzer_sparx.h"
00042 #include "util.h"
00043 #include "cmp.h"
00044 #include "sparx/lapackblas.h"
00045 #include "sparx/varimax.h"
00046 
00047 using namespace EMAN;
00048 
00049 namespace EMAN {
00050 
00051         const string PCAsmall::NAME = "pca";
00052         const string PCAlarge::NAME = "pca_large";
00053         const string varimax::NAME = "varimax";
00054         const string KMeansAnalyzer::NAME = "kmeans";
00055         const string SVDAnalyzer::NAME = "svd_gsl";
00056 
00057         template <> Factory < Analyzer >::Factory()
00058         {
00059                 force_add<PCAsmall>();
00060                 force_add<PCAlarge>();
00061                 force_add<varimax>();
00062                 force_add<KMeansAnalyzer>();
00063                 force_add<SVDAnalyzer>();
00064         }
00065 
00066 }
00067 
00068 int Analyzer::insert_images_list(vector<EMData *> image_list)
00069 {
00070         vector<EMData *>::const_iterator iter;
00071                 for(iter=image_list.begin(); iter!=image_list.end(); ++iter) {
00072                         insert_image(*iter);
00073                 }
00074         return 0;
00075 }
00076 
00077 void KMeansAnalyzer::set_params(const Dict & new_params)
00078 {
00079         params = new_params;
00080         if (params.has_key("ncls")) ncls = params["ncls"];
00081         if (params.has_key("maxiter"))maxiter = params["maxiter"];
00082         if (params.has_key("minchange"))minchange = params["minchange"];
00083         if (params.has_key("mininclass"))mininclass = params["mininclass"];
00084         if (params.has_key("slowseed"))slowseed = params["slowseed"];
00085         if (params.has_key("verbose"))verbose = params["verbose"];
00086         if (params.has_key("calcsigmamean")) calcsigmamean=params["calcsigmamean"];
00087 
00088 }
00089 
00090 vector<EMData *> KMeansAnalyzer::analyze()
00091 {
00092 if (ncls<=1) return vector<EMData *>();
00093 //srandom(time(0));
00094 
00095 // These are the class centers, start each with a random image
00096 int nptcl=images.size();
00097 int nclstot=ncls;
00098 if (calcsigmamean) centers.resize(nclstot*2);
00099 else centers.resize(nclstot);
00100 if (mininclass<1) mininclass=1;
00101 
00102 for (int i=0; i<nptcl; i++) images[i]->set_attr("is_ok_center",(int)5);  // if an image becomes part of too small a set, it will (eventually) be marked as a bad center
00103 
00104 if (slowseed) {
00105         if (ncls>25) slowseed=ncls/25+1;        // this becomes the number to seed in each step
00106 //      if (maxiter<ncls*3+20) maxiter=ncls*3+20;       // We need to make sure we have enough iterations to seed all of the classes
00107 //      ncls=2;
00108 }
00109 
00110 for (int i=0; i<ncls; i++) {
00111         // Fixed by d.woolford, Util.get_irand is inclusive (added a -1)
00112         centers[i]=images[Util::get_irand(0,nptcl-1)]->copy();
00113 
00114 }
00115 
00116 if (calcsigmamean) {
00117         for (int i=nclstot; i<nclstot*2; i++) centers[i]=new EMData(images[0]->get_xsize(),images[0]->get_ysize(),images[0]->get_zsize());
00118 }
00119 
00120 
00121 for (int i=0; i<maxiter; i++) {
00122         nchanged=0;
00123         reclassify();
00124         if (verbose) printf("iter %d>  %d (%d)\n",i,nchanged,ncls);
00125         if (nchanged<minchange && ncls==nclstot) break;
00126         update_centers();
00127 
00128         if (slowseed && i%3==2 && ncls<nclstot) {
00129                 for (int j=0; j<slowseed && ncls<nclstot; j++) {
00130                         centers[ncls]=0;
00131                         ncls++;
00132                 }
00133                 reseed();
00134         }
00135 }
00136 update_centers(calcsigmamean);
00137 
00138 return centers;
00139 }
00140 
00141 void KMeansAnalyzer::update_centers(int sigmas) {
00142 int nptcl=images.size();
00143 //int repr[ncls];
00144 int * repr = new int[ncls];
00145 
00146 for (int i=0; i<ncls; i++) {
00147         centers[i]->to_zero();
00148         if (sigmas) centers[i+ncls]->to_zero();
00149         repr[i]=0;
00150 }
00151 
00152 // compute new position for each center
00153 for (int i=0; i<nptcl; i++) {
00154         int cid=images[i]->get_attr("class_id");
00155         if ((int)images[i]->get_attr("is_ok_center")>0) {
00156                 centers[cid]->add(*images[i]);
00157                 if (sigmas) centers[cid+ncls]->addsquare(*images[i]);
00158                 repr[cid]++;
00159         }
00160 }
00161 
00162 for (int i=0; i<ncls; i++) {
00163         // If this class is too small
00164         if (repr[i]<mininclass) {
00165                 // find all of the particles in the class, and decrement their "is_ok_center" counter.
00166                 // when it reaches zero the particle will no longer participate in determining the location of a center
00167                 for (int j=0; j<nptcl; j++) {
00168                         if ((int)images[j]->get_attr("class_id")==i) images[i]->set_attr("is_ok_center",(int)images[i]->get_attr("is_ok_center")-1);
00169                 }
00170                 // Mark the center for reseeding
00171                 delete centers[i];
00172                 centers[i]=0;
00173                 repr[i]=0;
00174         }
00175         // finishes off the statistics we started computing above
00176         else {
00177                 centers[i]->mult((float)1.0/(float)(repr[i]));
00178                 centers[i]->set_attr("ptcl_repr",repr[i]);
00179                 if (sigmas) {
00180                         centers[i+ncls]->mult((float)1.0/(float)(repr[i]));             // sum of squares over n
00181                         centers[i+ncls]->subsquare(*centers[i]);                                        // subtract the mean value squared
00182                         centers[i+ncls]->process("math.sqrt");                                  // square root
00183                         centers[i+ncls]->mult((float)1.0/(float)sqrt((float)repr[i]));          // divide by sqrt(N) to get std. dev. of mean
00184                 }
00185 
00186         }
00187         if (verbose>1) printf("%d(%d)\t",i,(int)repr[i]);
00188 }
00189 
00190 if (verbose>1) printf("\n");
00191 
00192 reseed();
00193 
00194 delete [] repr;
00195 }
00196 
00197 // This will look for any unassigned points and reseed each inside the class with the broadest distribution widely distributed
00198 void KMeansAnalyzer::reseed() {
00199 int nptcl=images.size();
00200 int i,j;
00201 
00202 // if no classes need reseeding just return
00203 for (i=0; i<ncls; i++) {
00204         if (!centers[i]) break;
00205 }
00206 if (i==ncls) return;
00207 
00208 // make a list of all particles which could be centers
00209 vector<int> goodcen;
00210 for (int i=0; i<nptcl; i++) if ((int)images[i]->get_attr("is_ok_center")>0) goodcen.push_back(i);
00211 
00212 if (goodcen.size()==0) throw UnexpectedBehaviorException("Kmeans ran out of valid center particles with the provided parameters");
00213 
00214 // pick a random particle for the new seed
00215 for (i=0; i<ncls; i++) {
00216         if (centers[i]) continue;               // center doesn't need reseeding
00217         j=Util::get_irand(0,goodcen.size()-1);
00218         centers[i]=images[j]->copy();
00219         centers[i]->set_attr("ptcl_repr",1);
00220         printf("reseed %d -> %d\n",i,j);
00221 }
00222 
00223 
00224 }
00225 
00226 
00227 void KMeansAnalyzer::reclassify() {
00228 int nptcl=images.size();
00229 
00230 Cmp *c = Factory < Cmp >::get("sqeuclidean");
00231 for (int i=0; i<nptcl; i++) {
00232         float best=1.0e38f;
00233         int bestn=0;
00234         for (int j=0; j<ncls; j++) {
00235                 float d=c->cmp(images[i],centers[j]);
00236 //images[i]->cmp("sqeuclidean",centers[j]);
00237                 if (d<best) { best=d; bestn=j; }
00238         }
00239         int oldn=images[i]->get_attr_default("class_id",0);
00240         if (oldn!=bestn) nchanged++;
00241         images[i]->set_attr("class_id",bestn);
00242 }
00243 delete c;
00244 }
00245 
00246 #define covmat(i,j) covmat[ ((j)-1)*nx + (i)-1 ]
00247 #define imgdata(i)  imgdata[ (i)-1 ]
00248 int PCAsmall::insert_image(EMData * image)
00249 {
00250         if(mask==0)
00251                 throw NullPointerException("Null mask image pointer, set_params() first");
00252 
00253    EMData *maskedimage = Util::compress_image_mask(image,mask);
00254 
00255    int nx = maskedimage->get_xsize();
00256    float *imgdata = maskedimage->get_data();
00257    if (nx != ncov) {
00258       fprintf(stderr,"insert_image: something is wrong...\n");
00259       exit(1);
00260    }
00261 
00262    // there is a faster version of the following rank-1 update
00263    nimages++;
00264    for (int j = 1; j <= nx; j++)
00265        for (int i = 1; i<=nx; i++) {
00266            covmat(i,j) += imgdata(i)*imgdata(j);
00267    }
00268 
00269    EMDeletePtr(maskedimage);
00270    return 0;
00271 }
00272 #undef covmat
00273 
00274 #define eigvec(i,j) eigvec[(j)*ncov + (i)]
00275 vector<EMData*> PCAsmall::analyze()
00276 {
00277         float *eigvec;
00278         int status = 0;
00279 //              printf("start analyzing..., ncov = %d\n", ncov);
00280         eigval = (float*)calloc(ncov,sizeof(float));
00281         eigvec = (float*)calloc(ncov*ncov,sizeof(float));
00282         status = Util::coveig(ncov, covmat, eigval, eigvec);
00283 //       for (int i=1; i<=nvec; i++) printf("eigval = %11.4e\n",
00284 //            eigval[ncov-i]);
00285 
00286         // pack eigenvectors into the return imagelist
00287         EMData *eigenimage = new EMData();
00288         eigenimage->set_size(ncov,1,1);
00289         float *rdata = eigenimage->get_data();
00290         for (int j = 1; j<= nvec; j++) {
00291             for (int i = 0; i < ncov; i++) rdata[i] = eigvec(i,ncov-j);
00292 
00293                 EMData* recons_eigvec = Util::reconstitute_image_mask(eigenimage,mask);
00294                 recons_eigvec->set_attr( "eigval", eigval[j-1] );
00295             images.push_back(recons_eigvec);
00296         }
00297 
00298         free(eigvec);
00299         EMDeletePtr(eigenimage);
00300 
00301         return images;
00302 }
00303 #undef eigvec
00304 
00305 void PCAsmall::set_params(const Dict & new_params)
00306 {
00307         params = new_params;
00308         mask = params["mask"];
00309         nvec = params["nvec"];
00310 
00311         // count the number of pixels under the mask
00312         // (this is really ugly!!!)
00313         EMData *dummy = new EMData();
00314 
00315         int nx = mask->get_xsize();
00316         int ny = mask->get_ysize();
00317         int nz = mask->get_zsize();
00318 
00319         dummy->set_size(nx,ny,nz);
00320 
00321         EMData *dummy1d = Util::compress_image_mask(dummy,mask);
00322         ncov = dummy1d->get_xsize();
00323         EMDeletePtr(dummy);
00324         EMDeletePtr(dummy1d);
00325 
00326         // allocate and set up the covriance matrix
00327         nimages = 0;
00328         covmat = (float*)calloc(ncov*ncov,sizeof(float));
00329 }
00330 
00331 //------------------------------------------------------------------
00332 // for large-scale PCA incore
00333 
00334 int PCAlarge::insert_image(EMData * image)
00335 {
00336         if(mask==0)
00337                 throw NullPointerException("Null mask image pointer, set_params() first");
00338 
00339    EMData *maskedimage = Util::compress_image_mask(image,mask);
00340 
00341    FILE *fp;
00342    string scratchfile = string("maskedimages.scratch");
00343 
00344    fp = fopen(scratchfile.c_str(),"ab");
00345 
00346    int nx = maskedimage->get_xsize();
00347    float *imgdata = maskedimage->get_data();
00348    fwrite(imgdata, sizeof(float), nx, fp);
00349    nimages++;
00350 
00351    fclose(fp);
00352 
00353    EMDeletePtr(maskedimage);
00354 
00355    return 0;
00356 }
00357 
00358 void PCAlarge::set_params(const Dict & new_params)
00359 {
00360         params = new_params;
00361         mask = params["mask"];
00362         nvec = params["nvec"];
00363 
00364         // count the number of pixels under the mask
00365         // (this is really ugly!!!)
00366         EMData *dummy = new EMData();
00367 
00368         int nx = mask->get_xsize();
00369         int ny = mask->get_ysize();
00370         int nz = mask->get_zsize();
00371 
00372         dummy->set_size(nx,ny,nz);
00373 
00374         EMData *dummy1d = Util::compress_image_mask(dummy,mask);
00375 
00376         ncov = dummy1d->get_xsize();
00377 
00378         EMDeletePtr(dummy);
00379         EMDeletePtr(dummy1d);
00380         // no need to allocate the covariance matrix
00381         nimages = 0;
00382 }
00383 
00384 #define qmat(i,j)   qmat[((j)-1)*kstep + (i) -1]
00385 #define diag(i)     diag[(i)-1]
00386 #define rdata(i)    rdata[(i)-1]
00387 #define eigvec(i,j) eigvec[((j)-1)*ncov + (i)-1]
00388 #define eigval(i)   eigval[(i)-1]
00389 
00390 vector<EMData*> PCAlarge::analyze()
00391 {
00392         int status = 0;
00393         int ione = 1;
00394         float one = 1.0, zero = 0.0;
00395         char trans;
00396         float *eigvec;
00397         string scratchfile = string("maskedimages.scratch");
00398         char command[100];
00399 
00400         printf("start analyzing..., ncov = %d\n", ncov);
00401 
00402         float resnrm = 0.0;
00403 
00404         if ( nvec > nimages || nvec ==0 ) nvec = nimages;
00405         int nx = ncov;
00406 
00407         // the definition of kstep is purely a heuristic for right now
00408         int kstep = nvec + 20;
00409         if (kstep > nimages) kstep = nimages;
00410 
00411         float *diag    = new float[kstep];
00412         float *subdiag = new float[kstep-1];
00413         float *vmat    = new float[nx*kstep];
00414 
00415         // run kstep-step Lanczos factorization
00416         status = Lanczos(scratchfile, &kstep, diag, subdiag,
00417                          vmat, &resnrm);
00418 
00419         // remove scratch file
00420 #ifdef _WIN32
00421         if (_unlink(scratchfile.c_str()) == -1) {
00422                 fprintf(stderr,"PCAlarge: cannot remove scratchfile\n");
00423         }
00424 #else
00425         sprintf(command,"rm -f %s\n", scratchfile.c_str());
00426         status = system(command);
00427         if (status != 0) {
00428                 fprintf(stderr,"PCAlarge: cannot remove scratchfile\n");
00429         }
00430 #endif  //_WIN32
00431 
00432         char jobz[2] = "V";
00433         float *qmat  = new float[kstep*kstep];
00434         // workspace size will be optimized later
00435         int   lwork  = 100 + 4*kstep + kstep*kstep;
00436         int   liwork = 3+5*kstep;
00437 
00438         float *work  = new float[lwork];
00439         int   *iwork = new int[liwork];
00440         int   info = 0;
00441 
00442         // call LAPACK tridiagonal eigensolver
00443         sstevd_(jobz, &kstep, diag, subdiag, qmat, &kstep, work, &lwork,
00444                 iwork, &liwork, &info);
00445 
00446         // store eigenvalues
00447         eigval = (float*)calloc(ncov,sizeof(float));
00448         eigvec = (float*)calloc(ncov*nvec,sizeof(float));
00449 
00450         for (int j = 0; j < nvec; j++) {
00451             eigval[j] = diag(kstep-j);
00452         }
00453 
00454 //         for (int i=0; i<nvec; i++) printf("eigval = %11.4e\n",
00455 //             eigval[i]);
00456 
00457         // compute eigenvectors
00458         for (int j=1; j<=nvec; j++) {
00459             trans = 'N';
00460             sgemv_(&trans, &nx,  &kstep, &one, vmat, &nx, &qmat(1,kstep-j+1),
00461                    &ione, &zero, &eigvec(1,j), &ione);
00462         }
00463 
00464         // pack eigenvectors into the return imagelist
00465         EMData *eigenimage = new EMData();
00466         eigenimage->set_size(ncov,1,1);
00467         float *rdata = eigenimage->get_data();
00468         for (int j = 1; j<= nvec; j++) {
00469             for (int i = 1; i <= ncov; i++)
00470                 rdata(i) = eigvec(i,j);
00471 
00472             EMData* recons_eigvec = Util::reconstitute_image_mask(eigenimage,mask);
00473 
00474             recons_eigvec->set_attr( "eigval", eigval[j-1] );
00475 
00476             images.push_back( recons_eigvec );
00477         }
00478 
00479         free(eigvec);
00480         EMDeletePtr(eigenimage);
00481 
00482         return images;
00483 }
00484 #undef qmat
00485 #undef diag
00486 #undef rdata
00487 #undef eigvec
00488 #undef eigval
00489 
00490 #define TOL 1e-7
00491 #define V(i,j)      V[((j)-1)*imgsize + (i) - 1]
00492 #define v0(i)       v0[(i)-1]
00493 #define Av(i)       Av[(i)-1]
00494 #define subdiag(i)  subdiag[(i)-1]
00495 #define diag(i)     diag[(i)-1]
00496 #define hvec(i)     hvec[(i)-1]
00497 
00498 int PCAlarge::Lanczos(const string &maskedimages, int *kstep,
00499                       float  *diag, float *subdiag, float *V,
00500                       float  *beta)
00501 {
00502     /*
00503         Purpose: Compute a kstep-step Lanczos factorization
00504                  on the covariant matrix X*trans(X), where
00505                  X (imgstack) contains a set of images;
00506 
00507         Input:
00508            imgstack (vector <EMData*>) a set of images on which PCA is
00509                                        to be performed;
00510 
00511            kstep (int*) The maximum number of Lanczos iterations allowed.
00512                           If Lanczos terminates before kstep steps
00513                           is reached (an invariant subspace is found),
00514                           kstep returns the number of steps taken;
00515 
00516         Output:
00517            diag (float *) The projection of the covariant matrix into a
00518                           Krylov subspace of dimension at most kstep.
00519                           The projection is a tridiagonal matrix. The
00520                           diagonal elements of this matrix is stored in
00521                           the diag array.
00522 
00523            subdiag (float*) The subdiagonal elements of the projection
00524                             is stored here.
00525 
00526            V (float *)    an imgsize by kstep array that contains a
00527                           set of orthonormal Lanczos basis vectors;
00528 
00529            beta (float *) the residual norm of the factorization;
00530     */
00531     int i, iter;
00532 
00533     float alpha;
00534     int   ione = 1;
00535     float zero = 0.0, one = 1.0, mone = -1.0;
00536     int   status = 0;
00537 
00538     char trans;
00539     int  imgsize = 0;
00540     float *v0, *Av, *hvec, *htmp, *imgdata;
00541     FILE  *fp=NULL;
00542 
00543     if (nimages <= 0) {
00544         status = 2; // no image in the stack
00545         goto EXIT;
00546     }
00547 
00548     imgsize = ncov;
00549     if (nimages <= 0) {
00550         status = 3; // no image in the stack
00551         goto EXIT;
00552     }
00553 
00554     v0   = new float[imgsize];
00555     Av   = new float[imgsize];
00556     hvec = new float[*kstep];
00557     htmp = new float[*kstep];
00558     imgdata = new float[imgsize];
00559 
00560     if (v0 == NULL || Av == NULL || hvec == NULL ||
00561         htmp == NULL || imgdata == NULL) {
00562         fprintf(stderr, "Lanczos: failed to allocate v0,Av,hvec,htmp\n");
00563         status = -1;
00564         goto EXIT;
00565     }
00566 
00567     // may choose a random starting guess here
00568     for ( i = 1; i <= imgsize; i++)
00569     {
00570         v0(i) = 1.0;
00571         Av(i) = 0.0;
00572     }
00573 
00574     // normalize the starting vector
00575     *beta  = snrm2_(&imgsize, v0, &ione);
00576     for (i = 1; i<=imgsize; i++)
00577         V(i,1) = v0(i) / (*beta);
00578 
00579     // do Av <-- A*v0, where A is a cov matrix
00580     fp = fopen(maskedimages.c_str(),"rb");
00581     if (fp==NULL) {
00582         fprintf(stderr,"Lanczos: cannot open %s\n", maskedimages.c_str());
00583     }
00584     for (i = 0; i < nimages; i++) {
00585        fread(imgdata, sizeof(float), imgsize, fp);
00586        alpha = sdot_(&imgsize, imgdata, &ione, V, &ione);
00587        saxpy_(&imgsize, &alpha, imgdata, &ione, Av, &ione);
00588     }
00589     fclose(fp);
00590 
00591 
00592     // Av <--- Av - V(:,1)*V(:,1)'*Av
00593     diag(1) = sdot_(&imgsize, V, &ione, Av, &ione);
00594     alpha   = -diag(1);
00595     saxpy_(&imgsize, &alpha, V, &ione, Av, &ione);
00596 
00597     // main loop
00598     for ( iter = 2 ; iter <= *kstep ; iter++ ) {
00599         *beta = snrm2_(&imgsize, Av, &ione);
00600 
00601         if (*beta < TOL) {
00602             // found an invariant subspace, exit
00603             *kstep = iter;
00604             break;
00605         }
00606 
00607         subdiag(iter-1) = *beta;
00608         for ( i = 1 ; i <= imgsize ; i++ ) {
00609             V(i,iter) = Av(i) / (*beta);
00610         }
00611 
00612         // do Av <-- A*V(:,iter), where A is a cov matrix
00613         for (i = 0; i < imgsize; i++) Av[i] = 0;
00614         fp = fopen(maskedimages.c_str(),"rb");
00615         for (i = 0; i < nimages; i++) {
00616            fread(imgdata, sizeof(float), imgsize, fp);
00617            alpha = sdot_(&imgsize, imgdata, &ione, &V(1,iter), &ione);
00618            saxpy_(&imgsize, &alpha, imgdata, &ione, Av, &ione);
00619         }
00620         fclose(fp);
00621 
00622         // f <--- Av - V(:,1:iter)*V(:,1:iter)'*Av
00623         trans = 'T';
00624         status = sgemv_(&trans, &imgsize, &iter, &one, V, &imgsize, Av, &ione,
00625                         &zero , hvec    , &ione);
00626         trans = 'N';
00627         status = sgemv_(&trans, &imgsize, &iter, &mone, V, &imgsize, hvec,
00628                         &ione , &one    , Av, &ione);
00629 
00630         // one step of reorthogonalization
00631         trans = 'T';
00632         status = sgemv_(&trans, &imgsize, &iter, &one, V, &imgsize, Av, &ione,
00633                         &zero , htmp    , &ione);
00634         saxpy_(&iter, &one, htmp, &ione, hvec, &ione);
00635         trans = 'N';
00636         status = sgemv_(&trans, &imgsize, &iter, &mone, V, &imgsize, htmp,
00637                         &ione , &one    , Av, &ione);
00638         diag(iter) = hvec(iter);
00639     }
00640 
00641     EMDeleteArray(v0);
00642     EMDeleteArray(Av);
00643     EMDeleteArray(hvec);
00644     EMDeleteArray(htmp);
00645     EMDeleteArray(imgdata);
00646 
00647 EXIT:
00648     return status;
00649 
00650 }
00651 #undef v0
00652 #undef Av
00653 #undef V
00654 #undef hvec
00655 #undef diag
00656 #undef subdiag
00657 #undef TOL
00658 
00659 void varimax::set_params(const Dict & new_params)
00660 {
00661         params = new_params;
00662         m_mask = params["mask"];
00663 
00664         // count the number of pixels under the mask
00665         // (this is really ugly!!!)
00666         EMData *dummy = new EMData();
00667 
00668         int nx = m_mask->get_xsize();
00669         int ny = m_mask->get_ysize();
00670         int nz = m_mask->get_zsize();
00671 
00672         dummy->set_size(nx,ny,nz);
00673 
00674         EMData *dummy1d = Util::compress_image_mask(dummy,m_mask);
00675 
00676         m_nlen = dummy1d->get_xsize();
00677         m_nfac = 0;
00678 
00679         EMDeletePtr(dummy);
00680         EMDeletePtr(dummy1d);
00681 }
00682 
00683 int varimax::insert_image(EMData* image)
00684 {
00685         if(m_mask==0)
00686                 throw NullPointerException("Null mask image pointer, set_params() first");
00687 
00688     EMData* img1d = Util::compress_image_mask(image,m_mask);
00689 
00690     m_data.insert( m_data.end(), img1d->get_data(), img1d->get_data() + m_nlen );
00691 
00692     m_nfac++;
00693 
00694     Assert( (int)m_data.size() == m_nfac*m_nlen);
00695 
00696     return 0;
00697 }
00698 
00699 vector<EMData*> varimax::analyze()
00700 {
00701     int itmax = 10000;
00702     float eps = 1e-4f;
00703     int verbose = 1;
00704     float params[4];
00705     params[0] = 1.0;
00706     varmx( &m_data[0], m_nlen, m_nfac, IVARIMAX, params, NULL, itmax, eps, verbose);
00707 
00708     vector<EMData*> images;
00709 
00710     EMData* img1d = new EMData();
00711     img1d->set_size(m_nlen, 1, 1);
00712     for( int i=0; i < m_nfac; ++i )
00713     {
00714         float* imgdata = img1d->get_data();
00715 
00716         int offset = i * m_nlen;
00717         for( int i=0; i < m_nlen; ++i )
00718         {
00719             imgdata[i] = m_data[offset+i];
00720         }
00721 
00722         EMData* img = Util::reconstitute_image_mask(img1d,m_mask);
00723         images.push_back(img);
00724     }
00725 
00726     EMDeletePtr(img1d);
00727 
00728     return images;
00729 }
00730 
00731 int SVDAnalyzer::insert_image(EMData * image)
00732 {
00733         if (mask==0)
00734                 throw NullPointerException("Null mask image pointer, set_params() first");
00735 
00736         // count pixels under mask
00737         size_t totpix=mask->get_xsize()*mask->get_ysize()*mask->get_zsize();
00738         float  *d=image->get_data();
00739         float *md=mask ->get_data();
00740         for (size_t i=0,j=0; i<totpix; ++i) {
00741                 if (md[i]) {
00742                         gsl_matrix_set(A,j,nsofar,d[i]);
00743                         j++;
00744                 }
00745         }
00746         nsofar++;
00747 
00748    return 0;
00749 }
00750 #undef covmat
00751 
00752 #define eigvec(i,j) eigvec[(j)*ncov + (i)]
00753 vector<EMData*> SVDAnalyzer::analyze()
00754 {
00755 // Allocate the working space
00756 gsl_vector *work=gsl_vector_alloc(nimg);
00757 gsl_vector *S=gsl_vector_alloc(nimg);
00758 gsl_matrix *V=gsl_matrix_alloc(nimg,nimg);
00759 gsl_matrix *X=gsl_matrix_alloc(nimg,nimg);
00760 
00761 
00762 // Do the decomposition. All the real work is here
00763 gsl_linalg_SV_decomp_mod (A,X, V, S, work);
00764 //else gsl_linalg_SV_decomp_jacobi(A,V,S);
00765 
00766 vector<EMData*> ret;
00767 //unpack the results and write the output file
00768 float *md=mask->get_data();
00769 size_t totpix=mask->get_xsize()*mask->get_ysize()*mask->get_zsize();
00770 for (int k=0; k<nvec; k++) {
00771         EMData *img = new EMData;
00772         img->set_size(mask->get_xsize(),mask->get_ysize(),mask->get_zsize());
00773 
00774         float  *d=img->get_data();
00775         for (size_t i=0,j=0; i<totpix; ++i) {
00776                 if (md[i]) {
00777                         d[i]=(float)gsl_matrix_get(A,j,k);
00778                         j++;
00779                 }
00780         }
00781         img->set_attr( "eigval", gsl_vector_get(S,k));
00782         ret.push_back(img);
00783 }
00784 
00785 gsl_vector_free(work);
00786 gsl_vector_free(S);
00787 gsl_matrix_free(V);
00788 gsl_matrix_free(X);
00789 
00790 gsl_matrix_free(A);
00791 A=NULL;
00792 mask=NULL;
00793 
00794 return ret;
00795 }
00796 
00797 void SVDAnalyzer::set_params(const Dict & new_params)
00798 {
00799         params = new_params;
00800         mask = params["mask"];
00801         nvec = params["nvec"];
00802         nimg = params["nimg"];
00803 
00804         // count pixels under mask
00805         pixels=0;
00806         size_t totpix=mask->get_xsize()*mask->get_ysize()*mask->get_zsize();
00807         float *d=mask->get_data();
00808         for (size_t i=0; i<totpix; ++i) if (d[i]) ++pixels;
00809 
00810         printf("%d,%d\n",pixels,nimg);
00811         A=gsl_matrix_alloc(pixels,nimg);
00812         nsofar=0;
00813 }
00814 
00815 
00816 void EMAN::dump_analyzers()
00817 {
00818         dump_factory < Analyzer > ();
00819 }
00820 
00821 map<string, vector<string> > EMAN::dump_analyzers_list()
00822 {
00823         return dump_factory_list < Analyzer > ();
00824 }
00825 
00826 
00827 
00828 
00829 
00830 
00831 

Generated on Tue Jun 11 13:46:13 2013 for EMAN2 by  doxygen 1.3.9.1