bioem.cpp 37.3 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
360
361
362
363
364
365
366
367
368
369
370
371
		 if (param.param_device.writeCC)
        	{      int  cc=0;
                	 for (int cent_x = 0; cent_x < param.param_device.NumberPixels -param.param_device.CCdisplace ; cent_x = cent_x + param.param_device.CCdisplace)
                         {
                                 for (int cent_y = 0; cent_y < param.param_device.NumberPixels - param.param_device.CCdisplace ; cent_y = cent_y + param.param_device.CCdisplace)
                                 {
                                 bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
                                 pProbCC.forCC = 0.0;
 
                                 cc++;
                                 }
                         }
		}                 
372
	}
David Rohr's avatar
David Rohr committed
373
	// **************************************************************************************
374
375
	deviceStartRun();

David Rohr's avatar
David Rohr committed
376
	// ******************************** MAIN CYCLE ******************************************
David Rohr's avatar
David Rohr committed
377

378
	mycomplex_t* proj_mapsFFT;
379
	myfloat_t* conv_map = NULL;
380
	mycomplex_t* conv_mapFFT;
381
	myfloat_t sumCONV, sumsquareCONV;
382
383

	//allocating fftw_complex vector
384
385
	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);
386
	conv_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
387
	if (!FFTAlgo) conv_map = (myfloat_t*) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
388

David Rohr's avatar
David Rohr committed
389
	HighResTimer timer, timer2;
390

391
	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
392
393
394
395

	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
396

397
	for (int iOrientAtOnce = iOrientStart; iOrientAtOnce < iOrientEnd; iOrientAtOnce += nProjectionsAtOnce)
398
	{
David Rohr's avatar
David Rohr committed
399
400
		// ***************************************************************************************
		// ***** Creating Projection for given orientation and transforming to Fourier space *****
David Rohr's avatar
David Rohr committed
401
402
		if (DebugOutput >= 1) timer2.ResetStart();
		if (DebugOutput >= 2) timer.ResetStart();
403
404
405
406
407
408
409
		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);
410

411
		for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
412
		{
413
414
415
416
417
418
			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 ***
419

420
421
422
				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);
423

424
425
426
427
				// ***************************************************************************************
				// *** Comparing each calculated convoluted map with all experimental maps ***
				if (DebugOutput >= 2) timer.ResetStart();
				compareRefMaps(iOrient, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
428

429
430
431
432
433
434
435
436
437
438
439
440
441
442
				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)
443
			{
444
445
				printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient, timer2.GetCurrentElapsedTime(), mpi_rank);
				timer2.ResetStart();
446
			}
447
448
449
		}
	}
	//deallocating fftw_complex vector
450
	myfftw_free(proj_mapsFFT);
451
	myfftw_free(conv_mapFFT);
452
	if (!FFTAlgo) myfftw_free(conv_map);
David Rohr's avatar
David Rohr committed
453

454
455
	deviceFinishRun();

David Rohr's avatar
David Rohr committed
456
	// ************* Writing Out Probabilities ***************
457

David Rohr's avatar
David Rohr committed
458
	// *** Angular Probability ***
David Rohr's avatar
David Rohr committed
459

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

514
			if (mpi_rank == 0)
David Rohr's avatar
David Rohr committed
515
			{
516
517
518
519
520
521
				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
522
			}
David Rohr's avatar
David Rohr committed
523

524
525
526
527
			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
528
		}
529
530
531

		//Angle Reduction and Probability summation for individual angles
		if (param.param_device.writeAngles)
David Rohr's avatar
David Rohr committed
532
		{
533
534
535
536
537
			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
538
			{
539
					tmp1[i] = pProb.getProbMap(i).Constoadd;
David Rohr's avatar
David Rohr committed
540
			}
541
			MPI_Allreduce(tmp1, tmp2, count, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
David Rohr's avatar
David Rohr committed
542
543
544
545
546
			for (int i = 0;i < RefMap.ntotRefMap;i++)
			{
				for (int j = 0;j < param.nTotGridAngles;j++)
				{
					bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
547
548
549
550
551
552
553
554
555
556
557
558
559
560
					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
561
562
				}
			}
563
564
565
			delete[] tmp1;
			delete[] tmp2;
			delete[] tmp3;
David Rohr's avatar
David Rohr committed
566
		}
567
	}
David Rohr's avatar
David Rohr committed
568
#endif
569

David Rohr's avatar
David Rohr committed
570
	if (mpi_rank == 0)
