00001 #include "mpi.h"
00002 #include "stdlib.h"
00003
00004
00005 #include "emdata.h"
00006 #include "emassert.h"
00007 #include "projector.h"
00008
00009 #include "ali3d_unified_mpi.h"
00010 #include "project3d.h"
00011 #include "utilcomm.h"
00012 #include "fgcalc.h"
00013
00014
00015 #include <scitbx/lbfgsb/raw.h>
00016 #include <scitbx/array_family/shared.h>
00017
00018 using namespace EMAN;
00019 using namespace std;
00020
00021 namespace {
00022
00023 using scitbx::af::shared;
00024 using scitbx::fn::pow2;
00025 using namespace scitbx::lbfgsb::raw;
00026
00027 template <typename ElementType>
00028 ref1<ElementType>
00029 make_ref1(shared<ElementType>& a)
00030 {
00031 return ref1<ElementType>(a.begin(), a.size());
00032 }
00033 }
00034
00035
00036 #define grad(i) grad[(i)-1]
00037 int unified(MPI_Comm comm, EMData *volume, EMData **projdata,
00038 float *angleshift, int nloc, int max_iter, char *fname_base)
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
00057
00058 MPI_Comm_rank(comm,&mypid);
00059 MPI_Comm_size(comm,&ncpus);
00060
00061
00062 MPI_Allreduce(&nloc, &nang, 1, MPI_INT, MPI_SUM, comm);
00063
00064
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
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
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;
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
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
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;
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
00196 for ( int i = 0 ; i < 5 * nangloc ; ++i ) {
00197 angleshift[i] = x(nnz + 1 + nbase[mypid] + i);
00198 }
00199
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 }
00217
00218
00219
00220