bioem.cpp 32.6 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
	nProjectionsAtOnce = getenv("BIOEM_PROJECTIONS_AT_ONCE") == NULL ? 1 : atoi(getenv("BIOEM_PROJECTIONS_AT_ONCE"));
85
86
87
88
89
90
91
92
}

bioem::~bioem()
{
}

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

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

David Rohr's avatar
David Rohr committed
101
102
103
104
105
	if (mpi_rank == 0)
	{
		// *** Inizialzing default variables ***
		std::string infile, modelfile, mapfile;
		Model.readPDB = false;
106
		param.param_device.writeAngles = false;
David Rohr's avatar
David Rohr committed
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
206
		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
207
		}
David Rohr's avatar
David Rohr committed
208
		catch(std::exception& e)
209
		{
David Rohr's avatar
David Rohr committed
210
211
			cout << e.what() << "\n";
			return 1;
212
		}
David Rohr's avatar
David Rohr committed
213
214
215
216
217
218
219
220
221

		//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
222
		// ********************* Reading Parameter Input ***************************
223
		param.readParameters(infile.c_str());
David Rohr's avatar
David Rohr committed
224
225

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

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

David Rohr's avatar
David Rohr committed
233
#ifdef WITH_MPI
234
235
236
237
238
239
240
241
242
243
244
245
246
	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
247
		if (DebugOutput >= 2 && mpi_rank == 0) printf("MPI Broadcast of Input Data %f\n", timer.GetCurrentElapsedTime());
248
249
	}
	#endif
250

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

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

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

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

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

284
	return(0);
285
286
}

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

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

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

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

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

322
	return(0);
323
324
325
326
}

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

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

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

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

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

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

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

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

365
	mycomplex_t* proj_mapsFFT;
366
	myfloat_t* conv_map = NULL;
367
	mycomplex_t* conv_mapFFT;
368
	myfloat_t sumCONV, sumsquareCONV;
369
370

	//allocating fftw_complex vector
371
372
	const int ProjMapSize = (param.FFTMapSize + 64) & ~63;	//Make sure this is properly aligned for fftw..., Actually this should be ensureb by using FFTMapSize, but it is not due to a bug in CUFFT which cannot handle padding properly
	proj_mapsFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * ProjMapSize * nProjectionsAtOnce);
373
	conv_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
374
	if (!FFTAlgo) conv_map = (myfloat_t*) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
375

David Rohr's avatar
David Rohr committed
376
	HighResTimer timer, timer2;
377

378
	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
379
380
381
382

	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
383

384
	for (int iOrientAtOnce = iOrientStart; iOrientAtOnce < iOrientEnd; iOrientAtOnce += nProjectionsAtOnce)
385
	{
David Rohr's avatar
David Rohr committed
386
387
		// ***************************************************************************************
		// ***** Creating Projection for given orientation and transforming to Fourier space *****
David Rohr's avatar
David Rohr committed
388
389
		if (DebugOutput >= 1) timer2.ResetStart();
		if (DebugOutput >= 2) timer.ResetStart();
390
391
392
393
394
395
396
		int iTmpEnd = std::min(iOrientEnd, iOrientAtOnce + nProjectionsAtOnce);
#pragma omp parallel for
		for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
		{
			createProjection(iOrient, &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize]);
		}
		if (DebugOutput >= 2) printf("\tTime Projection %d: %f (rank %d)\n", iOrientAtOnce, timer.GetCurrentElapsedTime(), mpi_rank);
397

398
		for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
399
		{
400
401
402
403
404
405
			mycomplex_t* proj_mapFFT = &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize];
			// ***************************************************************************************
			// ***** **** Internal Loop over convolutions **** *****
			for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
			{
				// *** Calculating convolutions of projection map and crosscorrelations ***
406

407
408
409
				if (DebugOutput >= 2) timer.ResetStart();
				createConvolutedProjectionMap(iOrient, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
				if (DebugOutput >= 2) printf("\t\tTime Convolution %d %d: %f (rank %d)\n", iOrient, iConv, timer.GetCurrentElapsedTime(), mpi_rank);
410

411
412
413
414
				// ***************************************************************************************
				// *** Comparing each calculated convoluted map with all experimental maps ***
				if (DebugOutput >= 2) timer.ResetStart();
				compareRefMaps(iOrient, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
415

416
417
418
419
420
421
422
423
424
425
426
427
428
429
				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;

					printf("\t\tTime 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);
				}
			}
			if (DebugOutput >= 1)
430
			{
431
432
				printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient, timer2.GetCurrentElapsedTime(), mpi_rank);
				timer2.ResetStart();
433
			}
434
435
436
		}
	}
	//deallocating fftw_complex vector
437
	myfftw_free(proj_mapsFFT);
438
	myfftw_free(conv_mapFFT);
