#include "mpi.h"
#include "emdata.h"
#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include "HyBR_Cart.h"
#include "hybr.h"
#include "utilcomm_Cart.h"
#include "project3d_Cart.h"
#include "project3d.h"
Include dependency graph for HyBR_Cart.cpp:
Go to the source code of this file.
Functions | |
int | recons3d_HyBR_mpi_Cart (MPI_Comm comm_2d, MPI_Comm comm_row, MPI_Comm comm_col, EMData **images, float *angleshift, EMData *&xvol, int nangloc, int radius, int maxiter, std::string symmetry, int insolve) |
void | LBD_Cart (MPI_Comm comm_2d, MPI_Comm comm_row, MPI_Comm comm_col, float *angleshift, Vec3i volsize, int nraysloc, int nnzloc, Vec3i origin, int radius, int *ptrs, int myptrstart, int *cord, int nangloc, int nx, int ny, float *Ukloc, float **B, float **Vloc, int m, int n, int k, int maxiter, std::string symmetry) |
float | findomega (float *bhat, float *s, int m, int n, int insolve) |
float | gcvstopfun (float alpha, float *u, float *s, float beta, int k, int m, int n, int insolve) |
void | recons3d_CGLS_mpi_Cart (MPI_Comm comm_2d, MPI_Comm comm_row, MPI_Comm comm_col, EMData **images, float *angleshift, EMData *&xvol, int nangloc, int radius, int maxiter, std::string symmetry) |
|
Definition at line 528 of file HyBR_Cart.cpp. References free_vector(), myvector(), and sum_vector(). Referenced by recons3d_HyBR_mpi_Cart(). 00534 : Assume the 'optimal' regularization parameter to be the 00535 * smallest singular value. Then we take the derivative of the GCV 00536 * function with respect to alpha, evaluate it at alpha_opt, set the 00537 * derivative equal to zero and then solve for omega. 00538 * 00539 * Input: bhat - vector U'*b, where U = left singular vectors 00540 * s - vector containing the singular values 00541 * m,n - dimensions of the (projected) problem 00542 * insolve - inner solver 00543 * Output: omega - computed value for the omega parameter. 00544 */ 00545 { 00546 int i; 00547 float alpha, t0, t1, t3, t4, t5, v2, omega; 00548 float * s2 = myvector(n); 00549 float * tt = myvector(n); 00550 float * t2 = myvector(n); 00551 float * v1 = myvector(n); 00552 float * hold = myvector(m-n); 00553 float * hold2 = myvector(n); 00554 float * hold3 = myvector(n); 00555 float * hold4 = myvector(n); 00556 float * hold5 = myvector(n); 00557 float * hold6 = myvector(n); 00558 00559 if(insolve ==0){ 00560 int kopt = n; 00561 omega = (m*bhat[kopt-1]*bhat[kopt-1])/(kopt*bhat[kopt-1]*bhat[kopt-1] +2*bhat[kopt]*bhat[kopt]); 00562 } else { 00563 /* First assume the 'optimal' regularization parameter to be the smallest 00564 * singular value. 00565 */ 00566 alpha = s[n-1]; 00567 00568 /* 00569 * Compute the needed elements for the function. 00570 */ 00571 for (i = 0; i < m-n; i++) 00572 hold[i] = fabs(bhat[n+i])*fabs(bhat[n+i]); 00573 00574 t0 = sum_vector(hold, m-n); 00575 00576 for (i=0; i<n ; i++){ 00577 s2[i] = fabs(s[i]) * fabs(s[i]); 00578 tt[i] = (float)1.0 / (s2[i] + alpha*alpha); 00579 hold2[i] = s2[i]* tt[i]; 00580 } 00581 00582 t1 = sum_vector(hold2,n); 00583 for (i=0; i<n; i++){ 00584 t2[i] = fabs(bhat[i]*alpha*s[i]) * fabs(bhat[i]*alpha*s[i]); 00585 hold3[i] = t2[i] * fabs(tt[i]*tt[i]*tt[i]); 00586 hold4[i] = (s[i]*tt[i])*(s[i]*tt[i]); 00587 hold5[i] = fabs(alpha*alpha*bhat[i]*tt[i])*(fabs(alpha*alpha*bhat[i]*tt[i])); 00588 v1[i] = fabs(bhat[i]*s[i])*fabs(bhat[i]*s[i]); 00589 hold6[i]= v1[i] * fabs(tt[i]*tt[i]*tt[i]); 00590 } 00591 t3 = sum_vector(hold3,n); 00592 t4 = sum_vector(hold4,n); 00593 t5 = sum_vector(hold5,n); 00594 00595 v2 = sum_vector(hold6,n); 00596 /* 00597 * Now compute the omega value. 00598 */ 00599 omega = ((float)m*alpha*alpha*v2) / (t1*t3 + t4*(t5 + t0)); 00600 } 00601 free_vector(s2); 00602 free_vector(tt); 00603 free_vector(t2); 00604 free_vector(v1); 00605 free_vector(hold); 00606 free_vector(hold2); 00607 free_vector(hold3); 00608 free_vector(hold4); 00609 free_vector(hold5); 00610 free_vector(hold6); 00611 00612 return omega; 00613 }
|
|
Definition at line 617 of file HyBR_Cart.cpp. References free_vector(), myvector(), and sum_vector(). Referenced by recons3d_HyBR_mpi_Cart(). 00621 : alpha - regularization parameter at current iteration (k) 00622 * u - vector Ub'*(vector), where Ub = left singular vectors of BB 00623 * s - vector containing singular values of BB 00624 * beta - norm of b(right hand side) 00625 * k - current iteration 00626 * n - dimension of the original problem 00627 * insolve - inner solver 00628 * Output: G - computed value of the GCV function 00629 */ 00630 { 00631 float beta2, alpha2, num, den, G; 00632 float *s2 = myvector(k); 00633 float *t1 = myvector(k); 00634 float *t2 = myvector(k+1); 00635 float *t3 = myvector(k); 00636 int i; 00637 00638 beta2 = beta*beta; 00639 alpha2 = alpha*alpha; 00640 00641 if(insolve == 0){ 00642 for(i=(int) alpha; i<k+1; i++){ 00643 t2[i] = fabs(u[i])*fabs(u[i]); 00644 } 00645 G = (n*beta2*sum_vector(t2,k+1))/((m-alpha)*(m-alpha)); 00646 00647 }else{ 00648 for(i=0; i<k; i++) 00649 s2[i] = fabs(s[i]) * fabs(s[i]) ; 00650 00651 for (i=0; i<k; i++){ 00652 t1[i] = (float)1.0 / (s2[i] + alpha2); 00653 t2[i] = fabs(alpha2*u[i] * t1[i]) * fabs(alpha2*u[i] * t1[i]); 00654 t3[i] = s2[i] * t1[i]; 00655 } 00656 00657 num = n*beta2*(sum_vector(t2,k) + fabs(u[k])*fabs(u[k])); 00658 den = (((float)m - sum_vector(t3,k)))*(((float)m - sum_vector(t3, k))); 00659 00660 G = num / den; 00661 } 00662 free_vector(s2); 00663 free_vector(t1); 00664 free_vector(t2); 00665 free_vector(t3); 00666 00667 return G; 00668 }
|
|
Definition at line 379 of file HyBR_Cart.cpp. References bckpj3_Cart(), COL, copy_vector(), cord, dm, free_vector(), fwdpj3_Cart(), get_col(), ierr, make_proj_mat(), myptrstart, myvector(), nnzloc, nraysloc, nx, ny, phi, PI, ptrs, put_col(), ROW, sqrt(), theta, and EMAN::Vec3i. Referenced by recons3d_HyBR_mpi_Cart(). 00390 { 00391 float * vc = myvector(nnzloc); 00392 float * vc_loc = myvector(nnzloc); 00393 float * uc = myvector(nangloc*nx*ny); 00394 float * uc_loc = myvector(nangloc*nx*ny); 00395 float alpha, alpha_loc, beta, beta_loc; 00396 float * Vk = myvector(nnzloc); 00397 int ierr; 00398 float phi, theta, psi; 00399 float dm[8]; 00400 00401 MPI_Status mpistatus; 00402 int ROW = 0, COL = 1; 00403 int my2dpid; 00404 int mycoords[2], dims[2], periods[2]; 00405 int srpid, srcoords[2]; // Send/receive info 00406 00407 // Get dims and my coordinates 00408 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords); 00409 MPI_Comm_rank(comm_2d, &my2dpid); 00410 00411 int nsym = 0; 00412 Transform3D RA; 00413 Transform3D Tf; 00414 nsym = Tf.get_nsym(symmetry); 00415 Transform3D::EulerType EULER_SPIDER = Transform3D::SPIDER; 00416 Dict angdict; 00417 00418 get_col(Vloc, nnzloc, k-1, Vk); 00419 if (k == 0){ // BACKPROJECTION CODE HERE!!! 00420 for ( int i = 0 ; i < nangloc ; ++i ) { 00421 // retrieve the angles and shifts from angleshift 00422 RA = Transform3D(EULER_SPIDER, angleshift[5*i + 0], angleshift[5*i + 1], angleshift[5*i + 2]); 00423 dm[6] = angleshift[5*i + 3] * -1.0; 00424 dm[7] = angleshift[5*i + 4] * -1.0; 00425 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00426 // iterate over symmetries 00427 Tf = Tf.get_sym(symmetry, ns) * RA; 00428 angdict = Tf.get_rotation(EULER_SPIDER); 00429 00430 phi = (float) angdict["phi"] * PI/180.0; 00431 theta = (float) angdict["theta"] * PI/180.0; 00432 psi = (float) angdict["psi"] * PI/180.0; 00433 make_proj_mat(phi, theta, psi, dm); 00434 00435 ierr = bckpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, myptrstart, &Ukloc[nx*nx*i], vc_loc); 00436 } 00437 } 00438 ierr = MPI_Allreduce(vc_loc, vc, nnzloc, MPI_FLOAT, MPI_SUM, comm_col); 00439 00440 } 00441 else { //BACKPROJECTION CODE HERE!!!! 00442 for ( int i = 0 ; i < nangloc ; ++i ) { 00443 // retrieve the angles and shifts from angleshift 00444 RA = Transform3D(EULER_SPIDER, angleshift[5*i + 0], angleshift[5*i + 1], angleshift[5*i + 2]); 00445 dm[6] = angleshift[5*i + 3] * -1.0; 00446 dm[7] = angleshift[5*i + 4] * -1.0; 00447 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00448 // iterate over symmetries 00449 Tf = Tf.get_sym(symmetry, ns) * RA; 00450 angdict = Tf.get_rotation(EULER_SPIDER); 00451 00452 phi = (float) angdict["phi"] * PI/180.0; 00453 theta = (float) angdict["theta"] * PI/180.0; 00454 psi = (float) angdict["psi"] * PI/180.0; 00455 make_proj_mat(phi, theta, psi, dm); 00456 00457 ierr = bckpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, myptrstart, &Ukloc[nx*nx*i], vc_loc); 00458 } 00459 } 00460 ierr = MPI_Allreduce(vc_loc, vc, nnzloc, MPI_FLOAT, MPI_SUM, comm_col); 00461 00462 for(int i=0; i<nnzloc; i++) 00463 vc[i] = vc[i] - B[k][k-1]*Vk[i]; 00464 } 00465 00466 //alpha = norm(vc, n); 00467 alpha_loc = 0.0; 00468 for(int i=0; i<nnzloc; i++) 00469 alpha_loc += vc[i]* vc[i]; 00470 00471 MPI_Allreduce (&alpha_loc, &alpha, 1, MPI_FLOAT, MPI_SUM, comm_row); 00472 alpha = sqrt(alpha); 00473 00474 for (int i=0; i<nnzloc; i++) 00475 vc[i] = vc[i] / alpha; 00476 00477 // FORWARDPROJECTION CODE HERE!!!!! 00478 for ( int i = 0 ; i < nangloc ; ++i ) { 00479 // retrieve the angles and shifts from angleshift 00480 RA = Transform3D(EULER_SPIDER, angleshift[5*i + 0], angleshift[5*i + 1], angleshift[5*i + 2]); 00481 dm[6] = angleshift[5*i + 3] * -1.0; 00482 dm[7] = angleshift[5*i + 4] * -1.0; 00483 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00484 // iterate over symmetries 00485 Tf = Tf.get_sym(symmetry, ns) * RA; 00486 angdict = Tf.get_rotation(EULER_SPIDER); 00487 00488 phi = (float) angdict["phi"] * PI/180.0; 00489 theta = (float) angdict["theta"] * PI/180.0; 00490 psi = (float) angdict["psi"] * PI/180.0; 00491 make_proj_mat(phi, theta, psi, dm); 00492 00493 ierr = fwdpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, myptrstart, vc, &uc_loc[nx*nx*i]); 00494 } 00495 } 00496 // Now an all reduce along the rows 00497 ierr = MPI_Allreduce(uc_loc, uc, nangloc*nx*nx, MPI_FLOAT, MPI_SUM, comm_row); 00498 00499 for (int i=0; i<nangloc*nx*nx; i++) 00500 uc[i] = uc[i] - alpha*Ukloc[i]; 00501 00502 //beta = norm(uc, m); 00503 beta_loc = 0.0; 00504 for(int i=0; i<nangloc*nx*nx; i++) 00505 beta_loc += uc[i]* uc[i]; 00506 00507 MPI_Allreduce (&beta_loc, &beta, 1, MPI_FLOAT, MPI_SUM, comm_col); 00508 beta = sqrt(beta); 00509 00510 for (int i=0; i<nangloc*nx*nx; i++) 00511 uc[i] = uc[i] / beta; 00512 00513 copy_vector(uc, Ukloc,nangloc*nx*nx); 00514 put_col(Vloc, nnzloc, vc, k); 00515 00516 B[k][k] = alpha; 00517 B[k+1][k] = beta; 00518 00519 free_vector(vc); 00520 free_vector(vc_loc); 00521 free_vector(uc); 00522 free_vector(uc_loc); 00523 free_vector(Vk); 00524 } //end LBD_Cart
|
|
Definition at line 672 of file HyBR_Cart.cpp. References bckpj3_Cart(), COL, cord, dm, EMDeleteArray(), free_vector(), fwdpj3_Cart(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), getcb2sph(), getnnz(), ierr, images, make_proj_mat(), myptrstart, myvector(), nnz, nnzloc, nrays, nraysloc, nx, phi, PI, ptrs, q, ROW, EMAN::EMData::set_size(), setpart_gr1(), sph2cb(), sphpart(), sqrt(), theta, EMAN::EMData::to_zero(), EMAN::Vec3i, and x. 00679 : Conjugate Gradient for Least Squares 00680 * 00681 * Reference: A. Bjorck, "Numerical Methods for Least Squares Problems" 00682 * SIAM, 1996, pg. 289. 00683 * A method for solving large-scale, ill-posed inverse problems of the form: 00684 * b = A*x + noise 00685 * 00686 * This code was written specifically for the cryo-EM project,to be used 00687 * with MPI parallel implementation on a Cartesian virtual topology 00688 * 00689 * Input: 00690 * comm_2d, comm_row, comm_col : MPI communicators 00691 * images : 2D projection images (distributed) 00692 * angleshift : vector of angles (distributed) 00693 * nangloc : local number of 2D images 00694 * radius: used to determine spherical format 00695 * maxiter : maximum number of HyBR iterations 00696 * symmetry : symmetry parameter 00697 * 00698 * Output: xvol : CGLS solution 00699 */ 00700 { 00701 MPI_Status mpistatus; 00702 int ROW = 0, COL = 1; 00703 int my2dpid, ierr; 00704 int mycoords[2], dims[2], periods[2]; 00705 int srpid, srcoords[2]; // Send/receive info 00706 00707 // Get dims and my coordinates 00708 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords); 00709 MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology 00710 00711 int * psize; 00712 int * nbase; 00713 00714 int nangles; 00715 psize = new int[dims[ROW]]; 00716 nbase = new int[dims[ROW]]; 00717 MPI_Allreduce(&nangloc, &nangles, 1, MPI_INT, MPI_SUM, comm_col); 00718 00719 int nsym = 0; 00720 // get image size from first image 00721 int nx = images[0]->get_xsize(); 00722 00723 // make radius as large as possible if the user didn't provide one 00724 if ( radius == -1 ) radius = nx/2 - 1; 00725 00726 Vec3i volsize, origin; 00727 volsize[0] = nx; 00728 volsize[1] = nx; 00729 volsize[2] = nx; 00730 origin[0] = nx/2+1; 00731 origin[1] = nx/2+1; 00732 origin[2] = nx/2+1; 00733 00734 // vector of symmetrized angles 00735 std::vector<float> symangles(3,0.0); 00736 00737 // Now distribute the volume (in spherical format) among columns 00738 // of processors and use nnz to determine the splitting. 00739 int nrays, nnz; 00740 ierr = getnnz(volsize, radius, origin, &nrays, &nnz); 00741 00742 int nnzloc, nraysloc; 00743 int * ptrs = new int[nrays+1]; 00744 int * cord = new int[3*nrays]; 00745 int *nnzpart = new int[dims[COL]]; 00746 int *nnzbase = new int[dims[COL]+1]; 00747 int *ptrstart = new int[dims[COL]+1]; 00748 00749 ierr = getcb2sph(volsize, radius, origin, nnz, ptrs, cord); 00750 nnzloc = setpart_gr1(comm_2d, nnz, nnzpart, nnzbase); 00751 nraysloc = sphpart(comm_2d, nrays, ptrs, nnzbase, ptrstart); 00752 00753 int myptrstart = ptrstart[mycoords[COL]]; 00754 int nnzall[dims[COL]]; 00755 for (int i = 0; i<dims[COL]; i++) 00756 nnzall[i] = ptrs[ptrstart[i+1]] - ptrs[ptrstart[i]]; 00757 00758 nnzloc = nnzall[mycoords[COL]]; 00759 00760 int m = nx*nx*nangles; 00761 int n = nnz; 00762 00763 int k; 00764 float * trAb_loc = myvector(nnzloc); 00765 float * trAb = myvector(nnzloc); 00766 float * q_loc = myvector(nx*nx*nangloc); 00767 float * q = myvector(nx*nx*nangloc); 00768 float * s = myvector(nx*nx*nangloc); 00769 float * x = myvector(nnzloc); 00770 float * r_loc = myvector(nnzloc); 00771 float * r = myvector(nnzloc); 00772 float * p = myvector(nnzloc); 00773 float nrm_trAb_loc, nrm_trAb, tol, gamma_loc, gamma, beta, oldgamma, nq_loc, nq, alpha; 00774 float phi, theta, psi; 00775 float dm[8]; 00776 float * Rnrm = myvector(maxiter+1); 00777 00778 for (int i=0; i< nnzloc; i++){ 00779 trAb_loc [i] = 0.0; 00780 trAb [i] = 0.0; 00781 x[i] = 0.0; 00782 r_loc[i] = 0.0; 00783 r[i] = 0.0; 00784 p[i] = 0.0; 00785 } 00786 for (int i=0; i< nx*nx*nangloc; i++){ 00787 q_loc[i] = 0.0; 00788 q[i] = 0.0; 00789 s[i] = 0.0; 00790 } 00791 00792 EMData * current_image; 00793 Transform3D RA; 00794 Transform3D Tf; 00795 nsym = Tf.get_nsym(symmetry); 00796 Transform3D::EulerType EULER_SPIDER = Transform3D::SPIDER; 00797 Dict angdict; 00798 00799 float * image_data; 00800 00801 // trAb = A'*b; 00802 for (int i=0; i<nangloc; i++){ 00803 current_image = images[i]; 00804 image_data = current_image->get_data(); 00805 // retrieve the angles and shifts associated with each image from the array 00806 // angleshift. 00807 phi = angleshift[5*i + 0]; 00808 theta = angleshift[5*i + 1]; 00809 psi = angleshift[5*i + 2]; 00810 dm[6] = angleshift[5*i + 3] * -1.0; 00811 dm[7] = angleshift[5*i + 4] * -1.0; 00812 // make an instance of Transform3D with the angles 00813 RA = Transform3D(EULER_SPIDER, phi, theta, psi); 00814 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00815 // compose it with each symmetry rotation in turn 00816 // shifts (stored in dm[6] and dm[7] remain fixed 00817 Tf = Tf.get_sym(symmetry, ns) * RA; 00818 angdict = Tf.get_rotation(EULER_SPIDER); 00819 phi = (float) angdict["phi"] * PI/180.0; 00820 theta = (float) angdict["theta"] * PI/180.0; 00821 psi = (float) angdict["psi"] * PI/180.0; 00822 make_proj_mat(phi, theta, psi, dm); 00823 // accumulate the back-projected images in bvol_loc 00824 ierr = bckpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, 00825 myptrstart, image_data, trAb_loc); 00826 } 00827 } 00828 // Now an all reduce along the columns 00829 ierr = MPI_Allreduce (trAb_loc, trAb, nnzloc, MPI_FLOAT, MPI_SUM, comm_col); 00830 00831 nrm_trAb_loc = 0.0; 00832 for(int i=0; i<nnzloc; i++) 00833 nrm_trAb_loc += trAb[i]*trAb[i]; 00834 00835 ierr = MPI_Allreduce (&nrm_trAb_loc, &nrm_trAb, 1, MPI_FLOAT, MPI_SUM, comm_row); 00836 00837 nrm_trAb = sqrt(nrm_trAb); 00838 tol = sqrt(FLT_EPSILON)*nrm_trAb; 00839 00840 // s = A*x 00841 // The following code is needed if we have an initial guess for x. 00842 /*for (int i=0; i<nangs; i++){ 00843 phi = angles[5*i+0] * PI/180.0; 00844 theta = angles[5*i+1] * PI/180.0; 00845 psi = angles[5*i+2] * PI/180.0; 00846 dm[6] = angles[5*i+3] * -1.0; 00847 dm[7] = angles[5*i+4] * -1.0; 00848 00849 make_proj_mat(phi, theta, psi, dm); 00850 00851 ierr = fwdpj3(volsize, nrays, nnz, dm, origin, radius, ptrs, cord, x, &s[nx*ny*i]); 00852 } 00853 */ 00854 //s = b - s; 00855 for (int i=0; i<nangloc; i++){ 00856 current_image = images[i]; 00857 image_data = current_image->get_data(); 00858 for (int j=0; j<nx*nx; j++) 00859 s[i*nx*nx+j] = image_data[j]- s[i*nx*nx+j]; 00860 } 00861 00862 //r = A'*s; 00863 for (int i=0; i<nangloc; i++){ 00864 phi = angleshift[5*i + 0]; 00865 theta = angleshift[5*i + 1]; 00866 psi = angleshift[5*i + 2]; 00867 dm[6] = angleshift[5*i + 3] * -1.0; 00868 dm[7] = angleshift[5*i + 4] * -1.0; 00869 // make an instance of Transform3D with the angles 00870 RA = Transform3D(EULER_SPIDER, phi, theta, psi); 00871 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00872 // compose it with each symmetry rotation in turn 00873 // shifts (stored in dm[6] and dm[7] remain fixed 00874 Tf = Tf.get_sym(symmetry, ns) * RA; 00875 angdict = Tf.get_rotation(EULER_SPIDER); 00876 phi = (float) angdict["phi"] * PI/180.0; 00877 theta = (float) angdict["theta"] * PI/180.0; 00878 psi = (float) angdict["psi"] * PI/180.0; 00879 make_proj_mat(phi, theta, psi, dm); 00880 // accumulate the back-projected images in bvol_loc 00881 ierr = bckpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, 00882 myptrstart, &s[nx*nx*i], r_loc); 00883 } 00884 } 00885 // Now an all reduce along the columns 00886 ierr = MPI_Allreduce (r_loc, r, nnzloc, MPI_FLOAT, MPI_SUM, comm_col); 00887 00888 gamma_loc = 0.0; 00889 for(int i=0; i<nnzloc; i++) 00890 gamma_loc += r[i]*r[i]; 00891 00892 ierr = MPI_Allreduce (&gamma_loc, &gamma, 1, MPI_FLOAT, MPI_SUM, comm_row); 00893 Rnrm[0] = sqrt(gamma); 00894 if (my2dpid==0) printf("Iteration %d, Rnorm = %f\n", 1, Rnrm[0]/nrm_trAb); 00895 00896 k = 0; 00897 /* 00898 * Main CGLS loop. 00899 */ 00900 double t0; 00901 t0 = MPI_Wtime(); 00902 while (Rnrm[k]>tol && k<= maxiter-1){ 00903 k++; 00904 if(k==1){ 00905 for(int i=0;i<nnzloc; i++) 00906 p[i] = r[i]; 00907 }else{ 00908 beta = gamma / oldgamma; 00909 for(int i=0;i<nnzloc; i++) 00910 p[i] = r[i]+beta*p[i]; 00911 } 00912 // q = A*p; 00913 for (int i=0; i < nx*nx*nangloc; i++) 00914 q_loc[i] = 0.0; 00915 00916 for (int i=0; i < nangloc; i++){ 00917 // retrieve the angles and shifts from angleshift 00918 RA = Transform3D(EULER_SPIDER, angleshift[5*i + 0], angleshift[5*i + 1], angleshift[5*i + 2]); 00919 dm[6] = angleshift[5*i + 3] * -1.0; 00920 dm[7] = angleshift[5*i + 4] * -1.0; 00921 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00922 // iterate over symmetries 00923 Tf = Tf.get_sym(symmetry, ns) * RA; 00924 angdict = Tf.get_rotation(EULER_SPIDER); 00925 00926 phi = (float) angdict["phi"] * PI/180.0; 00927 theta = (float) angdict["theta"] * PI/180.0; 00928 psi = (float) angdict["psi"] * PI/180.0; 00929 make_proj_mat(phi, theta, psi, dm); 00930 00931 ierr = fwdpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, 00932 myptrstart, p, &q_loc[nx*nx*i]); 00933 } 00934 } 00935 // Now an all reduce along the rows 00936 ierr = MPI_Allreduce(q_loc, q, nangloc*nx*nx, MPI_FLOAT, MPI_SUM, comm_row); 00937 00938 nq_loc = 0.0; 00939 for(int i=0; i< nangloc*nx*nx; i++) 00940 nq_loc += q[i]*q[i]; 00941 00942 MPI_Allreduce (&nq_loc, &nq, 1, MPI_FLOAT, MPI_SUM, comm_col); 00943 00944 alpha = gamma/nq; 00945 00946 for(int i=0; i<nnzloc; i++){ 00947 x[i] = x[i]+alpha*p[i]; 00948 } 00949 00950 for(int i=0; i<nangloc*nx*nx; i++) 00951 s[i] = s[i]- alpha*q[i]; 00952 00953 //r = A'*s; 00954 for(int i=0; i< nnzloc; i++) 00955 r_loc[i] = 0.0; 00956 00957 for (int i=0; i<nangloc; i++){ 00958 // retrieve the angles and shifts from angleshift 00959 RA = Transform3D(EULER_SPIDER, angleshift[5*i + 0], angleshift[5*i + 1], angleshift[5*i + 2]); 00960 dm[6] = angleshift[5*i + 3] * -1.0; 00961 dm[7] = angleshift[5*i + 4] * -1.0; 00962 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00963 // iterate over symmetries 00964 Tf = Tf.get_sym(symmetry, ns) * RA; 00965 angdict = Tf.get_rotation(EULER_SPIDER); 00966 00967 phi = (float) angdict["phi"] * PI/180.0; 00968 theta = (float) angdict["theta"] * PI/180.0; 00969 psi = (float) angdict["psi"] * PI/180.0; 00970 make_proj_mat(phi, theta, psi, dm); 00971 00972 ierr = bckpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, 00973 myptrstart, &s[nx*nx*i], r_loc); 00974 } 00975 } 00976 ierr = MPI_Allreduce(r_loc, r, nnzloc, MPI_FLOAT, MPI_SUM, comm_col); 00977 00978 oldgamma = gamma; 00979 gamma_loc = 0.0; 00980 for(int i=0; i< nnzloc; i++) 00981 gamma_loc += r[i]*r[i]; 00982 00983 MPI_Allreduce (&gamma_loc, &gamma, 1, MPI_FLOAT, MPI_SUM, comm_row); 00984 00985 Rnrm[k] = sqrt(gamma); 00986 if(k!=1) 00987 if (my2dpid ==0) printf("Iteration %d, Rnorm = %f\n", k, Rnrm[k]/nrm_trAb); 00988 } 00989 if (my2dpid == 0)printf("Exit CGLS: time = %11.3e\n", MPI_Wtime()-t0); 00990 // Bring all parts of spherical volume together and turn it into cube format. 00991 if (mycoords[ROW] == 0 && mycoords[COL] != 0 ){ 00992 srcoords[ROW] = 0; 00993 srcoords[COL] = 0; 00994 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00995 MPI_Send(x, nnzloc, MPI_FLOAT, srpid, my2dpid, comm_2d); 00996 } 00997 00998 if (mycoords[ROW] == 0 && mycoords[COL] == 0 ){ 00999 xvol->set_size(nx, nx, nx); 01000 xvol->to_zero(); 01001 float * voldata = xvol->get_data(); 01002 float * xvol_sph = new float[nnz]; 01003 01004 for(int i=0; i< dims[COL]; i++){ 01005 if (i==0){ 01006 for(int i=0; i<nnzloc; i++) xvol_sph[i] = x[i]; 01007 } 01008 else { 01009 srcoords[ROW] = 0; 01010 srcoords[COL] = i; 01011 MPI_Cart_rank(comm_2d, srcoords, &srpid); 01012 MPI_Recv(&xvol_sph[ptrs[ptrstart[i]]-1], nnzall[i], MPI_FLOAT, srpid, srpid, 01013 comm_2d, &mpistatus); 01014 } 01015 } 01016 // unpack the spherical volume back out into the original EMData object 01017 ierr = sph2cb(xvol_sph, volsize, nrays, radius, nnz, ptrs, cord, voldata); 01018 EMDeleteArray(xvol_sph); 01019 } 01020 01021 EMDeleteArray(ptrs); 01022 EMDeleteArray(cord); 01023 EMDeleteArray(psize); 01024 EMDeleteArray(nbase); 01025 EMDeleteArray(nnzpart); 01026 EMDeleteArray(nnzbase); 01027 EMDeleteArray(ptrstart); 01028 01029 free_vector(x); 01030 free_vector(trAb_loc); 01031 free_vector(trAb); 01032 free_vector(s); 01033 free_vector(q_loc); 01034 free_vector(q); 01035 free_vector(r_loc); 01036 free_vector(r); 01037 free_vector(p); 01038 free_vector(Rnrm); 01039 }
|
|
Definition at line 19 of file HyBR_Cart.cpp. References COL, copy_vector(), cord, dm, EMDeleteArray(), findomega(), free_matrix(), free_vector(), gcvstopfun(), EMAN::EMData::get_data(), get_row(), EMAN::EMData::get_xsize(), getcb2sph(), getnnz(), ierr, images, LBD_Cart(), matrix(), matvec_mult(), matvec_multT(), min, myptrstart, myvector(), nnz, nnzloc, nrays, nraysloc, nx, ptrs, ROW, EMAN::EMData::set_size(), setpart_gr1(), sph2cb(), sphpart(), sqrt(), submatrix(), subvector(), sum_vector(), svd(), tikhonov(), EMAN::EMData::to_zero(), tsvd(), EMAN::Vec3i, and x. Referenced by main(). 00030 : 00031 * b = A*x + noise 00032 * The method combines an iterative Lanczos Bidiagonalization (LBD) Method 00033 * with a SVD-based regularization method to stabilize the semiconvergence 00034 * behavior that is characteristic of many ill-posed problems. 00035 * 00036 * Note: Here we use a "stripped-down" version where all options are set 00037 * to default. 00038 * 00039 * This code was written specifically for the cryo-EM project,to be used 00040 * with MPI parallel implementation on a Cartesian virtual topology 00041 * 00042 * Input: 00043 * comm_2d, comm_row, comm_col : MPI communicators 00044 * images : 2D projection images (distributed) 00045 * angleshift : vector of angles (distributed) 00046 * nangloc : local number of 2D images 00047 * radius: used to determine spherical format 00048 * maxiter : maximum number of HyBR iterations 00049 * symmetry : symmetry parameter 00050 * insolve : inner solver (0=TSVD, 1=Tikhonov) 00051 * 00052 * Output: xvol : HyBR solution 00053 */ 00054 { 00055 MPI_Status mpistatus; 00056 int ROW = 0, COL = 1; 00057 int my2dpid, ierr; 00058 int mycoords[2], dims[2], periods[2]; 00059 int srpid, srcoords[2]; // Send/receive info 00060 00061 // Get dims and my coordinates 00062 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords); 00063 MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology 00064 00065 int * psize; 00066 int * nbase; 00067 00068 int nangles; 00069 psize = new int[dims[ROW]]; 00070 nbase = new int[dims[ROW]]; 00071 MPI_Allreduce(&nangloc, &nangles, 1, MPI_INT, MPI_SUM, comm_col); 00072 00073 int nsym = 0; 00074 // get image size from first image 00075 int nx = images[0]->get_xsize(); 00076 if ( radius == -1 ) radius = nx/2 - 1; // make radius as large as possible if the user didn't provide one 00077 00078 Vec3i volsize, origin; 00079 volsize[0] = nx; 00080 volsize[1] = nx; 00081 volsize[2] = nx; 00082 origin[0] = nx/2+1; 00083 origin[1] = nx/2+1; 00084 origin[2] = nx/2+1; 00085 00086 // vector of symmetrized angles 00087 std::vector<float> symangles(3,0.0); 00088 00089 // Now distribute the volume (in spherical format) among columns 00090 // of processors and use nnz to determine the splitting. 00091 int nrays, nnz; 00092 ierr = getnnz(volsize, radius, origin, &nrays, &nnz); 00093 00094 int nnzloc, nraysloc; 00095 int * ptrs = new int[nrays+1]; 00096 int * cord = new int[3*nrays]; 00097 int *nnzpart = new int[dims[COL]]; 00098 int *nnzbase = new int[dims[COL]+1]; 00099 int *ptrstart = new int[dims[COL]+1]; 00100 00101 ierr = getcb2sph(volsize, radius, origin, nnz, ptrs, cord); 00102 nnzloc = setpart_gr1(comm_2d, nnz, nnzpart, nnzbase); 00103 nraysloc = sphpart(comm_2d, nrays, ptrs, nnzbase, ptrstart); 00104 00105 int myptrstart = ptrstart[mycoords[COL]]; 00106 int nnzall[dims[COL]]; 00107 for (int i = 0; i<dims[COL]; i++) 00108 nnzall[i] = ptrs[ptrstart[i+1]] - ptrs[ptrstart[i]]; 00109 00110 nnzloc = nnzall[mycoords[COL]]; 00111 00112 int m = nx*nx*nangles; 00113 int n = nnz; 00114 00115 int iterations_save, warning; 00116 float w, omega,curmin, g, alpha, beta_loc, beta; 00117 00118 float * Ukloc = myvector(nx*nx*nangloc); 00119 float ** Vloc = matrix(nnzloc, maxiter+1); 00120 float ** BB = matrix(maxiter+2, maxiter+1); 00121 00122 float ** Ub = matrix(maxiter+2, maxiter+2); 00123 float ** Vb = matrix(maxiter+1, maxiter+1); 00124 float * Sb = myvector(maxiter+1); 00125 00126 float * vec = myvector(maxiter+2); 00127 float * f = myvector(maxiter+2); 00128 float * Omegas = myvector(maxiter+1); 00129 float * GCV = myvector(maxiter+1); 00130 00131 float * x_save = myvector(nnzloc); 00132 float * x = myvector(nnzloc); 00133 00134 float * UbTvector; 00135 float ** BB2; 00136 float ** Ub2; 00137 float * Sb2, *f2, *Utb; 00138 float ** Vb2; 00139 float * vec2; 00140 float * row1Ub2 = myvector(maxiter+2); 00141 00142 for (int i=0; i< nnzloc; i++){ 00143 x[i] = 0.0; 00144 x_save[i] = 0.0; 00145 for(int j=0; j<maxiter+1; j++) 00146 Vloc[i][j] = 0.0; 00147 } 00148 00149 EMData * current_image; 00150 float phi, theta, psi; 00151 Transform3D RA; 00152 Transform3D Tf; 00153 nsym = Tf.get_nsym(symmetry); 00154 Transform3D::EulerType EULER_SPIDER = Transform3D::SPIDER; 00155 Dict angdict; 00156 00157 float * image_data; 00158 float dm[8]; 00159 00160 // beta = norm(b, m); 00161 beta_loc = 0.0; 00162 for (int i=0; i<nangloc; i++){ 00163 current_image = images[i]; 00164 image_data = current_image->get_data(); 00165 00166 for(int j=0; j<nx*nx; j++) 00167 beta_loc += image_data[j] * image_data[j]; 00168 } 00169 ierr = MPI_Allreduce (&beta_loc, &beta, 1, MPI_FLOAT, MPI_SUM, comm_col); 00170 beta = sqrt(beta); 00171 00172 for (int i=0; i<nangloc; i++){ 00173 current_image = images[i]; 00174 image_data = current_image->get_data(); 00175 00176 for(int j=0; j<nx*nx; j++) 00177 Ukloc[nx*nx*i+j] = image_data[j] / beta; 00178 } 00179 00180 vec[0] = beta; 00181 warning = 0; 00182 00183 double t0; 00184 t0 = MPI_Wtime(); 00185 for (int i = 0; i< maxiter+1; i++){ 00186 LBD_Cart(comm_2d, comm_row, comm_col, angleshift, volsize, nraysloc, nnzloc, origin, radius, ptrs, myptrstart, cord,nangloc, nx, nx, Ukloc, BB, Vloc, m, n, i, maxiter, symmetry); 00187 00188 if (i+1 >= 2){ // Use inner solver and adaptive GCV to solve projected problem 00189 BB2 = submatrix(BB, maxiter+2, maxiter+1, i+2, i+1); 00190 Ub2 = submatrix(Ub, maxiter+2, maxiter+2, i+2, i+2); 00191 Vb2 = submatrix(Vb, maxiter+1, maxiter+1, i+1,i+1); 00192 Sb2 = subvector(Sb, maxiter+1, i+1); 00193 f2 = subvector(f, maxiter+2 ,i+1); 00194 vec2 = subvector(vec, maxiter+2, i+2); 00195 00196 svd(BB2, i+2, i+1, Ub2, Sb2, Vb2); 00197 00198 /*Use the adaptive, modified GCV method*/ 00199 UbTvector = myvector(i+2); 00200 matvec_multT(Ub2,vec2, i+2, i+2, UbTvector); 00201 omega = findomega(UbTvector, Sb2, i+2, i+1, insolve); 00202 Omegas[i-1] = min((float)1.0, omega); 00203 00204 w = sum_vector(Omegas, i)/((float)(i+1) -(float)1.0); 00205 00206 alpha = -(float)1.0; 00207 //Solve the projected problem 00208 if(insolve == 0) 00209 tsvd(Ub2, Sb2, Vb2, vec2, i+2, i+1, w, f2, &alpha); 00210 else 00211 tikhonov(Ub2, Sb2, Vb2, vec2, i+2, i+1, w, f2, &alpha); 00212 00213 /*Compute the GCV value used to find the stopping criteria*/ 00214 get_row(Ub2, i+2, 0, row1Ub2); 00215 g = gcvstopfun(alpha, row1Ub2, Sb2, beta, i+1, m, n, insolve); 00216 00217 GCV[i-1] = g; 00218 00219 /*Determine if GCV wants us to stop*/ 00220 if (i+1 > 2){ 00221 /*-------- GCV curve is flat, we stop ------------------------*/ 00222 if (fabs((GCV[i-1]-GCV[i-2]))/GCV[0] < (pow((float)10.0, - (float)6.0))){ 00223 matvec_mult(Vloc, f2, nnzloc, i+1, x); 00224 // Return volume in cube format. 00225 if (mycoords[ROW] == 0 && mycoords[COL] != 0 ){ 00226 srcoords[ROW] = 0; 00227 srcoords[COL] = 0; 00228 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00229 00230 MPI_Send(x, nnzloc, MPI_FLOAT, srpid, my2dpid, comm_2d); 00231 } 00232 00233 if (mycoords[ROW] == 0 && mycoords[COL] == 0 ){ 00234 printf("Exit Code 1: Time = %11.3e\n", MPI_Wtime()-t0); 00235 xvol->set_size(nx, nx, nx); 00236 xvol->to_zero(); 00237 float * voldata = xvol->get_data(); 00238 float * xvol_sph = new float[nnz]; 00239 00240 for(int i=0; i< dims[COL]; i++){ 00241 if (i==0){ 00242 for(int i=0; i<nnzloc; i++) 00243 xvol_sph[i] = x[i]; 00244 }else{ 00245 srcoords[ROW] = 0; 00246 srcoords[COL] = i; 00247 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00248 00249 MPI_Recv(&xvol_sph[ptrs[ptrstart[i]]-1], nnzall[i], MPI_FLOAT, srpid, srpid, comm_2d, &mpistatus); 00250 } 00251 } 00252 00253 // unpack the spherical volume back out into the original EMData object 00254 ierr = sph2cb(xvol_sph, volsize, nrays, radius, nnz, ptrs, cord, voldata); 00255 EMDeleteArray(xvol_sph); 00256 } 00257 return 0; 00258 } 00259 /*--- Have warning: Avoid bumps by using a window of 4 its --*/ 00260 else if (warning && i > iterations_save + 3){ 00261 curmin = GCV[iterations_save]; 00262 for(int j = 1; j< 3; j++){ 00263 if (GCV[iterations_save + j]< curmin){ 00264 curmin = GCV[iterations_save+ j]; 00265 } 00266 } 00267 if (GCV[iterations_save] <= curmin){ 00268 /* We should have stopped at iterations_save.*/ 00269 copy_vector(x_save, x, nnzloc); 00270 00271 // Return volume in cube format. 00272 if (mycoords[ROW] == 0 && mycoords[COL] != 0 ){ 00273 srcoords[ROW] = 0; 00274 srcoords[COL] = 0; 00275 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00276 00277 MPI_Send(x, nnzloc, MPI_FLOAT, srpid, my2dpid, comm_2d); 00278 } 00279 00280 if (mycoords[ROW] == 0 && mycoords[COL] == 0 ){ 00281 printf("Exit Code 2: Time = %11.3e\n", MPI_Wtime()-t0); 00282 xvol->set_size(nx, nx, nx); 00283 xvol->to_zero(); 00284 float * voldata = xvol->get_data(); 00285 float * xvol_sph = new float[nnz]; 00286 00287 for(int i=0; i< dims[COL]; i++){ 00288 if (i==0){ 00289 for(int i=0; i<nnzloc; i++) 00290 xvol_sph[i] = x[i]; 00291 }else{ 00292 srcoords[ROW] = 0; 00293 srcoords[COL] = i; 00294 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00295 00296 MPI_Recv(&xvol_sph[ptrs[ptrstart[i]]-1], nnzall[i], MPI_FLOAT, srpid, srpid, comm_2d, &mpistatus); 00297 } 00298 } 00299 // unpack the spherical volume back out into the original EMData object 00300 ierr = sph2cb(xvol_sph, volsize, nrays, radius, nnz, ptrs, cord, voldata); 00301 EMDeleteArray(xvol_sph); 00302 } 00303 return 0; 00304 } 00305 else { 00306 /* It was just a bump... keep going*/ 00307 warning = 0; 00308 iterations_save = maxiter; 00309 } 00310 } 00311 /* ----- No warning: Check GCV function-----------------------*/ 00312 else if ( ! warning ){ 00313 if (GCV[i-2] < GCV[i-1]){ 00314 warning = 1; 00315 matvec_mult(Vloc, f2,nnzloc, i+1, x_save); 00316 iterations_save = i; 00317 } 00318 } 00319 00320 free_matrix(BB2); 00321 free_matrix(Ub2); 00322 free_matrix(Vb2); 00323 free_vector(Sb2); 00324 free_vector(vec2); 00325 } 00326 matvec_mult(Vloc, f2, nnzloc, i+1, x); 00327 } 00328 if (my2dpid ==0) printf("Iteration %d, alpha = %f, omega is %f\n", i+1, alpha, w); 00329 } //end HyBR main loop 00330 00331 // Bring all parts of spherical volume together and turn it into cube format. 00332 if (mycoords[ROW] == 0 && mycoords[COL] != 0 ){ 00333 srcoords[ROW] = 0; 00334 srcoords[COL] = 0; 00335 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00336 00337 MPI_Send(x, nnzloc, MPI_FLOAT, srpid, my2dpid, comm_2d); 00338 } 00339 00340 if (mycoords[ROW] == 0 && mycoords[COL] == 0 ){ 00341 printf("Exit Code 3: Time = %11.3e\n", MPI_Wtime()-t0); 00342 00343 xvol->set_size(nx, nx, nx); 00344 xvol->to_zero(); 00345 float * voldata = xvol->get_data(); 00346 float * xvol_sph = new float[nnz]; 00347 00348 for(int i=0; i< dims[COL]; i++){ 00349 if (i==0){ 00350 for(int i=0; i<nnzloc; i++) 00351 xvol_sph[i] = x[i]; 00352 }else{ 00353 srcoords[ROW] = 0; 00354 srcoords[COL] = i; 00355 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00356 00357 MPI_Recv(&xvol_sph[ptrs[ptrstart[i]]-1], nnzall[i], MPI_FLOAT, srpid, srpid, comm_2d, &mpistatus); 00358 } 00359 } 00360 // unpack the spherical volume back out into the original EMData object 00361 ierr = sph2cb(xvol_sph, volsize, nrays, radius, nnz, ptrs, cord, voldata); 00362 EMDeleteArray(xvol_sph); 00363 } 00364 00365 EMDeleteArray(ptrs); 00366 EMDeleteArray(cord); 00367 EMDeleteArray(psize); 00368 EMDeleteArray(nbase); 00369 EMDeleteArray(nnzpart); 00370 EMDeleteArray(nnzbase); 00371 EMDeleteArray(ptrstart); 00372 00373 return 0; // recons3d_HyBR_mpi_Cart 00374 }
|