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

utilcomm.cpp

Go to the documentation of this file.
00001 #include "mpi.h"
00002 #include "stdlib.h"
00003 #include "emdata.h"
00004 
00005 #include "ali3d_unified_mpi.h"
00006 #include "utilcomm.h"
00007 
00008 using namespace EMAN;
00009 
00010 int ReadVandBcast(MPI_Comm comm, EMData *volume, char *volfname)
00011 {
00012     int mypid, ierr, ndim, nx, ny, nz, mpierr=0;
00013     MPI_Comm_rank(comm,&mypid);
00014     FILE *fp=NULL;
00015 
00016     ierr = 0;   
00017     if (mypid == 0) {
00018         // check to see if the file exists
00019         fp = fopen(volfname,"r");
00020         if (!fp) {
00021             ierr = 1;
00022             printf("failed to open %s\n", volfname);
00023         }
00024         else {
00025             fclose(fp);
00026             volume->read_image(volfname);
00027             ndim = volume->get_ndim();
00028             nx   = volume->get_xsize();
00029             ny   = volume->get_ysize();
00030             nz   = volume->get_zsize();
00031             printf("ndim = %d, nx = %d, ny = %d, nz = %d\n", ndim, nx, ny, nz);
00032         }
00033     }
00034     mpierr = MPI_Bcast(&ierr, 1, MPI_INT, 0, comm);
00035     if (ierr==0) {
00036         mpierr = MPI_Bcast(&nx, 1, MPI_INT, 0, comm);
00037         mpierr = MPI_Bcast(&ny, 1, MPI_INT, 0, comm);
00038         mpierr = MPI_Bcast(&nz, 1, MPI_INT, 0, comm);
00039         if (mypid !=0) volume->set_size(nx,ny,nz);
00040         
00041         float * voldata = volume->get_data();
00042         mpierr = MPI_Bcast(voldata, nx*ny*nz, MPI_FLOAT, 0, comm);
00043         ierr = mpierr;
00044     }
00045     return ierr;
00046 }
00047 
00048 int ReadStackandDist(MPI_Comm comm, EMData ***images2D, char *stackfname, int *nloc)
00049 {
00050     int ncpus, mypid, ierr, mpierr=0;
00051     MPI_Status mpistatus;
00052     EMUtil *my_util;
00053     int nima; 
00054     FILE *fp=NULL;
00055 
00056     MPI_Comm_size(comm,&ncpus);
00057     MPI_Comm_rank(comm,&mypid);
00058 
00059     ierr = 0;
00060     if (mypid == 0) {
00061         fp = fopen(stackfname,"r");
00062         if (!fp) {
00063             ierr = 1;
00064             printf("failed to open %s\n", stackfname);
00065         }
00066         else {
00067             fclose(fp);
00068             nima = my_util->get_image_count(stackfname);
00069         }
00070     }
00071     mpierr = MPI_Bcast(&ierr, 1, MPI_INT, 0, comm);
00072     if (ierr == 0) {
00073         int *psize = new int[ncpus];
00074         int *nbase = new int[ncpus];
00075     
00076         mpierr = MPI_Bcast(&nima, 1, MPI_INT, 0, comm);
00077     
00078         *nloc = setpart(comm, nima, psize, nbase);
00079         *images2D = new EMData*[*nloc]; // NB!: whoever calls ReadStackandDist must delete this!
00080 
00081         EMData *img_ptr;
00082         int img_index;
00083         float *imgdata;
00084 
00085         // read the first image to get size
00086         img_ptr = new EMData();
00087         img_ptr->read_image(stackfname, 0);
00088         int nx = img_ptr->get_xsize();
00089         int ny = img_ptr->get_ysize();
00090 
00091         float s2x, s2y;
00092 
00093         if (mypid == 0) {
00094             printf("Master node reading and distributing %d images...\n", nima);
00095             for ( int ip = 0 ; ip < ncpus ; ++ip ) {
00096                 for ( int i = 0 ; i < psize[ip] ; ++i ) {
00097                     img_index = nbase[ip] + i;
00098                     if (ip != 0) {
00099                         img_ptr->read_image(stackfname, img_index);
00100                         // get a pointer to the image's data
00101                         imgdata = img_ptr->get_data();
00102                         // find the x/y shift values if it has them, otherwise set them to 0.0
00103                         try {
00104                             s2x = (*images2D)[i]->get_attr("s2x");
00105                         } catch ( std::exception& e ) {
00106                             s2x = 0.0;
00107                         }
00108                         try {
00109                             s2y = (*images2D)[i]->get_attr("s2y");
00110                         } catch ( std::exception& e ) {
00111                             s2y = 0.0;
00112                         }
00113                         // send these to processor ip
00114                         MPI_Send(imgdata, nx*ny, MPI_FLOAT, ip, ip, comm);
00115                         MPI_Send(&s2x, 1, MPI_FLOAT, ip, ip, comm);
00116                         MPI_Send(&s2y, 1, MPI_FLOAT, ip, ip, comm);
00117                     } else { // ip == 0                             
00118                         (*images2D)[i] = new EMData();
00119                         (*images2D)[i]->read_image(stackfname, img_index);
00120                         try {
00121                             s2x = (*images2D)[i]->get_attr("s2x");
00122                         } catch ( std::exception& e ) {
00123                             (*images2D)[i]->set_attr("s2x",0.0);
00124                         }
00125                         try {
00126                             s2y = (*images2D)[i]->get_attr("s2y");
00127                         } catch ( std::exception& e ) {
00128                             (*images2D)[i]->set_attr("s2y",0.0);
00129                         }
00130                     }
00131                 }
00132                 printf("finished reading data for processor %d\n", ip);
00133             }
00134         } else { // mypid != 0 : everyone else receives and reads in their data
00135             for ( int i = 0 ; i < psize[mypid] ; ++i ) {
00136                 (*images2D)[i] = new EMData();
00137                 (*images2D)[i]->set_size(nx, ny, 1);
00138                 imgdata = (*images2D)[i]->get_data();
00139                 MPI_Recv(imgdata, nx*ny, MPI_FLOAT, 0, mypid, comm, &mpistatus);
00140                 MPI_Recv(&s2x, 1, MPI_FLOAT, 0, mypid, comm, &mpistatus);
00141                 MPI_Recv(&s2y, 1, MPI_FLOAT, 0, mypid, comm, &mpistatus);
00142                 (*images2D)[i]->set_attr("s2x",s2x);
00143                 (*images2D)[i]->set_attr("s2y",s2y);
00144             }
00145             printf("received %d images for processor %d\n", *nloc, mypid);
00146         }
00147         if (mypid == 0) printf("finished reading and distributing data\n");
00148         ierr = mpierr; 
00149 
00150         EMDeletePtr(img_ptr);
00151         EMDeleteArray(psize);
00152         EMDeleteArray(nbase);
00153     }
00154     return ierr;
00155 }
00156 
00157 int ReadAngTrandDist(MPI_Comm comm, float *angleshift, char *paramfname, int nloc)
00158 {
00159     // read a parameter file that contains a list of angles and shifts and distribute 
00160     // them evenly among different processors
00161 
00162     int mypid, ierr=0, mpierr=0, nimgs=0, ncpus;
00163     double t0;
00164     MPI_Status mpistatus;
00165     FILE *fp=NULL;
00166 
00167     float *iobuffer = new float[5*nloc];
00168     if (!iobuffer) {
00169        fprintf(stderr,"failed to allocate buffer to read angles shifts\n");
00170        ierr = -1;
00171        goto EXIT;
00172     }
00173 
00174     MPI_Comm_rank(comm,&mypid);
00175     MPI_Comm_size(comm,&ncpus);
00176 
00177     t0 = MPI_Wtime();
00178 
00179     if (mypid ==0) {
00180        fp = fopen(paramfname,"r");
00181        if (!fp)  ierr = 1;
00182     }
00183     MPI_Bcast(&ierr, 1, MPI_INT, 0, comm);
00184  
00185     if ( ierr ) {
00186        if (mypid ==0) fprintf(stderr,"failed to open %s\n", paramfname);
00187     }
00188     else {
00189        if (mypid == 0) {
00190           for (int iproc = 0; iproc < ncpus; iproc++) {
00191              // figure out the number of images assigned to processor iproc
00192              if (iproc > 0) {
00193                 MPI_Recv(&nimgs, 1, MPI_INT, iproc, iproc, comm, &mpistatus);
00194                 // Read the next nimgs set of angles and shifts
00195                 for (int i = 0; i < nimgs; i++) {
00196                    fscanf(fp,"%f %f %f %f %f", 
00197                           &iobuffer[5*i+0],
00198                           &iobuffer[5*i+1],
00199                           &iobuffer[5*i+2],
00200                           &iobuffer[5*i+3],
00201                           &iobuffer[5*i+4]);
00202                 }
00203                 MPI_Send(iobuffer,5*nimgs,MPI_FLOAT,iproc,iproc,comm);
00204              }
00205              else {
00206                 for (int i = 0; i < nloc; i++) {
00207                    fscanf(fp,"%f %f %f %f %f", 
00208                           &angleshift[5*i+0],
00209                           &angleshift[5*i+1],
00210                           &angleshift[5*i+2],
00211                           &angleshift[5*i+3],
00212                           &angleshift[5*i+4]);
00213                 }
00214              }
00215           }
00216           fclose(fp);
00217        }
00218        else {
00219           // send image count to the master processor (mypid = 0)
00220           MPI_Send(&nloc, 1, MPI_INT, 0, mypid, comm);
00221           // Receive angleshifts
00222           MPI_Recv(angleshift, 5*nloc, MPI_FLOAT, 0, mypid, comm, &mpistatus);
00223        }
00224     }
00225     EMDeleteArray(iobuffer);
00226     if (mypid == 0)
00227        printf("I/O time for reading angles & shifts = %11.3e\n",
00228               MPI_Wtime() - t0);
00229 EXIT:
00230     return ierr;
00231 }
00232 
00233 int CleanStack(MPI_Comm comm, EMData ** image_stack, int nloc, int ri, Vec3i volsize, Vec3i origin)
00234 {
00235     int nx = volsize[0];
00236     int ny = volsize[1];
00237         
00238     float * rhs = new float[nx*ny];
00239     float * imgdata;
00240         
00241     // Calculate average "background" from all pixels strictly outside of radius
00242     double aba, abaloc; // average background
00243     aba = 0.0;
00244     abaloc = 0.0;
00245     int klp, klploc; // number of pixels in background
00246     klp = 0;
00247     klploc = 0;
00248         
00249     // calculate avg background in parallel
00250     for ( int i = 0 ; i < nloc ; ++i ) {
00251         imgdata = image_stack[i]->get_data();
00252         asta2(imgdata, nx, ny, ri, &abaloc, &klploc);
00253     }
00254         
00255     MPI_Allreduce(&abaloc, &aba, 1, MPI_DOUBLE, MPI_SUM, comm);
00256     MPI_Allreduce(&klploc, &klp, 1, MPI_INT, MPI_SUM, comm);
00257 
00258     aba /= klp;
00259 
00260     // subtract off the average background from pixels weakly inside of radius
00261     int x_summand, y_summand;
00262     int r_squared = ri * ri;
00263     for ( int i = 0 ; i < nloc ; ++i ) {
00264         imgdata = image_stack[i]->get_data();
00265         for ( int j = 0 ; j < nx ; ++j) {
00266             x_summand = (j - origin[0]) *  (j - origin[0]);
00267             for ( int k = 0 ; k < ny ; ++k ) {
00268                 y_summand = (k - origin[1]) *  (k - origin[1]);
00269                 if ( x_summand + y_summand <= r_squared) {
00270                     imgdata[j*ny + k] -= aba;
00271                 }
00272             }
00273         }
00274     }
00275     EMDeleteArray(rhs);
00276     return 0;
00277 }
00278 
00279 int setpart(MPI_Comm comm, int nang, int *psize, int *nbase)
00280 {
00281    int ncpus, mypid, nangloc, nrem;
00282 
00283    MPI_Comm_size(comm,&ncpus);
00284    MPI_Comm_rank(comm,&mypid);
00285 
00286    nangloc = nang/ncpus;
00287    nrem    = nang - ncpus*nangloc;
00288    if (mypid < nrem) nangloc++;
00289 
00290    for (int i = 0; i < ncpus; i++) {
00291       psize[i] = nang/ncpus;
00292       if (i < nrem) psize[i]++;
00293    }
00294  
00295    nbase[0] = 0; 
00296    for (int i = 1; i < ncpus; i++) {
00297       nbase[i] = nbase[i-1] + psize[i-1];
00298    }
00299    
00300    return nangloc;
00301 }
00302 

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