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

ali3d_unified_mpi.h File Reference

#include "mpi.h"

Include dependency graph for ali3d_unified_mpi.h:

Include dependency graph

This graph shows which files directly or indirectly include this file:

Included by dependency graph

Go to the source code of this file.

Defines

#define NUMBER   float
#define PI   3.141592653589793

Functions

int unified (MPI_Comm comm, EMData *volume, EMData **projdata, float *angleshift, int nang, int max_iter, char *fname_base)


Define Documentation

#define NUMBER   float
 

Definition at line 5 of file ali3d_unified_mpi.h.

Referenced by fcalc(), fgcalc(), and unified().

#define PI   3.141592653589793
 

Definition at line 8 of file ali3d_unified_mpi.h.


Function Documentation

int unified MPI_Comm  comm,
EMData volume,
EMData **  projdata,
float *  angleshift,
int  nang,
int  max_iter,
char *  fname_base
 

Definition at line 37 of file ali3d_unified_mpi.cpp.

References cb2sph(), cord, EMDeleteArray(), fgcalc(), EMAN::EMData::get_data(), EMAN::EMData::get_ndim(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), getnnz(), ierr, imgdata, make_ref1(), nnz, nrays, NUMBER, nx, ny, PI, ptrs, rhs, setpart(), sph2cb(), EMAN::Vec3i, EMAN::EMData::write_image(), and x.

Referenced by main().

