#include <mpi_angular_projection_matching.h>

Inheritance diagram for MpiProgAngularProjectionMatching:

Collaboration diagram for MpiProgAngularProjectionMatching:

Public Member Functions
void	read (int argc, char **argv)

	MpiProgAngularProjectionMatching ()

	~MpiProgAngularProjectionMatching ()

void	processAllImages ()

void	writeOutputFiles ()

bool	distributeJobs (size_t *imagesToSent, int node)

bool	requestJobs (std::vector< size_t > &imagesToProcess)

void	defineParams ()
	Define arguments accepted. More...

void	readParams ()
	Read arguments from command line. More...

void	produceSideInfo ()

void	computeChunks ()

void	computeChunkAngularDistance (int symmetry, int sym_order)

Public Member Functions inherited from ProgAngularProjectionMatching
void	show ()

void	run ()

void	rotationallyAlignOneImage (Matrix2D< double > &img, int imgno, int &opt_samplenr, double &opt_psi, bool &opt_flip, double &maxcorr)

void	translationallyAlignOneImage (MultidimArray< double > &img, const int &samplenr, const double &psi, const bool &opt_flip, double &opt_xoff, double &opt_yoff, double &maxcorr)

void	scaleAlignOneImage (MultidimArray< double > &img, const int &samplenr, const double &psi, const bool &opt_flip, const double &opt_xoff, const double &opt_yoff, const double &old_scale, double &opt_scale, double &maxcorr)

void	getCurrentReference (int refno, Polar_fftw_plans &local_plans)

void	processSomeImages (const std::vector< size_t > &imagesToProcess)

void	getCurrentImage (size_t imgid, Image< double > &img)

void	destroyAndClean ()

Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)

const char *	getParam (const char param, const char subparam, int arg=0)

int	getIntParam (const char *param, int arg=0)

int	getIntParam (const char param, const char subparam, int arg=0)

double	getDoubleParam (const char *param, int arg=0)

double	getDoubleParam (const char param, const char subparam, int arg=0)

float	getFloatParam (const char *param, int arg=0)

float	getFloatParam (const char param, const char subparam, int arg=0)

void	getListParam (const char *param, StringVector &list)

int	getCountParam (const char *param)

bool	checkParam (const char *param)

bool	existsParam (const char *param)

void	addParamsLine (const String &line)

void	addParamsLine (const char *line)

ParamDef *	getParamDef (const char *param) const

virtual void	quit (int exit_code=0) const

virtual int	tryRun ()

void	initProgress (size_t total, size_t stepBin=60)

void	setProgress (size_t value=0)

void	endProgress ()

void	processDefaultComment (const char param, const char left)

void	setDefaultComment (const char param, const char comment)

virtual void	initComments ()

void	setProgramName (const char *name)

void	addUsageLine (const char *line, bool verbatim=false)

void	clearUsage ()

void	addExampleLine (const char *example, bool verbatim=true)

void	addSeeAlsoLine (const char *seeAlso)

void	addKeywords (const char *keywords)

const char *	name () const

virtual void	usage (int verb=0) const

virtual void	usage (const String &param, int verb=2)

int	version () const

virtual void	show () const

virtual void	read (int argc, const char **argv, bool reportErrors=true)

virtual void	read (int argc, char **argv, bool reportErrors=true)

void	read (const String &argumentsLine)

	XmippProgram ()

	XmippProgram (int argc, const char **argv)

virtual	~XmippProgram ()

Additional Inherited Members
Public Attributes inherited from ProgAngularProjectionMatching
FileName	fn_exp

FileName	fn_ref

FileName	fn_out

FileName	fn_ctf

MetaDataDb	DFexp

MetaDataDb	DFo

size_t	dim

size_t	paddim

double	pad

double	max_shift

int	Ri

int	Ro

double	avail_memory

int	max_nr_refs_in_memory

int	max_nr_imgs_in_memory

int	total_nr_refs

int	counter_refs_in_memory

std::vector< int >	pointer_allrefs2refsinmem

std::vector< int >	pointer_refsinmem2allrefs

std::vector< size_t >	convert_refno_to_stack_position

std::vector< size_t >	ids

Polar< std::complex< double > > *	fP_ref

Polar< std::complex< double > > *	fP_img

Polar< std::complex< double > > *	fPm_img

MultidimArray< double > *	proj_ref

Polar_fftw_plans	global_plans

double *	stddev_ref

double *	stddev_img

Sampling	mysampling

bool	loop_forward_refs

int	search5d_shift

int	search5d_step

std::vector< int >	search5d_xoff

