bioem.cpp 29.8 KB
Newer Older
Pilar Cossio's avatar
License  
Pilar Cossio committed
1
2
3
4
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        < BioEM software for Bayesian inference of Electron Microscopy images>
            Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer.
            Max Planck Institute of Biophysics, Frankfurt, Germany.
5

Pilar Cossio's avatar
License  
Pilar Cossio committed
6
7
8
9
                See license statement for terms of distribution.

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

10
11
12
13
14
15
16
17
#include <mpi.h>

#define MPI_CHK(expr) \
	if (expr != MPI_SUCCESS) \
	{ \
		fprintf(stderr, "Error in MPI function %s: %d\n", __FILE__, __LINE__); \
	}

18
19
20
21
22
23
24
25
26
#include <fstream>
#include <boost/program_options.hpp>
#include <iostream>
#include <algorithm>
#include <iterator>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cmath>
27

28
#ifdef WITH_OPENMP
29
#include <omp.h>
30
#endif
31
32
33
34
35
36
37
38
39
40

#include <fftw3.h>
#include <math.h>
#include "cmodules/timer.h"

#include "param.h"
#include "bioem.h"
#include "model.h"
#include "map.h"

41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#ifdef BIOEM_USE_NVTX
#include "nvToolsExt.h"

const uint32_t colors[] = { 0x0000ff00, 0x000000ff, 0x00ffff00, 0x00ff00ff, 0x0000ffff, 0x00ff0000, 0x00ffffff };
const int num_colors = sizeof(colors)/sizeof(colors[0]);

#define cuda_custom_timeslot(name,cid) { \
	int color_id = cid; \
	color_id = color_id%num_colors;\
	nvtxEventAttributes_t eventAttrib = {0}; \
	eventAttrib.version = NVTX_VERSION; \
	eventAttrib.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE; \
	eventAttrib.colorType = NVTX_COLOR_ARGB; \
	eventAttrib.color = colors[color_id]; \
	eventAttrib.messageType = NVTX_MESSAGE_TYPE_ASCII; \
	eventAttrib.message.ascii = name; \
	nvtxRangePushEx(&eventAttrib); \
}
#define cuda_custom_timeslot_end nvtxRangePop();
#else
#define cuda_custom_timeslot(name,cid)
#define cuda_custom_timeslot_end
#endif
64

65
66
67
68
69
70
71
72
73
74
75
#include "bioem_algorithm.h"

using namespace boost;
namespace po = boost::program_options;

using namespace std;

// A helper function of Boost
template<class T>
ostream& operator<<(ostream& os, const vector<T>& v)
{
76
77
	copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
	return os;
78
79
80
81
}

bioem::bioem()
{
82
	FFTAlgo = getenv("FFTALGO") == NULL ? 1 : atoi(getenv("FFTALGO"));
83
	DebugOutput = getenv("BIOEM_DEBUG_OUTPUT") == NULL ? 2 : atoi(getenv("BIOEM_DEBUG_OUTPUT"));
84
85
86
87
88
89
90
91
}

bioem::~bioem()
{
}