571
	{
David Rohr's avatar
David Rohr committed
572
573
574
575
576
		ofstream angProbfile;
		if(param.param_device.writeAngles)
		{
			angProbfile.open ("ANG_PROB");
		}
577
		ofstream ccProbfile;
578
579
580
581
		if(param.param_device.writeCC)
		{
			ccProbfile.open ("CROSS_CORRELATION");
		}
David Rohr's avatar
David Rohr committed
582
583
584
585
586
587
588

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

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

David Rohr's avatar
David Rohr committed
592
			outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
593

David Rohr's avatar
David Rohr committed
594
595
596
597
598
599
600
601
602
			// *** 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
603
604
605
606
			outputProbFile << pProbMap.max.max_prob_cent_y << " " ;
	                outputProbFile << pProbMap.max.max_prob_norm  << " " ;
        	        outputProbFile << pProbMap.max.max_prob_mu ;
                	outputProbFile << "\n";
607

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

David Rohr's avatar
David Rohr committed
610
			if(param.param_device.writeAngles)
611
			{
David Rohr's avatar
David Rohr committed
612
613
614
				for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
				{
					bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
615

David Rohr's avatar
David Rohr committed
616
617
					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";
				}
618
			}
619
620
			if(param.param_device.writeCC)
			{
Pilar Cossio's avatar
Pilar Cossio committed
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
			int  cc=0;
                        int halfPix,rx,ry;
	                myfloat_t localcc[ (param.param_device.NumberPixels+1) * (param.param_device.NumberPixels+1) ];

                        halfPix = param.param_device.NumberPixels / 2 ;
			// Ordering the centers of the Cross Correlation

		for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels -param.param_device.CCdisplace ; rx = rx + param.param_device.CCdisplace)
                         {
                                 for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels - param.param_device.CCdisplace ; ry = ry + param.param_device.CCdisplace)
                                 {
				localcc[ rx * param.param_device.NumberPixels + ry ] = 0.0;
                                }
                        }

636
637
638
639
			for (int cent_x = 0; cent_x < param.param_device.NumberPixels -param.param_device.CCdisplace ; cent_x = cent_x + param.param_device.CCdisplace)
                         {
                                 for (int cent_y = 0; cent_y < param.param_device.NumberPixels - param.param_device.CCdisplace ; cent_y = cent_y + param.param_device.CCdisplace)
                                 {
Pilar Cossio's avatar
Pilar Cossio committed
640
				//localcc[ rx * param.param_device.NumberPixels + ry ] = 0.0;
641
						bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
Pilar Cossio's avatar
Pilar Cossio committed
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
						if(cent_x < halfPix && cent_y < halfPix){
					//	ccProbfile << " " << iRefMap << " " << (myfloat_t) halfPix  - cent_x << " " << halfPix - cent_y << " " << pProbCC.forCC <<"\n";
						rx = halfPix  - cent_x;
                                                ry = halfPix  - cent_y;}
						if(cent_x >= halfPix && cent_y < halfPix){
                                          //      ccProbfile << " " << iRefMap << " " << (myfloat_t) 3 * halfPix  - cent_x << " " << halfPix - cent_y << " " << pProbCC.forCC <<"\n"; 
						rx = 3 * halfPix  - cent_x;
                                                ry = halfPix  - cent_y;}
                                                if(cent_x < halfPix && cent_y >= halfPix){
                                          //      ccProbfile << " " << iRefMap << " " << (myfloat_t) halfPix  - cent_x << " " << 3 * halfPix - cent_y << " " << pProbCC.forCC <<"\n";
						rx = halfPix  - cent_x;
                                                ry = 3 * halfPix  - cent_y;}
                                                 if(cent_x >= halfPix && cent_y >= halfPix){
                                        //        ccProbfile << " " << iRefMap << " " << 3* halfPix  - cent_x << " " << 3 * halfPix - cent_y << " " << pProbCC.forCC <<"\n";
						rx = 3 * halfPix  - cent_x;
                                                ry = 3 * halfPix  - cent_y;}
//						cout << " TT " << cent_x << " " << rx << " " << cent_y << " " << ry << " " <<  pProbCC.forCC << "\n";
						localcc[ rx * param.param_device.NumberPixels + ry ] = pProbCC.forCC;
                                  cc++;
	 			}
                                  //              ccProbfile << "\n";
663
			}
Pilar Cossio's avatar
Pilar Cossio committed
664
665
666
667
668
669
670
671
672
		         for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels -param.param_device.CCdisplace ; rx = rx + param.param_device.CCdisplace)
                         {
                                 for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels - param.param_device.CCdisplace ; ry = ry + param.param_device.CCdisplace)
                                 {
                                    ccProbfile << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ; 
				}
				ccProbfile << "\n";
			}			
  
673
                       }
674
		}