std::vector< int >	search5d_yoff

MultidimArray< double >	Mctf

bool	phase_flipped

int	threads

int	numOrientations

size_t	nr_trans

barrier_t	thread_barrier

bool	do_scale

WriteModeMetaData	do_overwrite

double	scale_step

double	scale_nsteps

int	progress_bar_step

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes inherited from XmippProgram
int	errorCode

ProgramDef *	progDef
	Program definition and arguments parser. More...

std::map< String, CommentList >	defaultComments

int	argc
	Original command line arguments. More...

const char **	argv

Detailed Description

Program to parallelize the ML 2D alignment program

Definition at line 36 of file mpi_angular_projection_matching.h.

Constructor & Destructor Documentation

◆ MpiProgAngularProjectionMatching()

MpiProgAngularProjectionMatching::MpiProgAngularProjectionMatching ( )

Constructor

Definition at line 34 of file mpi_angular_projection_matching.cpp.

 {
     imagesBuffer = nullptr;
     last_chunk = nullptr;
 }

◆ ~MpiProgAngularProjectionMatching()

MpiProgAngularProjectionMatching::~MpiProgAngularProjectionMatching ( )

Destructor

Definition at line 40 of file mpi_angular_projection_matching.cpp.

 {
     if (node->isMaster())
     {
         delete[] imagesBuffer;
         delete[] last_chunk;
     }
     delete node; //this calls MPI_Finalize
 }

Member Function Documentation

◆ computeChunkAngularDistance()

void MpiProgAngularProjectionMatching::computeChunkAngularDistance	(	int	symmetry,
		int	sym_order
	)

Definition at line 321 of file mpi_angular_projection_matching.cpp.

 {
     double non_reduntant_area_of_sphere =
         chunk_mysampling.SL.nonRedundantProjectionSphere(symmetry,
                 sym_order);
     double number_cpus = (double) node->size - 1;
     //NEXT ONE IS SAMPLING NOT ANOTHERSAMPLING
     double neighborhood_radius = fabs(acos(mysampling.cos_neighborhood_radius));
     //NEXT ONE IS SAMPLING NOT ANOTHERSAMPLING
     if (mysampling.cos_neighborhood_radius < -1.001)
         neighborhood_radius = 0;
     int counter = 0;
     while (1)
     {
         if (counter++ > 1000)
         {
             chunk_angular_distance = 0.001;
             std::cerr << "****************************************************"
             << std::endl;
             std::cerr << "* WARNING: The neighbourhood does not fit in memory "
             << std::endl;
             std::cerr << "****************************************************"
             << std::endl;
             break;
         }
         double area_chunk = non_reduntant_area_of_sphere / number_cpus;
         //area chunk is area of spheric casket=2 PI h
         chunk_angular_distance = acos(1 - area_chunk / (2 * PI));
         double area_chunck_neigh = 2 * PI
                                    * (1 - cos(chunk_angular_distance + neighborhood_radius));
         //double area_chunck= 2 * PI *( 1 - cos(chunk_angular_distance));
         //let us see how many references from the reference library fit
         //in area_chunk, that is divide area_chunk between the voronoi
         //region of the sampling points of the reference library
         double areaVoronoiRegionReferenceLibrary = 2 *( 3 *(  acos(
                     //NEXT ONE IS SAMPLING NOT ANOTHERSAMPLING
                     cos(mysampling.sampling_rate_rad)/(1+cos(mysampling.sampling_rate_rad)) )  ) - PI);
         auto number_of_images_that_fit_in_a_chunck_neigh =(int)
             ceil(area_chunck_neigh / areaVoronoiRegionReferenceLibrary);
         //#define DEBUG
 #ifdef DEBUG
 
         std::cerr << "\n\ncounter " << counter << std::endl;
         std::cerr << "area_chunk " << area_chunk << std::endl;
         std::cerr << "2*chunk_angular_distance " << 2*chunk_angular_distance << std::endl;
         //NEXT ONE IS SAMPLING NOT ANOTHERSAMPLING
         std::cerr << "sampling_rate_rad " << mysampling.sampling_rate_rad
         << " " << mysampling.sampling_rate_rad*180/PI
         << std::endl;
         std::cerr << "neighborhood_radius " << neighborhood_radius
         << std::endl;
         std::cerr << "areaVoronoiRegionReferenceLibrary " << areaVoronoiRegionReferenceLibrary << std::endl;
         std::cerr << "number_of_images_that_fit_in_a_chunck_neigh " << number_of_images_that_fit_in_a_chunck_neigh << std::endl;
         std::cerr << "number_cpus " << number_cpus << std::endl;
         std::cerr << "max_nr_imgs_in_memory " << max_nr_imgs_in_memory << std::endl;
 #endif
 #undef DEBUG
 
         if (number_of_images_that_fit_in_a_chunck_neigh > max_nr_imgs_in_memory)
             number_cpus = 1.2 * number_cpus;
         else
             break;
     }
     //chunk_angular_distance -= neighborhood_radius;
     chunk_angular_distance *= 2.0;
 
     //#define DEBUG
 #ifdef DEBUG
 
     std::cerr << "chunk_angular_distance " << chunk_angular_distance
     << std::endl
     << "neighborhood_radius " << neighborhood_radius
     << std::endl;
 #endif
 #undef DEBUG
     //chuck should not be bigger than a triangle in the icosahedra
     if (chunk_angular_distance >= 0.5 * cte_w)
         chunk_angular_distance = 0.5 * cte_w;
     chunk_angular_distance *= (180. / PI);
     //#define DEBUG
 #ifdef DEBUG
 
     std::cerr << "chunk_angular_distance_degrees " << chunk_angular_distance
     << std::endl;
 #endif
 #undef DEBUG
 
 }

