#include "mpi.h"
#include "stdlib.h"
#include "emdata.h"
#include "emassert.h"
#include "projector.h"
#include "ali3d_unified_mpi.h"
#include "project3d.h"
#include "utilcomm.h"
#include "fgcalc.h"
#include <scitbx/lbfgsb/raw.h>
#include <scitbx/array_family/shared.h>
Include dependency graph for ali3d_unified_mpi.cpp:
Go to the source code of this file.
Namespaces | |
namespace | scitbx::lbfgsb::raw |
Defines | |
#define | grad(i) grad[(i)-1] |
Functions | |
template<typename ElementType> | |
ref1< ElementType > | make_ref1 (shared< ElementType > &a) |
int | unified (MPI_Comm comm, EMData *volume, EMData **projdata, float *angleshift, int nloc, int max_iter, char *fname_base) |
|
Definition at line 36 of file ali3d_unified_mpi.cpp. Referenced by fgcalc(), and recons3d_sirt_mpi(). |
|
Definition at line 29 of file ali3d_unified_mpi.cpp. Referenced by unified(). 00030 {
00031 return ref1<ElementType>(a.begin(), a.size());
00032 }
|
|
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 }
|