675

David Rohr's avatar
David Rohr committed
676
677
678
679
		if(param.param_device.writeAngles)
		{
			angProbfile.close();
		}
680

681
682
683
684
		if(param.param_device.writeCC)
		{
			ccProbfile.close();
		}
685

David Rohr's avatar
David Rohr committed
686
		outputProbFile.close();
687
	}
688

689
	return(0);
690
691
}

692
int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
693
{
David Rohr's avatar
David Rohr committed
694
695
	//***************************************************************************************
	//***** BioEM routine for comparing reference maps to convoluted maps *****
696
	if (FFTAlgo)
697
	{
David Rohr's avatar
David Rohr committed
698
		//With FFT Algorithm
699
		#pragma omp parallel for schedule(dynamic, 1)
700
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
701
		{
702
			const int num = omp_get_thread_num();
703
			calculateCCFFT(iRefMap, iOrient, iConv, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
704
705
706
		}
	}
	else
707
	{
David Rohr's avatar
David Rohr committed
708
		//Without FFT Algorithm
709
		#pragma omp parallel for schedule(dynamic, 1)
710
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
711
		{
712
			compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, conv_map, pProb, param.param_device, RefMap);
713
714
715
716
717
		}
	}
	return(0);
}

718
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)
719
{
David Rohr's avatar
David Rohr committed
720
721
	//***************************************************************************************
	//***** Calculating cross correlation in FFTALGOrithm *****
Pilar Cossio's avatar
Pilar Cossio committed
722

723
	const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize];
724
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
725
	{
726
727
		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];
728
729
	}

730
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
731

732
	doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
Pilar Cossio's avatar
Pilar Cossio committed
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747

#ifdef PILAR_DEBUG
         if (param.param_device.writeCC)
                {      int  cc=0;
                         for (int cent_x = 0; cent_x < param.param_device.NumberPixels -param.param_device.CCdisplace ; cent_x = cent_x + param.param_device.CCdisplace)
                         {
                                 for (int cent_y = 0; cent_y < param.param_device.NumberPixels - param.param_device.CCdisplace ; cent_y = cent_y + param.param_device.CCdisplace)
                                 {
                                       cout << "CHECKCC " << " " << cent_x << " " << cent_y <<" " << lCC[cent_x * param.param_device.NumberPixels + cent_y] / (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels ) << "\n";
                                 cc++;
                                 }
                         }
                }
#endif

748
}
749

750
int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
751
{
David Rohr's avatar
David Rohr committed
752
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
753
754
	// ****  BioEM Create Projection routine in Euler angle predefined grid******************
	// ********************* and turns projection into Fourier space ************************
David Rohr's avatar
David Rohr committed
755
	// **************************************************************************************
756

757
758
	cuda_custom_timeslot("Projection", 0);

759
760
	myfloat3_t RotatedPointsModel[Model.nPointsModel];
	myfloat_t rotmat[3][3];
761
	myfloat_t alpha, gam, beta;
762
	myfloat_t* localproj;
763

764
	localproj = param.fft_scratch_real[omp_get_thread_num()];
765
	memset(localproj, 0, param.param_device.NumberPixels * param.param_device.NumberPixels * sizeof(*localproj));
766

767
768
769
	alpha = param.angles[iMap].pos[0];
	beta = param.angles[iMap].pos[1];
	gam = param.angles[iMap].pos[2];
770

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

David Rohr's avatar
David Rohr committed
773
	// ********** Creat Rotation with pre-defiend grid of orientations**********
774
775
776
777
778
779
780
781
782
783
784
	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++)
785
	{
786
787
788
		RotatedPointsModel[n].pos[0] = 0.0;
		RotatedPointsModel[n].pos[1] = 0.0;
		RotatedPointsModel[n].pos[2] = 0.0;
789
	}
790
	for(int n = 0; n < Model.nPointsModel; n++)
791
	{
792
		for(int k = 0; k < 3; k++)
793
		{
794
			for(int j = 0; j < 3; j++)
795
			{
796
				RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.points[n].point.pos[j];
797
798
799
800
801
802
			}
		}
	}

	int i, j;

Pilar Cossio's avatar
Pilar Cossio committed
803
	//********** Projection with radius ***************