439
	if (!FFTAlgo) myfftw_free(conv_map);
David Rohr's avatar
David Rohr committed
440

441
442
	deviceFinishRun();

David Rohr's avatar
David Rohr committed
443
	// ************* Writing Out Probabilities ***************
444

David Rohr's avatar
David Rohr committed
445
	// *** Angular Probability ***
David Rohr's avatar
David Rohr committed
446

David Rohr's avatar
David Rohr committed
447
#ifdef WITH_MPI
448
	if (mpi_size > 1)
David Rohr's avatar
David Rohr committed
449
	{
450
451
		if (DebugOutput >= 1 && mpi_rank == 0) timer.ResetStart();
		//Reduce Constant and summarize probabilities
David Rohr's avatar
David Rohr committed
452
		{
453
454
455
			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
456
457
			for (int i = 0;i < RefMap.ntotRefMap;i++)
			{
458
					tmp1[i] = pProb.getProbMap(i).Constoadd;
David Rohr's avatar
David Rohr committed
459
			}
460
			MPI_Allreduce(tmp1, tmp2, RefMap.ntotRefMap, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
David Rohr's avatar
David Rohr committed
461
462
			for (int i = 0;i < RefMap.ntotRefMap;i++)
			{
463
464
465
466
467
468
469
470
471
472
473
				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
474
				{
475
476
					bioem_Probability_map& pProbMap = pProb.getProbMap(i);
					tmpi1[i] = tmp2[i] <= pProbMap.Constoadd ? mpi_rank : -1;
David Rohr's avatar
David Rohr committed
477
				}
478
479
				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
480
				{
481
					if (tmpi2[i] == -1)
David Rohr's avatar
David Rohr committed
482
					{
483
						if (mpi_rank == 0) printf("Error: Could not find highest probability\n");
David Rohr's avatar
David Rohr committed
484
					}
485
					else if (tmpi2[i] != 0) //Skip if rank 0 already has highest probability
David Rohr's avatar
David Rohr committed
486
					{
487
488
489
490
491
492
493
494
						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
495
496
					}
				}
497
498
				delete[] tmpi1;
				delete[] tmpi2;
David Rohr's avatar
David Rohr committed
499
			}
David Rohr's avatar
David Rohr committed
500

501
			if (mpi_rank == 0)
David Rohr's avatar
David Rohr committed
502
			{
503
504
505
506
507
508
				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
509
			}
David Rohr's avatar
David Rohr committed
510

511
512
513
514
			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
515
		}
516
517
518

		//Angle Reduction and Probability summation for individual angles
		if (param.param_device.writeAngles)
David Rohr's avatar
David Rohr committed
519
		{
520
521
522
523
524
			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
525
			{
526
					tmp1[i] = pProb.getProbMap(i).Constoadd;
David Rohr's avatar
David Rohr committed
527
			}
528
			MPI_Allreduce(tmp1, tmp2, count, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
David Rohr's avatar
David Rohr committed
529
530
531
532
533
			for (int i = 0;i < RefMap.ntotRefMap;i++)
			{
				for (int j = 0;j < param.nTotGridAngles;j++)
				{
					bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
534
535
536
537
538
539
540
541
542
543
544
545
546
547
					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
548
549
				}
			}
550
551
552
			delete[] tmp1;
			delete[] tmp2;
			delete[] tmp3;
David Rohr's avatar
David Rohr committed
553
		}
554
	}
David Rohr's avatar
David Rohr committed
555
#endif
556

David Rohr's avatar
David Rohr committed
557
	if (mpi_rank == 0)
558
	{
David Rohr's avatar
David Rohr committed
559
560
561
562
563
		ofstream angProbfile;
		if(param.param_device.writeAngles)
		{
			angProbfile.open ("ANG_PROB");
		}
564
		ofstream ccProbfile;
565
566
567
568
		if(param.param_device.writeCC)
		{
			ccProbfile.open ("CROSS_CORRELATION");
		}
David Rohr's avatar
David Rohr committed
569
570
571
572
573
574
575

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

David Rohr's avatar
David Rohr committed
577
			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";
578

David Rohr's avatar
David Rohr committed
579
			outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
580

David Rohr's avatar
David Rohr committed
581
582
583
584
585
586
587
588
589
			// *** 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 << " ";
Pilar Cossio's avatar
Pilar Cossio committed
590
591
592
593
			outputProbFile << pProbMap.max.max_prob_cent_y << " " ;
	                outputProbFile << pProbMap.max.max_prob_norm  << " " ;
        	        outputProbFile << pProbMap.max.max_prob_mu ;
                	outputProbFile << "\n";
594

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

David Rohr's avatar
David Rohr committed
597
			if(param.param_device.writeAngles)
598
			{
David Rohr's avatar
David Rohr committed
599
600
601
				for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
				{
					bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
602

David Rohr's avatar
David Rohr committed
603
604
					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";
				}
605
			}
606
607
608
609
610
611
612
613
614
615
616
617
			if(param.param_device.writeCC)
			{
				int  cc=0;
				for (int cent_x = 0; cent_x < param.param_device.NumberPixels; cent_x = cent_x + param.param_device.CCdisplace)
				{
					for (int cent_y = 0; cent_y < param.param_device.NumberPixels; cent_y = cent_y + param.param_device.CCdisplace)
					{
						bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
						ccProbfile << " " << iRefMap << " " << cent_x << " " << cent_y << " " << log(pProbCC.forCC) <<"\n";
					}
				}

618
			}
619
		}
620

David Rohr's avatar
David Rohr committed
621
622
623
624
		if(param.param_device.writeAngles)
		{
			angProbfile.close();
		}
625

626
627
628
629
		if(param.param_device.writeCC)
		{
			ccProbfile.close();
		}
630

David Rohr's avatar
David Rohr committed
631
		outputProbFile.close();
632
	}
633

634
	return(0);
635
636
}