◆ computeChunks()

void MpiProgAngularProjectionMatching::computeChunks ( )

These two function will be executed only by master

Definition at line 246 of file mpi_angular_projection_matching.cpp.

 {
     size_t max_number_of_images_in_around_a_sampling_point = 0;
     //process the symmetry file
     if (!chunk_mysampling.SL.isSymmetryGroup(fn_sym, symmetry, sym_order))
         REPORT_ERROR(ERR_NUMERICAL,
                      (String)"mpi_angular_proj_match::prerun Invalid symmetry: " + fn_sym);
     chunk_mysampling.SL.readSymmetryFile(fn_sym);
     // find a value for chunk_angular_distance if != -1
     if (chunk_angular_distance == -1)
         computeChunkAngularDistance(symmetry, sym_order);
     //store symmetry matrices, this is faster than computing them
     //each time we need them
     chunk_mysampling.fillLRRepository();
     int remaining_points = 0;
 
     while (remaining_points == 0)
     {
         //first set sampling rate
         chunk_mysampling.setSampling(chunk_angular_distance);
         //create sampling points in the whole sphere
         //chunk_mysampling.computeSamplingPoints(false, 0, 180);
         chunk_mysampling.computeSamplingPoints(false, 180, 0);
         //precompute product between symmetry matrices
         //and experimental data
         chunk_mysampling.fillExpDataProjectionDirectionByLR(fn_exp);
         //remove redundant sampling points: symmetry
         chunk_mysampling.removeRedundantPoints(symmetry, sym_order);
         remaining_points =
             chunk_mysampling.no_redundant_sampling_points_angles.size();
         if (chunk_angular_distance > 2)
             chunk_angular_distance -= 1;
         else
             chunk_angular_distance /= 2;
         //if(remaining_points ==0)
         if (verbose)
             std::cout << "New chunk_angular_distance " << chunk_angular_distance
             << std::endl;
         if (chunk_angular_distance < 0)
             REPORT_ERROR(ERR_VALUE_INCORRECT,
                          "Can't compute chunk_angular_distance");
     }
     //remove sampling points too far away from experimental data
     chunk_mysampling.removePointsFarAwayFromExperimentalData();
     //for each sampling point find the experimental images
     //closer to that point than to any other
     chunk_mysampling.findClosestExperimentalPoint();
     //print number of points per node
     chunk_number = chunk_mysampling.my_exp_img_per_sampling_point.size();
     for (int j = 0; j < chunk_number; j++)
     {
         if (max_number_of_images_in_around_a_sampling_point
             < chunk_mysampling.my_exp_img_per_sampling_point[j].size())
             max_number_of_images_in_around_a_sampling_point =
                 chunk_mysampling.my_exp_img_per_sampling_point[j].size();
     }
     if (verbose)
     {
         std::cout << "number of subsets: " << chunk_number << std::endl
         << "biggest subset (EXPERIMENTAL images per chunk): "
         << max_number_of_images_in_around_a_sampling_point << std::endl
         << "maximun number of references in memory: "
         << max_nr_refs_in_memory << std::endl;
     }
     //alloc memory for buffer
     if (mpi_job_size == -1)
         mpi_job_size = (int)ceil((double)DFexp.size()/(node->size - 1));
 
     //Distribution related variables
     chunk_index = 0;
     last_chunk = new int[node->size];
     for (size_t i = 1; i < node->size; ++i)
         last_chunk[i] = -1;//by default no chunk assigned yet
 }