int bioem::configure(int ac, char* av[])
{
David Rohr's avatar
David Rohr committed
92
93
94
95
96
	// **************************************************************************************
	// **** Configuration Routine using boost for extracting parameters, models and maps ****
	// **************************************************************************************
	// ****** And Precalculating necessary grids, map crosscorrelations and kernels  ********
	// *************************************************************************************
97

David Rohr's avatar
David Rohr committed
98
99
	HighResTimer timer;

David Rohr's avatar
David Rohr committed
100
101
102
103
104
	if (mpi_rank == 0)
	{
		// *** Inizialzing default variables ***
		std::string infile, modelfile, mapfile;
		Model.readPDB = false;
105
		param.param_device.writeAngles = false;
David Rohr's avatar
David Rohr committed
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
		param.dumpMap = false;
		param.loadMap = false;
		RefMap.readMRC = false;
		RefMap.readMultMRC = false;

		// *************************************************************************************
		cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
		// *************************************************************************************

		// ********************* Command line reading input with BOOST ************************

		try {
			po::options_description desc("Command line inputs");
			desc.add_options()
			("Inputfile", po::value<std::string>(), "(Mandatory) Name of input parameter file")
			("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file")
			("Particlesfile", po::value< std::string>(), "(Mandatory) Name of paricles file")
			("ReadPDB", "(Optional) If reading model file in PDB format")
			("ReadMRC", "(Optional) If reading particle file in MRC format")
			("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
			("DumpMaps", "(Optional) Dump maps after they were red from maps file")
			("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
			("help", "(Optional) Produce help message")
			;


			po::positional_options_description p;
			p.add("Inputfile", -1);
			p.add("Modelfile", -1);
			p.add("Particlesfile", -1);
			p.add("ReadPDB", -1);
			p.add("ReadMRC", -1);
			p.add("ReadMultipleMRC", -1);
			p.add("DumpMaps", -1);
			p.add("LoadMapDump", -1);

			po::variables_map vm;
			po::store(po::command_line_parser(ac, av).
					  options(desc).positional(p).run(), vm);
			po::notify(vm);

			if((ac < 6)) {
				std::cout << desc << std::endl;
				return 1;
			}
			if (vm.count("help")) {
				cout << "Usage: options_description [options]\n";
				cout << desc;
				return 1;
			}

			if (vm.count("Inputfile"))
			{
				cout << "Input file is: ";
				cout << vm["Inputfile"].as< std::string >() << "\n";
				infile = vm["Inputfile"].as< std::string >();
			}
			if (vm.count("Modelfile"))
			{
				cout << "Model file is: "
					 << vm["Modelfile"].as<  std::string  >() << "\n";
				modelfile = vm["Modelfile"].as<  std::string  >();
			}

			if (vm.count("ReadPDB"))
			{
				cout << "Reading model file in PDB format.\n";
				Model.readPDB = true;
			}

			if (vm.count("ReadMRC"))
			{
				cout << "Reading particle file in MRC format.\n";
				RefMap.readMRC=true;
			}

			if (vm.count("ReadMultipleMRC"))
			{
				cout << "Reading Multiple MRCs.\n";
				RefMap.readMultMRC=true;
			}

			if (vm.count("DumpMaps"))
			{
				cout << "Dumping Maps after reading from file.\n";
				param.dumpMap = true;
			}

			if (vm.count("LoadMapDump"))
			{
				cout << "Loading Map dump.\n";
				param.loadMap = true;
			}

			if (vm.count("Particlesfile"))
			{
				cout << "Paricle file is: "
					 << vm["Particlesfile"].as< std::string >() << "\n";
				mapfile = vm["Particlesfile"].as< std::string >();
			}
David Rohr's avatar
David Rohr committed
206
		}
David Rohr's avatar
David Rohr committed
207
		catch(std::exception& e)
208
		{
David Rohr's avatar
David Rohr committed
209
210
			cout << e.what() << "\n";
			return 1;
211
		}
David Rohr's avatar
David Rohr committed
212
213
214
215
216
217
218
219
220

		//check for consitency in multiple MRCs
		if(RefMap.readMultMRC && not(RefMap.readMRC))
		{
			cout << "For Multiple MRCs command --ReadMRC is necesary too";
			exit(1);
		}

		if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
David Rohr's avatar
David Rohr committed
221
		// ********************* Reading Parameter Input ***************************
222
		param.readParameters(infile.c_str());
David Rohr's avatar
David Rohr committed
223
224

		// ********************* Reading Model Input ******************************
225
		Model.readModel(modelfile.c_str());
David Rohr's avatar
David Rohr committed
226
227

		// ********************* Reading Particle Maps Input **********************
228
		RefMap.readRefMaps(param, mapfile.c_str());
David Rohr's avatar
David Rohr committed
229
		if (DebugOutput >= 2 && mpi_rank == 0) printf("Reading Input Data %f\n", timer.GetCurrentElapsedTime());
David Rohr's avatar
David Rohr committed
230
	}
231

David Rohr's avatar
David Rohr committed
232
#ifdef WITH_MPI
233
234
235
236
237
238
239
240
241
242
243
244
245
	if (mpi_size > 1)
	{
		if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
		MPI_Bcast(&param, sizeof(param), MPI_BYTE, 0, MPI_COMM_WORLD);
		//refCtf, CtfParam, angles automatically filled by precalculare function below

		MPI_Bcast(&Model, sizeof(Model), MPI_BYTE, 0, MPI_COMM_WORLD);
		if (mpi_rank != 0) Model.points = (bioem_model::bioem_model_point*) mallocchk(sizeof(bioem_model::bioem_model_point) * Model.nPointsModel);
		MPI_Bcast(Model.points, sizeof(bioem_model::bioem_model_point) * Model.nPointsModel, MPI_BYTE, 0, MPI_COMM_WORLD);

		MPI_Bcast(&RefMap, sizeof(RefMap), MPI_BYTE, 0, MPI_COMM_WORLD);
		if (mpi_rank != 0) RefMap.maps = (myfloat_t*) mallocchk(RefMap.refMapSize * sizeof(myfloat_t) * RefMap.ntotRefMap);
		MPI_Bcast(RefMap.maps, RefMap.refMapSize * sizeof(myfloat_t) * RefMap.ntotRefMap, MPI_BYTE, 0, MPI_COMM_WORLD);
David Rohr's avatar
David Rohr committed
246
		if (DebugOutput >= 2 && mpi_rank == 0) printf("MPI Broadcast of Input Data %f\n", timer.GetCurrentElapsedTime());
247
248
	}
	#endif
249

David Rohr's avatar
David Rohr committed
250
	// ****************** Precalculating Necessary Stuff *********************
251
252

	if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
253
	param.PrepareFFTs();
254
255
256
257
258
	if (DebugOutput >= 2 && mpi_rank == 0)
	{
		printf("Time Prepare FFTs %f\n", timer.GetCurrentElapsedTime());
		timer.ResetStart();
	}
259
	precalculate();
David Rohr's avatar
David Rohr committed
260

261
262
	if (getenv("BIOEM_DEBUG_BREAK"))
	{
263
264
265
		const int cut = atoi(getenv("BIOEM_DEBUG_BREAK"));
		if (param.nTotGridAngles > cut) param.nTotGridAngles = cut;
		if (param.nTotCTFs > cut) param.nTotCTFs = cut;
266
	}
David Rohr's avatar
David Rohr committed
267

268
269
270
271
272
	if (DebugOutput >= 2 && mpi_rank == 0)
	{
		printf("Time Precalculate %f\n", timer.GetCurrentElapsedTime());
		timer.ResetStart();
	}
273
	pProb.init(RefMap.ntotRefMap, param.nTotGridAngles, *this);
274

275
276
277
278
279
	if (DebugOutput >= 2 && mpi_rank == 0)
	{
		printf("Time Init Probabilities %f\n", timer.GetCurrentElapsedTime());
		timer.ResetStart();
	}
280
	deviceInit();
281
	if (DebugOutput >= 2 && mpi_rank == 0) printf("Time Device Init %f\n", timer.GetCurrentElapsedTime());
282

283
	return(0);
284
285
}

286
287
288
void bioem::cleanup()
{
	//Deleting allocated pointers
289
	free_device_host(pProb.ptr);
290
291
292
	RefMap.freePointers();
}

293
294
int bioem::precalculate()
{
David Rohr's avatar
David Rohr committed
295
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
296
	// **Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels**
David Rohr's avatar
David Rohr committed
297
	// **************************************************************************************
298
	HighResTimer timer;
299

300
	// Generating Grids of orientations
301
	if (DebugOutput >= 3) timer.ResetStart();
302
	param.CalculateGridsParam();
303

304
305
306
307
308
	if (DebugOutput >= 3)
	{
		printf("\tTime Precalculate Grids Param: %f\n", timer.GetCurrentElapsedTime());
		timer.ResetStart();
	}
309
310
	// Precalculating CTF Kernels stored in class Param
	param.CalculateRefCTF();
311

312
313
314
315
316
	if (DebugOutput >= 3)
	{
		printf("\tTime Precalculate CTFs: %f\n", timer.GetCurrentElapsedTime());
		timer.ResetStart();
	}
317
318
	//Precalculate Maps
	RefMap.precalculate(param, *this);
319
	if (DebugOutput >= 3) printf("\tTime Precalculate Maps: %f\n", timer.GetCurrentElapsedTime());
320

321
	return(0);
322
323
324
325
}

int bioem::run()
{
David Rohr's avatar
David Rohr committed
326
327
328
	// **************************************************************************************
	// **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
329

330
331
	if (param.printModel)
	{
David Rohr's avatar
David Rohr committed
332

333
334
335
		//....
		return(0);
	}
336

David Rohr's avatar
David Rohr committed
337
338
	// **** If we want to control the number of threads -> omp_set_num_threads(XX); ******
	// ****************** Declarying class of Probability Pointer  *************************
339

David Rohr's avatar
David Rohr committed
340
	if (mpi_rank == 0) printf("\tInitializing Probabilities\n");
341
342
343
	// Inizialzing Probabilites to zero and constant to -Infinity
	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
344
345
346
347
		bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

		pProbMap.Total = 0.0;
		pProbMap.Constoadd = -9999999;
348
		if (param.param_device.writeAngles)
349
		{
350
351
352
			for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
			{
				bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
353

354
355
356
				pProbAngle.forAngles = 0.0;
				pProbAngle.ConstAngle = -99999999;
			}
357
358
		}
	}
David Rohr's avatar
David Rohr committed
359
	// **************************************************************************************
360
361
	deviceStartRun();

David Rohr's avatar
David Rohr committed
362
	// ******************************** MAIN CYCLE ******************************************
David Rohr's avatar
David Rohr committed
363

David Rohr's avatar
David Rohr committed
364
	// *** Declaring Private variables for each thread *****
365
	mycomplex_t* proj_mapFFT;
366
	myfloat_t* conv_map = new myfloat_t[param.param_device.NumberPixels * param.param_device.NumberPixels];
367
	mycomplex_t* conv_mapFFT;
368
	myfloat_t sumCONV, sumsquareCONV;
369
370

	//allocating fftw_complex vector
371
372
	proj_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
	conv_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
373

David Rohr's avatar
David Rohr committed
374
	HighResTimer timer, timer2;
375

376
	if (DebugOutput >= 1 && mpi_rank == 0) printf("\tMain Loop GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)², Pixels %d², OMP Threads %d, MPI Ranks %d\n", param.nTotGridAngles, param.nTotCTFs, RefMap.ntotRefMap, 2 * param.param_device.maxDisplaceCenter + param.param_device.GridSpaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, omp_get_max_threads(), mpi_size);
David Rohr's avatar
David Rohr committed
377
378
379
380

	const int iOrientStart = (int) ((long long int) mpi_rank * param.nTotGridAngles / mpi_size);
	int iOrientEnd = (int) ((long long int) (mpi_rank + 1) * param.nTotGridAngles / mpi_size);
	if (iOrientEnd > param.nTotGridAngles) iOrientEnd = param.nTotGridAngles;
David Rohr's avatar
David Rohr committed
381

David Rohr's avatar
David Rohr committed
382
	for (int iOrient = iOrientStart; iOrient < iOrientEnd; iOrient++)
383
	{
David Rohr's avatar
David Rohr committed
384
385
		// ***************************************************************************************
		// ***** Creating Projection for given orientation and transforming to Fourier space *****
David Rohr's avatar
David Rohr committed
386
387
		if (DebugOutput >= 1) timer2.ResetStart();
		if (DebugOutput >= 2) timer.ResetStart();
388
		createProjection(iOrient, proj_mapFFT);
David Rohr's avatar
David Rohr committed
389
		if (DebugOutput >= 2) printf("Time Projection %d: %f (rank %d)\n", iOrient, timer.GetCurrentElapsedTime(), mpi_rank);
390

David Rohr's avatar
David Rohr committed
391
392
		// ***************************************************************************************
		// ***** **** Internal Loop over convolutions **** *****
393
394
		for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
		{
David Rohr's avatar
David Rohr committed
395
			// *** Calculating convolutions of projection map and crosscorrelations ***
396

397
			if (DebugOutput >= 2) timer.ResetStart();
398
			createConvolutedProjectionMap(iOrient, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
David Rohr's avatar
David Rohr committed
399
			if (DebugOutput >= 2) printf("Time Convolution %d %d: %f (rank %d)\n", iOrient, iConv, timer.GetCurrentElapsedTime(), mpi_rank);
400

David Rohr's avatar
David Rohr committed
401
402
			// ***************************************************************************************
			// *** Comparing each calculated convoluted map with all experimental maps ***
403
			if (DebugOutput >= 2) timer.ResetStart();
404
			compareRefMaps(iOrient, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
405

406
407
408
409
410
411
412
413
414
415
			if (DebugOutput >= 2)
			{
				const double compTime = timer.GetCurrentElapsedTime();
				const int nShifts = 2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1;
				const double nFlops = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
									  (((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 5. + 25.) / compTime;
				const double nGBs = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
									(((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 2. + 8.) * (double) sizeof(myfloat_t) / compTime;
				const double nGBs2 = (double) RefMap.ntotRefMap * ((double) param.param_device.NumberPixels * (double) param.param_device.NumberPixels + 8.) * (double) sizeof(myfloat_t) / compTime;

David Rohr's avatar
David Rohr committed
416
				printf("Time Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s) (rank %d)\n", iOrient, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000., mpi_rank);
417
			}
418
		}
David Rohr's avatar
David Rohr committed
419
		if (DebugOutput >= 1) printf("Total time for orientation %d: %f (rank %d)\n", iOrient, timer2.GetCurrentElapsedTime(), mpi_rank);
420
421
422
423
	}
	//deallocating fftw_complex vector
	myfftw_free(proj_mapFFT);
	myfftw_free(conv_mapFFT);
424
	delete[] conv_map;
David Rohr's avatar
David Rohr committed
425

426
427
	deviceFinishRun();

David Rohr's avatar
David Rohr committed
428
	// ************* Writing Out Probabilities ***************
429

David Rohr's avatar
David Rohr committed
430
	// *** Angular Probability ***
David Rohr's avatar
David Rohr committed
431

David Rohr's avatar
David Rohr committed
432
#ifdef WITH_MPI
433
	if (mpi_size > 1)
David Rohr's avatar
David Rohr committed
434
	{
435
436
		if (DebugOutput >= 1 && mpi_rank == 0) timer.ResetStart();
		//Reduce Constant and summarize probabilities
David Rohr's avatar
David Rohr committed
437
		{
438
439
440
			myfloat_t* tmp1 = new myfloat_t[RefMap.ntotRefMap];
			myfloat_t* tmp2 = new myfloat_t[RefMap.ntotRefMap];
			myfloat_t* tmp3 = new myfloat_t[RefMap.ntotRefMap];
David Rohr's avatar
David Rohr committed
441
442
			for (int i = 0;i < RefMap.ntotRefMap;i++)
			{
443
					tmp1[i] = pProb.getProbMap(i).Constoadd;
David Rohr's avatar
David Rohr committed
444
			}
445
			MPI_Allreduce(tmp1, tmp2, RefMap.ntotRefMap, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
David Rohr's avatar
David Rohr committed
446
447
			for (int i = 0;i < RefMap.ntotRefMap;i++)
			{
448
449
450
451
452
453
454
455
456
457
458
				bioem_Probability_map& pProbMap = pProb.getProbMap(i);
				tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
			}
			MPI_Reduce(tmp1, tmp3, RefMap.ntotRefMap, MY_MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD);

			//Find MaxProb
			MPI_Status mpistatus;
			{
				int* tmpi1 = new int[RefMap.ntotRefMap];
				int* tmpi2 = new int[RefMap.ntotRefMap];
				for (int i = 0;i < RefMap.ntotRefMap;i++)
David Rohr's avatar
David Rohr committed
459
				{
460
461
					bioem_Probability_map& pProbMap = pProb.getProbMap(i);
					tmpi1[i] = tmp2[i] <= pProbMap.Constoadd ? mpi_rank : -1;
David Rohr's avatar
David Rohr committed
462
				}
463
464
				MPI_Allreduce(tmpi1, tmpi2, RefMap.ntotRefMap, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
				for (int i = 0;i < RefMap.ntotRefMap;i++)
David Rohr's avatar
David Rohr committed
465
				{
466
					if (tmpi2[i] == -1)
David Rohr's avatar
David Rohr committed
467
					{
468
						if (mpi_rank == 0) printf("Error: Could not find highest probability\n");
David Rohr's avatar
David Rohr committed
469
					}
470
					else if (tmpi2[i] != 0) //Skip if rank 0 already has highest probability
David Rohr's avatar
David Rohr committed
471
					{
472
473
474
475
476
477
478
479
						if (mpi_rank == 0)
						{
							MPI_Recv(&pProb.getProbMap(i).max, sizeof(pProb.getProbMap(i).max), MPI_BYTE, tmpi2[i], i, MPI_COMM_WORLD, &mpistatus);
						}
						else if (mpi_rank == tmpi2[i])
						{
							MPI_Send(&pProb.getProbMap(i).max, sizeof(pProb.getProbMap(i).max), MPI_BYTE, 0, i, MPI_COMM_WORLD);
						}
David Rohr's avatar
David Rohr committed
480
481
					}
				}
482
483
				delete[] tmpi1;
				delete[] tmpi2;
David Rohr's avatar
David Rohr committed
484
			}
David Rohr's avatar
David Rohr committed
485

486
			if (mpi_rank == 0)
David Rohr's avatar
David Rohr committed
487
			{
488
489
490
491
492
493
				for (int i = 0;i < RefMap.ntotRefMap;i++)
				{
						bioem_Probability_map& pProbMap = pProb.getProbMap(i);
						pProbMap.Total = tmp3[i];
						pProbMap.Constoadd = tmp2[i];
				}
David Rohr's avatar
David Rohr committed
494
			}
David Rohr's avatar
David Rohr committed
495

496
497
498
499
			delete[] tmp1;
			delete[] tmp2;
			delete[] tmp3;
			if (DebugOutput >= 1 && mpi_rank == 0 && mpi_size > 1) printf("Time MPI Reduction: %f\n", timer.GetCurrentElapsedTime());
David Rohr's avatar
David Rohr committed
500
		}
501
502
503

		//Angle Reduction and Probability summation for individual angles
		if (param.param_device.writeAngles)
David Rohr's avatar
David Rohr committed
504
		{
505
506
507
508
509
			const int count = RefMap.ntotRefMap * param.nTotGridAngles;
			myfloat_t* tmp1 = new myfloat_t[count];
			myfloat_t* tmp2 = new myfloat_t[count];
			myfloat_t* tmp3 = new myfloat_t[count];
			for (int i = 0;i < RefMap.ntotRefMap;i++)
David Rohr's avatar
David Rohr committed
510
			{
511
					tmp1[i] = pProb.getProbMap(i).Constoadd;
David Rohr's avatar
David Rohr committed
512
			}
513
			MPI_Allreduce(tmp1, tmp2, count, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
David Rohr's avatar
David Rohr committed
514
515
516
517
518
			for (int i = 0;i < RefMap.ntotRefMap;i++)
			{
				for (int j = 0;j < param.nTotGridAngles;j++)
				{
					bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
519
520
521
522
523
524
525
526
527
528
529
530
531
532
					tmp1[i * param.nTotGridAngles + j] = pProbAngle.forAngles * exp(pProbAngle.ConstAngle - tmp2[i * param.nTotGridAngles + j]);
				}
			}
			MPI_Reduce(tmp1, tmp3, count, MY_MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD);
			if (mpi_rank == 0)
			{
				for (int i = 0;i < RefMap.ntotRefMap;i++)
				{
					for (int j = 0;j < param.nTotGridAngles;j++)
					{
						bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
						pProbAngle.forAngles = tmp3[i * param.nTotGridAngles + j];
						pProbAngle.ConstAngle = tmp2[i * param.nTotGridAngles + j];
					}
David Rohr's avatar
David Rohr committed
533
534
				}
			}
535
536
537
			delete[] tmp1;
			delete[] tmp2;
			delete[] tmp3;
David Rohr's avatar
David Rohr committed
538
		}
539
	}
David Rohr's avatar
David Rohr committed
540
#endif
541

David Rohr's avatar
David Rohr committed
542
	if (mpi_rank == 0)
543
	{
David Rohr's avatar
David Rohr committed
544
545
546
547
548
549
550
551
552
553
554
555
		ofstream angProbfile;
		if(param.param_device.writeAngles)
		{
			angProbfile.open ("ANG_PROB");
		}

		ofstream outputProbFile;
		outputProbFile.open ("Output_Probabilities");
		for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
		{
			// **** Total Probability ***
			bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);
556

David Rohr's avatar
David Rohr committed
557
			outputProbFile << "RefMap " << iRefMap << " Probability  "  << log(pProbMap.Total) + pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProbMap.Constoadd  << "\n";
558

David Rohr's avatar
David Rohr committed
559
			outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
560

David Rohr's avatar
David Rohr committed
561
562
563
564
565
566
567
568
569
570
571
			// *** Param that maximize probability****
			outputProbFile << (pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5) * (log(2 * M_PI) + 1) + log(param.param_device.volu)) << " ";
			outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[0] << " ";
			outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[1] << " ";
			outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[2] << " ";
			outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[0] << " ";
			outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << " ";
			outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " ";
			outputProbFile << pProbMap.max.max_prob_cent_x << " ";
			outputProbFile << pProbMap.max.max_prob_cent_y;
			outputProbFile << "\n";
572

David Rohr's avatar
David Rohr committed
573
			// *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap);
574

David Rohr's avatar
David Rohr committed
575
			if(param.param_device.writeAngles)
576
			{
David Rohr's avatar
David Rohr committed
577
578
579
				for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
				{
					bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
580

David Rohr's avatar
David Rohr committed
581
582
					angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << log(pProbAngle.forAngles) + pProbAngle.ConstAngle + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
				}
583
584
			}
		}
585

David Rohr's avatar
David Rohr committed
586
587
588
589
590
		if(param.param_device.writeAngles)
		{
			angProbfile.close();
		}
		outputProbFile.close();
591
	}
592

593
	return(0);
594
595
}

596
int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
597
{
David Rohr's avatar
David Rohr committed
598
599
	//***************************************************************************************
	//***** BioEM routine for comparing reference maps to convoluted maps *****
600
	if (FFTAlgo)
601
	{
David Rohr's avatar
David Rohr committed
602
		//With FFT Algorithm
603
604
		#pragma omp parallel for
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
605
		{
606
			const int num = omp_get_thread_num();
607
			calculateCCFFT(iRefMap, iOrient, iConv, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
608
609
610
		}
	}
	else
611
	{
David Rohr's avatar
David Rohr committed
612
		//Without FFT Algorithm
613
		#pragma omp parallel for
614
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
615
		{
616
			compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, conv_map, pProb, param.param_device, RefMap);
617
618
619
620
621
		}
	}
	return(0);
}

622
inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC)
623
{
David Rohr's avatar
David Rohr committed
624
625
	//***************************************************************************************
	//***** Calculating cross correlation in FFTALGOrithm *****
Pilar Cossio's avatar
Pilar Cossio committed
626

627
	const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize];
628
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
629
	{
630
631
		localCCT[i][0] = localConvFFT[i][0] * RefMapFFT[i][0] + localConvFFT[i][1] * RefMapFFT[i][1];
		localCCT[i][1] = localConvFFT[i][1] * RefMapFFT[i][0] - localConvFFT[i][0] * RefMapFFT[i][1];
632
633
	}

634
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
635

636
	doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
637
}
638

639
int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
640
{
David Rohr's avatar
David Rohr committed
641
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
642
643
	// ****  BioEM Create Projection routine in Euler angle predefined grid******************
	// ********************* and turns projection into Fourier space ************************
David Rohr's avatar
David Rohr committed
644
	// **************************************************************************************
645

646
647
	cuda_custom_timeslot("Projection", 0);

648
649
	myfloat3_t RotatedPointsModel[Model.nPointsModel];
	myfloat_t rotmat[3][3];
650
	myfloat_t alpha, gam, beta;
651
	myfloat_t* localproj;
652

653
	localproj = param.fft_scratch_real[omp_get_thread_num()];
654
	memset(localproj, 0, param.param_device.NumberPixels * param.param_device.NumberPixels * sizeof(*localproj));
655

656
657
658
	alpha = param.angles[iMap].pos[0];
	beta = param.angles[iMap].pos[1];
	gam = param.angles[iMap].pos[2];
659

David Rohr's avatar
David Rohr committed
660
	// **** To see how things are going: cout << "Id " << omp_get_thread_num() <<  " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***
661

David Rohr's avatar
David Rohr committed
662
	// ********** Creat Rotation with pre-defiend grid of orientations**********
663
664
665
666
667
668
669
670
671
672
673
	rotmat[0][0] = cos(gam) * cos(alpha) - cos(beta) * sin(alpha) * sin(gam);
	rotmat[0][1] = cos(gam) * sin(alpha) + cos(beta) * cos(alpha) * sin(gam);
	rotmat[0][2] = sin(gam) * sin(beta);
	rotmat[1][0] = -sin(gam) * cos(alpha) - cos(beta) * sin(alpha) * cos(gam);
	rotmat[1][1] = -sin(gam) * sin(alpha) + cos(beta) * cos(alpha) * cos(gam);
	rotmat[1][2] = cos(gam) * sin(beta);
	rotmat[2][0] = sin(beta) * sin(alpha);
	rotmat[2][1] = -sin(beta) * cos(alpha);
	rotmat[2][2] = cos(beta);

	for(int n = 0; n < Model.nPointsModel; n++)
674
	{
675
676
677
		RotatedPointsModel[n].pos[0] = 0.0;
		RotatedPointsModel[n].pos[1] = 0.0;
		RotatedPointsModel[n].pos[2] = 0.0;
678
	}
679
	for(int n = 0; n < Model.nPointsModel; n++)
680
	{
681
		for(int k = 0; k < 3; k++)
682
		{
683
			for(int j = 0; j < 3; j++)
684
			{
685
				RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.points[n].point.pos[j];
686
687
688
689
690
691
			}
		}
	}

	int i, j;

David Rohr's avatar
David Rohr committed
692
	// ************ Projection over the Z axis********************
693
	for(int n = 0; n < Model.nPointsModel; n++)
694
695
	{
		//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
696
697
		i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
		j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
698

699
700
701
702
703
704
		if (i < 0 || j < 0 || i >= param.param_device.NumberPixels || j >= param.param_device.NumberPixels)
		{
			if (DebugOutput >= 3) cout << "Model Point out of map: " << i << ", " << j << "\n";
			continue;
		}

705
		localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density / Model.NormDen;
706
707
	}

David Rohr's avatar
David Rohr committed
708
	// **** Output Just to check****
709
#ifdef PILAR_DEBUG
710
	if(iMap == 10)
711
712
713
714
715
716
	{
		ofstream myexamplemap;
		ofstream myexampleRot;
		myexamplemap.open ("MAP_i10");
		myexampleRot.open ("Rot_i10");
		myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
717
		for(int k = 0; k < param.param_device.NumberPixels; k++)
718
		{
719
			for(int j = 0; j < param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j << " " << localproj[k * param.param_device.NumberPixels + j];
720
721
		}
		myexamplemap << " \n";
722
		for(int n = 0; n < Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
723
724
725
		myexamplemap.close();
		myexampleRot.close();
	}
726
#endif
727

David Rohr's avatar
David Rohr committed
728
729
	// ***** Converting projection to Fourier Space for Convolution later with kernel****
	// ********** Omp Critical is necessary with FFTW*******
730
	myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localproj, mapFFT);
731

732
733
	cuda_custom_timeslot_end;

734
735
736
	return(0);
}

737
int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
738
{
David Rohr's avatar
David Rohr committed
739
740
	// **************************************************************************************
	// ****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
David Rohr's avatar
David Rohr committed
741
742
	// **************** calculated Projection with convoluted precalculated Kernel***********
	// *************** and Backtransforming it to real Space ********************************
David Rohr's avatar
David Rohr committed
743
	// **************************************************************************************
744

745
746
	cuda_custom_timeslot("Convolution", 1);

747
	mycomplex_t* tmp = param.fft_scratch_complex[omp_get_thread_num()];
748

David Rohr's avatar
David Rohr committed
749
	// **** Multiplying FFTmap with corresponding kernel ****
750
	const mycomplex_t* refCTF = &param.refCTF[iConv * param.FFTMapSize];
751
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
752
	{
753
754
755
		localmultFFT[i][0] = lproj[i][0] * refCTF[i][0] + lproj[i][1] * refCTF[i][1];
		localmultFFT[i][1] = lproj[i][1] * refCTF[i][0] - lproj[i][0] * refCTF[i][1];
		// cout << "GG " << i << " " << j << " " << refCTF[i][0] << " " << refCTF[i][1] <<" " <<lproj[i][0] <<" " <<lproj[i][1] << "\n";
756
757
	}

David Rohr's avatar
David Rohr committed
758
	// *** Calculating Cross-correlations of cal-convoluted map with its self *****
759
	sumC = localmultFFT[0][0];
David Rohr's avatar
David Rohr committed
760

761
	sumsquareC = 0;
762
	if (FFTAlgo)
763
	{
764
765
		int jloopend = param.param_device.NumberFFTPixels1D;
		if ((param.param_device.NumberPixels & 1) == 0) jloopend--;
766
767
		for(int i = 0; i < param.param_device.NumberPixels; i++)
		{
768
			for (int j = 1;j < jloopend;j++)
769
770
771
772
773
774
			{
				int k = i * param.param_device.NumberFFTPixels1D + j;
				sumsquareC += (localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1]) * 2;
			}
			int k = i * param.param_device.NumberFFTPixels1D;
			sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1];
775
776
			if ((param.param_device.NumberPixels & 1) == 0)
			{
777
				k += param.param_device.NumberFFTPixels1D - 1;
778
779
				sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1];
			}
780
		}
David Rohr's avatar
David Rohr committed
781

782
783
		myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
		sumsquareC = sumsquareC / norm2;
784
	}
785
786
787
788
789
790
791
	else
	{
		//FFTW_C2R will destroy the input array, so we have to work on a copy here
		memcpy(tmp, localmultFFT, sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);

		// **** Bringing convoluted Map to real Space ****
		myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, Mapconv);
792

793
794
795
796
797
798
799
800
801
		for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
		{
			sumsquareC += Mapconv[i] * Mapconv[i];
		}

		myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
		myfloat_t norm4 = norm2 * norm2;
		sumsquareC = sumsquareC / norm4;
	}
David Rohr's avatar
David Rohr committed
802

803
804
	cuda_custom_timeslot_end;

805
	return(0);
806
807
}

808
int bioem::calcross_cor(myfloat_t* localmap, myfloat_t& sum, myfloat_t& sumsquare)
809
{
David Rohr's avatar
David Rohr committed
810
	// *********************** Routine to calculate Cross correlations***********************
811

812
813
	sum = 0.0;
	sumsquare = 0.0;
814
815
816
817
818
	for (int i = 0; i < param.param_device.NumberPixels; i++)
	{
		for (int j = 0; j < param.param_device.NumberPixels; j++)
		{
			// Calculate Sum of pixels
819
			sum += localmap[i * param.param_device.NumberPixels + j];
820
			// Calculate Sum of pixels squared
821
			sumsquare += localmap[i * param.param_device.NumberPixels + j] * localmap[i * param.param_device.NumberPixels + j];
822
823
824
		}
	}
	return(0);
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
}

int bioem::deviceInit()
{
	return(0);
}

int bioem::deviceStartRun()
{
	return(0);
}

int bioem::deviceFinishRun()
{
	return(0);
}
841
842
843
844
845
846
847
848
849
850

void* bioem::malloc_device_host(size_t size)
{
	return(mallocchk(size));
}

void bioem::free_device_host(void* ptr)
{
	free(ptr);
}