637
int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
638
{
David Rohr's avatar
David Rohr committed
639
640
	//***************************************************************************************
	//***** BioEM routine for comparing reference maps to convoluted maps *****
641
	if (FFTAlgo)
642
	{
David Rohr's avatar
David Rohr committed
643
		//With FFT Algorithm
644
		#pragma omp parallel for schedule(dynamic, 1)
645
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
646
		{
647
			const int num = omp_get_thread_num();
648
			calculateCCFFT(iRefMap, iOrient, iConv, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
649
650
651
		}
	}
	else
652
	{
David Rohr's avatar
David Rohr committed
653
		//Without FFT Algorithm
654
		#pragma omp parallel for schedule(dynamic, 1)
655
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
656
		{
657
			compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, conv_map, pProb, param.param_device, RefMap);
658
659
660
661
662
		}
	}
	return(0);
}

663
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)
664
{
David Rohr's avatar
David Rohr committed
665
666
	//***************************************************************************************
	//***** Calculating cross correlation in FFTALGOrithm *****
Pilar Cossio's avatar
Pilar Cossio committed
667

668
	const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize];
669
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
670
	{
671
672
		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];
673
674
	}

675
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
676

677
	doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
678
}
679

680
int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
681
{
David Rohr's avatar
David Rohr committed
682
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
683
684
	// ****  BioEM Create Projection routine in Euler angle predefined grid******************
	// ********************* and turns projection into Fourier space ************************
David Rohr's avatar
David Rohr committed
685
	// **************************************************************************************
686

687
688
	cuda_custom_timeslot("Projection", 0);

689
690
	myfloat3_t RotatedPointsModel[Model.nPointsModel];
	myfloat_t rotmat[3][3];
691
	myfloat_t alpha, gam, beta;
692
	myfloat_t* localproj;
693

694
	localproj = param.fft_scratch_real[omp_get_thread_num()];
695
	memset(localproj, 0, param.param_device.NumberPixels * param.param_device.NumberPixels * sizeof(*localproj));
696

697
698
699
	alpha = param.angles[iMap].pos[0];
	beta = param.angles[iMap].pos[1];
	gam = param.angles[iMap].pos[2];
700

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

David Rohr's avatar
David Rohr committed
703
	// ********** Creat Rotation with pre-defiend grid of orientations**********
704
705
706
707
708
709
710
711
712
713
714
	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++)
715
	{
716
717
718
		RotatedPointsModel[n].pos[0] = 0.0;
		RotatedPointsModel[n].pos[1] = 0.0;
		RotatedPointsModel[n].pos[2] = 0.0;
719
	}
720
	for(int n = 0; n < Model.nPointsModel; n++)
721
	{
722
		for(int k = 0; k < 3; k++)
723
		{
724
			for(int j = 0; j < 3; j++)
725
			{
726
				RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.points[n].point.pos[j];
727
728
729
730
731
732
			}
		}
	}

	int i, j;

Pilar Cossio's avatar
Pilar Cossio committed
733
	//********** Projection with radius ***************