◆ defineParams()

void MpiProgAngularProjectionMatching::defineParams ( )

virtual

Define arguments accepted.

Reimplemented from ProgAngularProjectionMatching.

Definition at line 66 of file mpi_angular_projection_matching.cpp.

 {
     ProgAngularProjectionMatching::defineParams();
     addParamsLine(
         "  [--mpi_job_size <size=10>]   : Number of images sent to a cpu in a single job ");
     addParamsLine("                                : 10 may be a good value");
     addParamsLine(
         "                                : if  -1 the computer will fill the value for you");
     addParamsLine(
         "  [--chunk_angular_distance <dist=-1>]  : sample the projection sphere with this ");
     addParamsLine(
         "                                 :using the voronoi regions");
     addParamsLine(
         "  [--sym <cn=\"c1\">]            : One of the 17 possible symmetries in");
     addParamsLine(
         "                                 :single particle electronmicroscopy");
     addParamsLine(
         "                                 :i.e.  ci, cs, cn, cnv, cnh, sn, dn, dnv,");
     addParamsLine(
         "                                 :dnh, t, td, th, o, oh, i1 (default MDB), i2, i3, i4, ih");
     addParamsLine(
         "                                 :i1h (default MDB), i2h, i3h, i4h");
     addParamsLine(
         "                                : where n may change from 1 to 99");
 }

◆ distributeJobs()

bool MpiProgAngularProjectionMatching::distributeJobs	(	size_t *	imagesToSent,
		int	node
	)

Function to distribute jobs to slaves. will return false if no more images to process. Will try to sent neighbours images to same node.

Definition at line 149 of file mpi_angular_projection_matching.cpp.

 {
     int node_index = last_chunk[node];
     static bool reached_last_chunk = false;
 
     if (node_index == -1
         || chunk_mysampling.my_exp_img_per_sampling_point[node_index].empty())
     {
 
         if (!reached_last_chunk)
         {
             node_index = chunk_index;
             chunk_index = (chunk_index + 1) % chunk_number;
             reached_last_chunk = (chunk_index == 0);
         }
         else
         {
             int i = 0;
             for (;
                  i < chunk_number
                  && chunk_mysampling.my_exp_img_per_sampling_point[chunk_index].empty();
                  ++i, chunk_index = (chunk_index + 1) % chunk_number)
                 ;
 
             //            std::cerr << formatString("DEBUG: master: assigned %d chuck to node %d, i = %d", chunk_index, node, i);
             if (i == chunk_number)
             {
                 imagesToSent[0] = 0;
                 return false;
             }
             last_chunk[node] = node_index = chunk_index;
             chunk_index = (chunk_index + 1) % chunk_number;
         }
     }
 
     size_t assigned_images = std::min((size_t)mpi_job_size, chunk_mysampling.my_exp_img_per_sampling_point[node_index].size());
     imagesToSent[0] = assigned_images;
 
     //std::cerr << formatString("DEBUG: master: assigned %lu images to node %d\n", assigned_images, node);
 
     for (size_t i = 1; i <= assigned_images; ++i)
     {
         imagesToSent[i] =
             chunk_mysampling.my_exp_img_per_sampling_point[node_index].back() + FIRST_IMAGE;
         //        std::cerr << " " << imagesToSent[i];
         chunk_mysampling.my_exp_img_per_sampling_point[node_index].pop_back();
     }
 
     //    std::cerr << std::endl;
 
     return true;
 }

◆ processAllImages()

void MpiProgAngularProjectionMatching::processAllImages ( )

virtual

Override virtual function implementations

Reimplemented from ProgAngularProjectionMatching.

Definition at line 102 of file mpi_angular_projection_matching.cpp.

 {
     if (node->isMaster())    //master distribute jobs
     {
 
         size_t finishingWorkers = 0;
         size_t processedImages = 0, finishedImages = 0, totalImages = DFexp.size();
         MPI_Status status;
 
         if (verbose)
         {
             progress_bar_step = XMIPP_MAX(1, totalImages / 80);
             init_progress_bar(totalImages);
         }
 
         while (finishingWorkers < node->size - 1)
         {
             //Receive a request of job from a worker
             //the number of images processed should be sent in imagesBuffer[0]
             MPI_Recv(&processedImages, 1, XMIPP_MPI_SIZE_T, MPI_ANY_SOURCE,
                      TAG_JOB_REQUEST, MPI_COMM_WORLD, &status);
             finishedImages += processedImages;
 
             if (!distributeJobs(imagesBuffer, status.MPI_SOURCE))
             {
                 ++finishingWorkers;
             }
             MPI_Send(imagesBuffer, imagesBuffer[0] + 1, XMIPP_MPI_SIZE_T,
                      status.MPI_SOURCE, TAG_JOB_REPLY, MPI_COMM_WORLD);
 
             if (verbose && processedImages > 0)
                 progress_bar(finishedImages);
         }
 
     } else            //slaves really work
     {
         std::vector<size_t> imagesToProcess;
         while (requestJobs (imagesToProcess))
         {
             processSomeImages(imagesToProcess);
         }
     }
 
     // Synchronize all nodes.
     //node->barrierWait();
 }