00039 {
00040     int mypid, ncpus;
00041     MPI_Status mpistatus;
00042 
00043     float *voldata, *imgdata;
00044     int   *ptrs, *cord;
00045     int   ndim, nx, ny, nz, ierr, nang, iglb, nangloc;
00046     float psi, theta, phi, sx, sy;
00047     int   nnz, nrays, ri;
00048     EMData *imgbuf;
00049     int    *psize, *nbase;
00050     char  filename[100];
00051     Vec3i volsize;
00052     Vec3i origin;
00053 
00054     float *rhs;
00055 
00056     // get parallel configuration
00057 
00058     MPI_Comm_rank(comm,&mypid);
00059     MPI_Comm_size(comm,&ncpus);
00060 
00061     // nang is the total number of images
00062     MPI_Allreduce(&nloc, &nang, 1, MPI_INT, MPI_SUM, comm);
00063 
00064     // we need psize and nbase to place distributed angleshift into the x array for optimization
00065     psize = new int[ncpus];
00066     nbase = new int[ncpus];
00067     nangloc = setpart(comm, nang, psize, nbase);
00068     if (nangloc != nloc) {
00069        printf("   ali3d_unified: nloc does not match with nangloc on Proc %d , exit...\n", mypid);
00070        ierr = MPI_Finalize();
00071        exit(1);
00072     }
00073     for (int i = 0; i < ncpus; i++) {
00074        psize[i] = psize[i]*5;
00075        nbase[i] = nbase[i]*5;
00076     }
00077 
00078     ndim = volume->get_ndim();
00079     nx   = volume->get_xsize();
00080     ny   = volume->get_ysize();
00081     nz   = volume->get_zsize();
00082 
00083     volsize[0] = nx;
00084     volsize[1] = ny;
00085     volsize[2] = nz;
00086     origin[0] = nx/2 + 1;
00087     origin[1] = ny/2 + 1;
00088     origin[2] = nz/2 + 1;
00089 
00090     voldata = volume->get_data();
00091 
00092     // the following will have to be cleaned up later...
00093     rhs = new float[nx*ny*nloc];
00094         
00095     for ( int i = 0 ; i < nloc ; ++i ) {
00096         imgdata = projdata[i]->get_data();
00097         for ( int j = 0 ; j < nx * ny ; ++j ) {
00098             rhs[nx*ny*i + j] = *(imgdata + j);
00099         }
00100     }
00101         
00102     double aba;
00103     ri = nx/2 - 1;
00104 
00105     ierr = getnnz(volsize, ri, origin, &nrays, &nnz);
00106     if (mypid ==0) printf("    nnz = %d\n", nnz);
00107         
00108     Dict  myparams;
00109 
00110     // LBFGS-B setup
00111     std::string task, csave;
00112     shared<bool> lsave_(4);
00113     ref1<bool> lsave = make_ref1(lsave_);
00114     int n=nnz+5*nang;
00115     int m=5;
00116     int iprint;
00117     shared<int> nbd_(n);
00118     ref1<int> nbd = make_ref1(nbd_);
00119     shared<int> iwa_(3*n);
00120     ref1<int> iwa = make_ref1(iwa_);
00121     shared<int> isave_(44);
00122     ref1<int> isave = make_ref1(isave_);
00123     NUMBER f, factr, pgtol;
00124     shared<NUMBER> x_(n);
00125     ref1<NUMBER> x = make_ref1(x_);
00126     shared<NUMBER> l_(n);
00127     ref1<NUMBER> l = make_ref1(l_);
00128     shared<NUMBER> u_(n);
00129     ref1<NUMBER> u = make_ref1(u_);
00130     shared<NUMBER> g_(n);
00131     ref1<NUMBER> g = make_ref1(g_);
00132     shared<NUMBER> dsave_(29);
00133     ref1<NUMBER> dsave = make_ref1(dsave_);
00134     shared<NUMBER> wa_(2*m*n+4*n+12*m*m+12*m);
00135     ref1<NUMBER> wa = make_ref1(wa_);
00136     NUMBER t1, t2;
00137     iprint = 0; // don't print out the whole initial guess
00138     factr=1.0e+1;
00139     pgtol=1.0e-5;
00140     for(int i=1;i<=n;i++) {
00141         nbd(i)=0;
00142     }
00143 
00144     // pack spherically masked volume into x
00145     ptrs = new int[nrays+1];
00146     cord = new int[3*nrays];
00147     ierr = cb2sph(voldata, volsize, ri, origin, nnz, ptrs, cord, &(x(1))); 
00148 
00149     // pack orientation parameters angleshift into x
00150     ierr = MPI_Allgatherv(angleshift, 5*nloc, MPI_FLOAT, &(x(nnz + 1)), psize, nbase, MPI_FLOAT, comm);
00151 
00152     for (int i = 0; i<nang; i++) {
00153         x((nnz+5*i+0)+1) *= (PI/180.0);
00154         x((nnz+5*i+1)+1) *= (PI/180.0);
00155         x((nnz+5*i+2)+1) *= (PI/180.0);
00156         x((nnz+5*i+3)+1) *= -1.0;
00157         x((nnz+5*i+4)+1) *= -1.0;               
00158     }
00159 
00160    
00161     if (mypid == 0) {
00162         printf(
00163             "\n"
00164             "     Structure refinement by the unified approach...\n"
00165             "\n");
00166     }
00167 
00168     int numfg = 0;
00169     char res_fname_base[64];
00170     task = "START";
00171  lbl_111:
00172     setulb(
00173         n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint,
00174         csave,lsave,isave,dsave);
00175     iprint = 1; // print out info after every iteration 
00176     if (task.substr(0,2) == "FG") {
00177 
00178         if (numfg >= max_iter) goto Terminate;
00179         sprintf(res_fname_base, "%s_fg%d" , fname_base, numfg);
00180         ierr =  fgcalc(comm, &(x(1)), volsize, nnz, nrays, origin, ri, 
00181                        ptrs, cord, &(x(nnz + 1)), nang, rhs, aba, &f, &(g(1)), res_fname_base); 
00182 
00183         numfg++;
00184         goto lbl_111;
00185     }
00186     if (task.substr(0,5) == "NEW_X") goto lbl_111;
00187  Terminate:
00188     sph2cb(&(x(1)), volsize, nrays, ri, nnz, ptrs, cord, voldata);
00189 
00190     EMUtil::ImageType WRITE_SPI = EMUtil::IMAGE_SINGLE_SPIDER;
00191     char out_fname[64];
00192     sprintf(out_fname, "vol_r_%s.spi", fname_base);
00193     if (mypid == 0) volume->write_image(out_fname, 0, WRITE_SPI);
00194 
00195     // copy new angles and shifts to angleshift from x
00196     for ( int i = 0 ; i < 5 * nangloc ; ++i ) {
00197         angleshift[i] = x(nnz + 1 + nbase[mypid] + i);
00198     }
00199     // and convert back to ali3d_d's format
00200     for ( int i = 0 ; i < nangloc ; ++i ) {
00201         angleshift[5*i + 0] *= (180.0/PI);
00202         angleshift[5*i + 1] *= (180.0/PI);
00203         angleshift[5*i + 2] *= (180.0/PI);
00204         angleshift[5*i + 3] *= -1.0;
00205         angleshift[5*i + 4] *= -1.0;            
00206     }
00207 
00208     EMDeleteArray(rhs);
00209     EMDeleteArray(ptrs);
00210     EMDeleteArray(cord);
00211     EMDeleteArray(psize);
00212     EMDeleteArray(nbase);
00213 
00214     return 0;
00215         
00216 }


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