734
735
736
737
	if(param.doaaradius)
	{
		int irad;
		myfloat_t dist, rad2;
Pilar Cossio's avatar
Pilar Cossio committed
738

739
740
741
742
743
		for(int n = 0; n < Model.nPointsModel; n++)
		{
			//Getting Centers of Sphere
			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);
Pilar Cossio's avatar
Pilar Cossio committed
744
745
			//Getting the radius
			irad=int( Model.points[n].radius / param.pixelSize );
746
			rad2= Model.points[n].radius * Model.points[n].radius;
Pilar Cossio's avatar
Pilar Cossio committed
747

748
749
			//Projecting over the radius
			for(int ii= i - irad; ii < i + irad; ii++)
Pilar Cossio's avatar
Pilar Cossio committed
750
751
752
			{	
				for(int jj = j - irad; jj < j + irad; jj++)
				{
753
754
					dist= ( (myfloat_t) (ii-i)*(ii-i)+(jj-j)*(jj-j) ) *  param.pixelSize ;
					if( dist < rad2 )
Pilar Cossio's avatar
Pilar Cossio committed
755
					{
756
757
						localproj[ii * param.param_device.NumberPixels + jj] += param.pixelSize * 2 * sqrt( rad2 - dist ) * Model.points[n].density
							/ Model.NormDen * 3 / (4 * M_PI * rad2 *  Model.points[n].radius); 	
Pilar Cossio's avatar
Pilar Cossio committed
758
					}
759
760
				}
			}
Pilar Cossio's avatar
Pilar Cossio committed
761
762
763
		}
	}
	else 
764
	{
765
		// ************ Simple Projection over the Z axis********************
Pilar Cossio's avatar
Pilar Cossio committed
766
767
		for(int n = 0; n < Model.nPointsModel; n++)
		{
768
769
770
			//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
			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);
771

772
773
774
775
776
			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;
			}
777

778
779
			localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density / Model.NormDen;
		}
780
781
	}

David Rohr's avatar
David Rohr committed
782
	// **** Output Just to check****
783
#ifdef PILAR_DEBUG
784
	if(iMap == 10)
785
786
787
788
789
790
	{
		ofstream myexamplemap;
		ofstream myexampleRot;
		myexamplemap.open ("MAP_i10");
		myexampleRot.open ("Rot_i10");
		myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
791
		for(int k = 0; k < param.param_device.NumberPixels; k++)
792
		{
793
			for(int j = 0; j < param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j << " " << localproj[k * param.param_device.NumberPixels + j];
794
795
		}
		myexamplemap << " \n";
796
		for(int n = 0; n < Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
797
798
799
		myexamplemap.close();
		myexampleRot.close();
	}
800
#endif
801

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

806
807
	cuda_custom_timeslot_end;

808
809
810
	return(0);
}

811
int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
812
{
David Rohr's avatar
David Rohr committed
813
814
	// **************************************************************************************
	// ****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
David Rohr's avatar
David Rohr committed
815
816
	// **************** calculated Projection with convoluted precalculated Kernel***********
	// *************** and Backtransforming it to real Space ********************************
David Rohr's avatar
David Rohr committed
817
	// **************************************************************************************
818

819
820
	cuda_custom_timeslot("Convolution", 1);

821
	mycomplex_t* tmp = param.fft_scratch_complex[omp_get_thread_num()];
822

David Rohr's avatar
David Rohr committed
823
	// **** Multiplying FFTmap with corresponding kernel ****
824
	const mycomplex_t* refCTF = &param.refCTF[iConv * param.FFTMapSize];
825
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
826
	{
827
828
829
		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";
830
831
	}

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

835
	sumsquareC = 0;
836
	if (FFTAlgo)
837
	{
838
839
		int jloopend = param.param_device.NumberFFTPixels1D;
		if ((param.param_device.NumberPixels & 1) == 0) jloopend--;
840
841
		for(int i = 0; i < param.param_device.NumberPixels; i++)
		{
842
			for (int j = 1;j < jloopend;j++)
843
844
845
846
847
848
			{
				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];
849
850
			if ((param.param_device.NumberPixels & 1) == 0)
			{
851
				k += param.param_device.NumberFFTPixels1D - 1;
852
853
				sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1];
			}
854
		}
David Rohr's avatar
David Rohr committed
855

856
857
		myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
		sumsquareC = sumsquareC / norm2;
858
	}
859
860
861
862
863
864
865
	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);
866

867
868
869
870
871
872
873
874
875
		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
876

877
878
	cuda_custom_timeslot_end;

879
	return(0);
880
881
}

882
int bioem::calcross_cor(myfloat_t* localmap, myfloat_t& sum, myfloat_t& sumsquare)
883
{
David Rohr's avatar
David Rohr committed
884
	// *********************** Routine to calculate Cross correlations***********************
885

886
887
	sum = 0.0;
	sumsquare = 0.0;
888
889
890
891
892
	for (int i = 0; i < param.param_device.NumberPixels; i++)
	{
		for (int j = 0; j < param.param_device.NumberPixels; j++)
		{
			// Calculate Sum of pixels
893
			sum += localmap[i * param.param_device.NumberPixels + j];
894
			// Calculate Sum of pixels squared
895
			sumsquare += localmap[i * param.param_device.NumberPixels + j] * localmap[i * param.param_device.NumberPixels + j];
896
897
898
		}
	}
	return(0);
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
}

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

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

int bioem::deviceFinishRun()
{
	return(0);
}
915
916
917
918
919
920
921
922
923
924

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