◆ produceSideInfo()

void MpiProgAngularProjectionMatching::produceSideInfo ( )

virtual

Redefine produceSideInfo

Reimplemented from ProgAngularProjectionMatching.

Definition at line 237 of file mpi_angular_projection_matching.cpp.

 {
     ProgAngularProjectionMatching::produceSideInfo();
 
     if (node->isMaster())
         computeChunks();
     node->barrierWait();
 }

◆ read()

void MpiProgAngularProjectionMatching::read	(	int	argc,
		char **	argv
	)

Redefine read

Definition at line 50 of file mpi_angular_projection_matching.cpp.

 {
     node = new MpiNode(argc, argv);
     // Master should read first
     if (node->isMaster())
         ProgAngularProjectionMatching::read(argc, (const char **)argv);
     node->barrierWait();
     if (!node->isMaster())
     {
 
         verbose = 0;//disable verbose for slaves
         ProgAngularProjectionMatching::read(argc, (const char **)argv);
     }
 }

◆ readParams()

void MpiProgAngularProjectionMatching::readParams ( )

virtual

Read arguments from command line.

Reimplemented from ProgAngularProjectionMatching.

Definition at line 93 of file mpi_angular_projection_matching.cpp.

 {
     ProgAngularProjectionMatching::readParams();
     mpi_job_size = getIntParam("--mpi_job_size");
     imagesBuffer = new size_t[mpi_job_size + 1];
     chunk_angular_distance = getDoubleParam("--chunk_angular_distance");
     fn_sym = getParam("--sym");
 }

◆ requestJobs()

bool MpiProgAngularProjectionMatching::requestJobs ( std::vector< size_t > & imagesToProcess )

Function of slaves nodes to ask jobs to master. will return false if no more images to process.

Definition at line 203 of file mpi_angular_projection_matching.cpp.

 {
     static size_t numberOfImages = 0;
     MPI_Status status;
 
     //Sent last receveiced numberOfImages, which should be processed by this worker
     MPI_Send(&numberOfImages, 1, XMIPP_MPI_SIZE_T, 0, TAG_JOB_REQUEST,
              MPI_COMM_WORLD);
     //Get new job to do
     MPI_Recv(imagesBuffer, mpi_job_size + 1, XMIPP_MPI_SIZE_T, 0,
              TAG_JOB_REPLY, MPI_COMM_WORLD, &status);
     numberOfImages = imagesBuffer[0];
 
     if (numberOfImages == 0)
         return false;
 
     imagesToProcess.clear();
     for (size_t i = 1; i <= numberOfImages; ++i)
         imagesToProcess.push_back(imagesBuffer[i]);
     return true;
 }

◆ writeOutputFiles()

void MpiProgAngularProjectionMatching::writeOutputFiles ( )

virtual

Write out results to disk This function should be override in MPI class, only master should write.

Reimplemented from ProgAngularProjectionMatching.

Definition at line 226 of file mpi_angular_projection_matching.cpp.

 {
     node->gatherMetadatas(DFo, fn_out);
     if (node->isMaster())
     {
         MetaDataDb mdAux;
         mdAux.sort(DFo,MDL_IMAGE);
         mdAux.write(fn_out, do_overwrite);
     }
 }

The documentation for this class was generated from the following files:

xmipp/libraries/parallel/mpi_angular_projection_matching.h
xmipp/libraries/parallel/mpi_angular_projection_matching.cpp

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ MpiProgAngularProjectionMatching()

◆ ~MpiProgAngularProjectionMatching()

Member Function Documentation

◆ computeChunkAngularDistance()

◆ computeChunks()

◆ defineParams()

◆ distributeJobs()

◆ processAllImages()

◆ produceSideInfo()

◆ read()

◆ readParams()

◆ requestJobs()

◆ writeOutputFiles()