804
805
806
807
	if(param.doaaradius)
	{
		int irad;
		myfloat_t dist, rad2;
Pilar Cossio's avatar
Pilar Cossio committed
808

809
810
811
812
813
		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
814
815
			//Getting the radius
			irad=int( Model.points[n].radius / param.pixelSize );
816
			rad2= Model.points[n].radius * Model.points[n].radius;
Pilar Cossio's avatar
Pilar Cossio committed
817

818
819
			//Projecting over the radius
			for(int ii= i - irad; ii < i + irad; ii++)
Pilar Cossio's avatar
Pilar Cossio committed
820
821
822
			{	
				for(int jj = j - irad; jj < j + irad; jj++)
				{
823
824
					dist= ( (myfloat_t) (ii-i)*(ii-i)+(jj-j)*(jj-j) ) *  param.pixelSize ;
					if( dist < rad2 )
Pilar Cossio's avatar
Pilar Cossio committed
825
					{
826
827
						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
828
					}
829
830
				}
			}
Pilar Cossio's avatar
Pilar Cossio committed
831
832
833
		}
	}
	else 
834
	{
835
		// ************ Simple Projection over the Z axis********************
Pilar Cossio's avatar
Pilar Cossio committed
836
837
		for(int n = 0; n < Model.nPointsModel; n++)
		{
838
839
840
			//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);
841

842
843
844
845
846
			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;
			}
847

848
849
			localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density / Model.NormDen;
		}
850
851
	}

David Rohr's avatar
David Rohr committed
852
	// **** Output Just to check****
853
#ifdef PILAR_DEBUG
854
	if(iMap == 10)
855
856
857
858
859
860
	{
		ofstream myexamplemap;
		ofstream myexampleRot;
		myexamplemap.open ("MAP_i10");
		myexampleRot.open ("Rot_i10");
		myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
861
		for(int k = 0; k < param.param_device.NumberPixels; k++)
862
		{
863
			for(int j = 0; j < param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j << " " << localproj[k * param.param_device.NumberPixels + j];
864
865
		}
		myexamplemap << " \n";
866
		for(int n = 0; n < Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
867
868
869
		myexamplemap.close();
		myexampleRot.close();
	}
870
#endif
871

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

876
877
	cuda_custom_timeslot_end;

878
879
880
	return(0);
}

881
int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
882
{
David Rohr's avatar
David Rohr committed
883
884
	// **************************************************************************************
	// ****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
David Rohr's avatar
David Rohr committed
885
886
	// **************** calculated Projection with convoluted precalculated Kernel***********
	// *************** and Backtransforming it to real Space ********************************
David Rohr's avatar
David Rohr committed
887
	// **************************************************************************************
888

889
890
	cuda_custom_timeslot("Convolution", 1);

891
	mycomplex_t* tmp = param.fft_scratch_complex[omp_get_thread_num()];
892

David Rohr's avatar
David Rohr committed
893
	// **** Multiplying FFTmap with corresponding kernel ****
894
	const mycomplex_t* refCTF = &param.refCTF[iConv * param.FFTMapSize];
895
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
896
	{
897
898
899
		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";
900
901
	}

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

905
	sumsquareC = 0;
906
	if (FFTAlgo)
907
	{
908
909
		int jloopend = param.param_device.NumberFFTPixels1D;
		if ((param.param_device.NumberPixels & 1) == 0) jloopend--;
910
911
		for(int i = 0; i < param.param_device.NumberPixels; i++)
		{
912
			for (int j = 1;j < jloopend;j++)
913
914
915
916
917
918
			{
				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];
919
920
			if ((param.param_device.NumberPixels & 1) == 0)
			{
921
				k += param.param_device.NumberFFTPixels1D - 1;
922
923
				sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1];
			}
924
		}
David Rohr's avatar
David Rohr committed
925

926
927
		myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
		sumsquareC = sumsquareC / norm2;
928
	}
929
930
931
932
933
934
935
	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);
936

937
938
939
940
941
942
943
944
945
		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
946

947
948
	cuda_custom_timeslot_end;

949
	return(0);
950
951
}

952
int bioem::calcross_cor(myfloat_t* localmap, myfloat_t& sum, myfloat_t& sumsquare)
953
{
David Rohr's avatar
David Rohr committed
954
	// *********************** Routine to calculate Cross correlations***********************
955

956
957
	sum = 0.0;
	sumsquare = 0.0;
958
959
960
961
962
	for (int i = 0; i < param.param_device.NumberPixels; i++)
	{
		for (int j = 0; j < param.param_device.NumberPixels; j++)
		{
			// Calculate Sum of pixels
963
			sum += localmap[i * param.param_device.NumberPixels + j];
964
			// Calculate Sum of pixels squared
965
			sumsquare += localmap[i * param.param_device.NumberPixels + j] * localmap[i * param.param_device.NumberPixels + j];
966
967
968
		}
	}
	return(0);
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
}

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

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

int bioem::deviceFinishRun()
{
	return(0);
}
985
986
987
988
989
990
991
992
993
994

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

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