bioem.cpp 21.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#include <fstream>
#include <boost/program_options.hpp>
#include <iostream>
#include <algorithm>
#include <iterator>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cmath>
#include <omp.h>

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

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

#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)
{
32
33
	copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
	return os;
34
35
36
37
}

bioem::bioem()
{
38
	FFTAlgo = getenv("FFTALGO") == NULL ? 0 : atoi(getenv("FFTALGO"));
39
40
41
42
}

bioem::~bioem()
{
David Rohr's avatar
David Rohr committed
43

44
45
46
47
}

int bioem::configure(int ac, char* av[])
{
David Rohr's avatar
David Rohr committed
48
49
50
51
52
	// **************************************************************************************
	// **** Configuration Routine using boost for extracting parameters, models and maps ****
	// **************************************************************************************
	// ****** And Precalculating necessary grids, map crosscorrelations and kernels  ********
	// *************************************************************************************
53

David Rohr's avatar
David Rohr committed
54
	// *** Inizialzing default variables ***
55
56
57
	std::string infile, modelfile, mapfile;
	Model.readPDB = false;
	param.writeAngles = false;
58
59
	param.dumpMap = false;
	param.loadMap = false;
David Rohr's avatar
David Rohr committed
60
61
	RefMap.readMRC = false;
	RefMap.readMultMRC = false;
62

David Rohr's avatar
David Rohr committed
63
	// *************************************************************************************
64
	cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
David Rohr's avatar
David Rohr committed
65
	// *************************************************************************************
66

David Rohr's avatar
David Rohr committed
67
	// ********************* Command line reading input with BOOST ************************
68

69
70
	try {
		po::options_description desc("Command line inputs");
David Rohr's avatar
David Rohr committed
71
72
73
74
75
76
77
78
79
80
81
		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")
		;
David Rohr's avatar
David Rohr committed
82

83
84
85
86
87
88

		po::positional_options_description p;
		p.add("Inputfile", -1);
		p.add("Modelfile", -1);
		p.add("Particlesfile", -1);
		p.add("ReadPDB", -1);
David Rohr's avatar
David Rohr committed
89
90
		p.add("ReadMRC", -1);
		p.add("ReadMultipleMRC", -1);
91
92
93
		p.add("DumpMaps", -1);
		p.add("LoadMapDump", -1);

94
95
96
97
98
99
100
		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;
101
			return 1;
102
103
104
105
		}
		if (vm.count("help")) {
			cout << "Usage: options_description [options]\n";
			cout << desc;
106
			return 1;
107
108
109
110
111
		}

		if (vm.count("Inputfile"))
		{
			cout << "Input file is: ";
112
113
			cout << vm["Inputfile"].as< std::string >() << "\n";
			infile = vm["Inputfile"].as< std::string >();
114
115
116
117
118
		}
		if (vm.count("Modelfile"))
		{
			cout << "Model file is: "
				 << vm["Modelfile"].as<  std::string  >() << "\n";
119
			modelfile = vm["Modelfile"].as<  std::string  >();
120
121
122
123
124
		}

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

David Rohr's avatar
David Rohr committed
128
129
130
131
132
133
134
135
136
137
138
		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;
		}
139
140
141
142

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

		if (vm.count("LoadMapDump"))
		{
			cout << "Loading Map dump.\n";
149
			param.loadMap = true;
150
151
152
153
154
155
		}

		if (vm.count("Particlesfile"))
		{
			cout << "Paricle file is: "
				 << vm["Particlesfile"].as< std::string >() << "\n";
156
			mapfile = vm["Particlesfile"].as< std::string >();
157
158
159
160
161
162
163
		}
	}
	catch(std::exception& e)
	{
		cout << e.what() << "\n";
		return 1;
	}
Pilar Cossio's avatar
Pilar Cossio committed
164
165
166
167
168
        //check for consitency in multiple MRCs
        if(  RefMap.readMultMRC && not(RefMap.readMRC) ){
         cout << "For Multiple MRCs command --ReadMRC is necesary too";
         exit(1);
        }
David Rohr's avatar
David Rohr committed
169
	// ********************* Reading Parameter Input ***************************
170
171
172
173
	// copying inputfile to param class
	param.fileinput = infile.c_str();
	param.readParameters();

David Rohr's avatar
David Rohr committed
174
	// ********************* Reading Model Input ******************************
175
176
177
178
	// copying modelfile to model class
	Model.filemodel = modelfile.c_str();
	Model.readModel();

David Rohr's avatar
David Rohr committed
179
180
	// ********************* Reading Particle Maps Input **********************
	// ********* HERE: PROBLEM if maps dont fit on the memory!! ***************
181
	// copying mapfile to ref map class
182
	param.filemap = mapfile.c_str();
183
184
	RefMap.readRefMaps(param);

David Rohr's avatar
David Rohr committed
185
	// ****************** Precalculating Necessary Stuff *********************
186
	precalculate();
David Rohr's avatar
David Rohr committed
187

188
189
190
191
192
	if (getenv("BIOEM_DEBUG_BREAK"))
	{
		param.nTotGridAngles = atoi(getenv("BIOEM_DEBUG_BREAK"));
		param.nTotCTFs = atoi(getenv("BIOEM_DEBUG_BREAK"));
	}
David Rohr's avatar
David Rohr committed
193

194
	pProb.init(RefMap.ntotRefMap, param.nTotGridAngles, *this);
195

196
197
	deviceInit();

198
	return(0);
199
200
}

201
202
203
void bioem::cleanup()
{
	//Deleting allocated pointers
204
	free_device_host(pProb.ptr);
205
206
207
	RefMap.freePointers();
}

208
209
int bioem::precalculate()
{
David Rohr's avatar
David Rohr committed
210
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
211
	// **Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels**
David Rohr's avatar
David Rohr committed
212
	// **************************************************************************************
213

214
215
	// Generating Grids of orientations
	param.CalculateGridsParam();
216

217
218
	// Precalculating CTF Kernels stored in class Param
	param.CalculateRefCTF();
219

220
221
	//Precalculate Maps
	RefMap.precalculate(param, *this);
222

223
	return(0);
224
225
226
227
}

int bioem::run()
{
David Rohr's avatar
David Rohr committed
228
229
230
	// **************************************************************************************
	// **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
	// **************************************************************************************
231

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

235
	printf("\tInitializing Probabilities\n");
236
237
238
	// Inizialzing Probabilites to zero and constant to -Infinity
	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
239
240
241
242
243
		bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

		pProbMap.Total = 0.0;
		pProbMap.Constoadd = -9999999;
		pProbMap.max_prob = -9999999;
244
		for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
245
		{
246
247
248
249
			bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);

			pProbAngle.forAngles = 0.0;
			pProbAngle.ConstAngle = -99999999;
250
251
		}
	}
David Rohr's avatar
David Rohr committed
252
	// **************************************************************************************
253
	deviceStartRun();
254
255
256
257
258
259
260
261
262
263
	{
		const int count = omp_get_max_threads();
		localCCT = new mycomplex_t*[count];
		lCC = new myfloat_t*[count];
		for (int i = 0;i < count;i++)
		{
			localCCT[i] = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
			lCC[i] = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
		}
	}
264

David Rohr's avatar
David Rohr committed
265
	// ******************************** MAIN CYCLE ******************************************
David Rohr's avatar
David Rohr committed
266

David Rohr's avatar
David Rohr committed
267
	// *** Declaring Private variables for each thread *****
268
	mycomplex_t* proj_mapFFT;
269
	myfloat_t* conv_map = new myfloat_t[param.param_device.NumberPixels * param.param_device.NumberPixels];
270
	mycomplex_t* conv_mapFFT;
271
	myfloat_t sumCONV, sumsquareCONV;
272
273

	//allocating fftw_complex vector
274
275
	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);
276
277
278
279

	HighResTimer timer;

	printf("\tMain Loop (GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)²), Pixels %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);
Pilar Cossio's avatar
Pilar Cossio committed
280
//	printf("\tInner Loop Count (%d %d %d) %lld\n", param.param_device.maxDisplaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, (long long int) (param.param_device.NumberPixels * param.param_device.NumberPixels * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1) * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1)));
281
	for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
282
	{
David Rohr's avatar
David Rohr committed
283
284
		// ***************************************************************************************
		// ***** Creating Projection for given orientation and transforming to Fourier space *****
285
		timer.ResetStart();
286
287
		createProjection(iOrient, proj_mapFFT);
		printf("Time Projection %d: %f\n", iOrient, timer.GetCurrentElapsedTime());
288

David Rohr's avatar
David Rohr committed
289
290
		// ***************************************************************************************
		// ***** **** Internal Loop over convolutions **** *****
291
292
		for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
		{
293
			printf("\t\tConvolution %d %d\n", iOrient, iConv);
David Rohr's avatar
David Rohr committed
294
			// *** Calculating convolutions of projection map and crosscorrelations ***
295

296
			timer.ResetStart();
297
298
			createConvolutedProjectionMap(iOrient, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
			printf("Time Convolution %d %d: %f\n", iOrient, iConv, timer.GetCurrentElapsedTime());
299

David Rohr's avatar
David Rohr committed
300
301
			// ***************************************************************************************
			// *** Comparing each calculated convoluted map with all experimental maps ***
302
			timer.ResetStart();
303
			compareRefMaps(iOrient, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
304

305
306
307
			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 *
308
								  (((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;
309
			const double nGBs = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
310
								(((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;
311
312
			const double nGBs2 = (double) RefMap.ntotRefMap * ((double) param.param_device.NumberPixels * (double) param.param_device.NumberPixels + 8.) * (double) sizeof(myfloat_t) / compTime;

313
			printf("Time Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s)\n", iOrient, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000.);
314
315
316
317
318
		}
	}
	//deallocating fftw_complex vector
	myfftw_free(proj_mapFFT);
	myfftw_free(conv_mapFFT);
319
	delete[] conv_map;
David Rohr's avatar
David Rohr committed
320

321
	deviceFinishRun();
322
323
324
325
326
327
328
329
330
331
	{
		const int count = omp_get_max_threads();
		for (int i = 0;i < count;i++)
		{
			myfftw_free(localCCT[i]);
			myfftw_free(lCC[i]);
		}
		delete[] localCCT;
		delete[] lCC;
	}
332

David Rohr's avatar
David Rohr committed
333
	// ************* Writing Out Probabilities ***************
334

David Rohr's avatar
David Rohr committed
335
	// *** Angular Probability ***
336

337
338
339
340
	// if(param.writeAngles){
	ofstream angProbfile;
	angProbfile.open ("ANG_PROB_iRefMap");
	// }
341

342
343
344
345
	ofstream outputProbFile;
	outputProbFile.open ("Output_Probabilities");
	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
David Rohr's avatar
David Rohr committed
346
		// **** Total Probability ***
347
348
349
		bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

		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";
350
351
352

		outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";

David Rohr's avatar
David Rohr committed
353
		// *** Param that maximize probability****
354
355
356
357
358
359
360
361
362
		outputProbFile << (pProbMap.max_prob + 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_prob_orient].pos[0] << " ";
		outputProbFile << param.angles[pProbMap.max_prob_orient].pos[1] << " ";
		outputProbFile << param.angles[pProbMap.max_prob_orient].pos[2] << " ";
		outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[0] << " ";
		outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[1] << " ";
		outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[2] << " ";
		outputProbFile << pProbMap.max_prob_cent_x << " ";
		outputProbFile << pProbMap.max_prob_cent_y;
363
		outputProbFile << "\n";
364

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

367
		if(param.writeAngles)
368
		{
369
			for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
370
			{
371
				bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
372

373
				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";
374
375
376
			}
		}
	}
377

378
379
	angProbfile.close();
	outputProbFile.close();
380

381
	return(0);
382
383
}

384
int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
385
{
David Rohr's avatar
David Rohr committed
386
387
	//***************************************************************************************
	//***** BioEM routine for comparing reference maps to convoluted maps *****
388
	if (FFTAlgo)
389
	{
David Rohr's avatar
David Rohr committed
390
		//With FFT Algorithm
391
392
		#pragma omp parallel for
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
393
		{
394
			const int num = omp_get_thread_num();
395
			calculateCCFFT(iRefMap, iOrient, iConv, sumC, sumsquareC, localmultFFT, localCCT[num], lCC[num]);
396
397
398
		}
	}
	else
399
	{
David Rohr's avatar
David Rohr committed
400
		//Without FFT Algorithm
401
		#pragma omp parallel for
402
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
403
		{
404
			compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, conv_map, pProb, param.param_device, RefMap);
405
406
407
408
409
		}
	}
	return(0);
}

410
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)
411
{
David Rohr's avatar
David Rohr committed
412
413
	//***************************************************************************************
	//***** Calculating cross correlation in FFTALGOrithm *****
Pilar Cossio's avatar
Pilar Cossio committed
414

415
	const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize];
416
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
417
	{
418
419
		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];
420
421
	}

422
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
423

424
	doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
425
}
426

427
int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
428
{
David Rohr's avatar
David Rohr committed
429
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
430
431
	// ****  BioEM Create Projection routine in Euler angle predefined grid******************
	// ********************* and turns projection into Fourier space ************************
David Rohr's avatar
David Rohr committed
432
	// **************************************************************************************
433
434
435

	myfloat3_t RotatedPointsModel[Model.nPointsModel];
	myfloat_t rotmat[3][3];
436
	myfloat_t alpha, gam, beta;
437
	myfloat_t* localproj;
438

439
	localproj = lCC[omp_get_thread_num()];
440
	memset(localproj, 0, param.param_device.NumberPixels * param.param_device.NumberPixels * sizeof(*localproj));
441

442
443
444
	alpha = param.angles[iMap].pos[0];
	beta = param.angles[iMap].pos[1];
	gam = param.angles[iMap].pos[2];
445

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

David Rohr's avatar
David Rohr committed
448
	// ********** Creat Rotation with pre-defiend grid of orientations**********
449
450
451
452
453
454
455
456
457
458
459
	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++)
460
	{
461
462
463
		RotatedPointsModel[n].pos[0] = 0.0;
		RotatedPointsModel[n].pos[1] = 0.0;
		RotatedPointsModel[n].pos[2] = 0.0;
464
	}
465
	for(int n = 0; n < Model.nPointsModel; n++)
466
	{
467
		for(int k = 0; k < 3; k++)
468
		{
469
			for(int j = 0; j < 3; j++)
470
			{
471
				RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.PointsModel[n].pos[j];
472
473
474
475
476
477
			}
		}
	}

	int i, j;

David Rohr's avatar
David Rohr committed
478
	// ************ Projection over the Z axis********************
479
	for(int n = 0; n < Model.nPointsModel; n++)
480
481
	{
		//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
482
483
		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);
484

485
		localproj[i * param.param_device.NumberPixels + j] += Model.densityPointsModel[n] / Model.NormDen;
486
487
	}

David Rohr's avatar
David Rohr committed
488
	// **** Output Just to check****
489
	if(iMap == 10)
490
491
492
493
494
495
	{
		ofstream myexamplemap;
		ofstream myexampleRot;
		myexamplemap.open ("MAP_i10");
		myexampleRot.open ("Rot_i10");
		myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
496
		for(int k = 0; k < param.param_device.NumberPixels; k++)
497
		{
498
			for(int j = 0; j < param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j << " " << localproj[k * param.param_device.NumberPixels + j];
499
500
		}
		myexamplemap << " \n";
501
		for(int n = 0; n < Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
502
503
504
505
		myexamplemap.close();
		myexampleRot.close();
	}

David Rohr's avatar
David Rohr committed
506
507
	// ***** Converting projection to Fourier Space for Convolution later with kernel****
	// ********** Omp Critical is necessary with FFTW*******
508
	myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localproj, mapFFT);
509
510
511
512

	return(0);
}

513
int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
514
{
David Rohr's avatar
David Rohr committed
515
516
	// **************************************************************************************
	// ****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
David Rohr's avatar
David Rohr committed
517
518
	// **************** calculated Projection with convoluted precalculated Kernel***********
	// *************** and Backtransforming it to real Space ********************************
David Rohr's avatar
David Rohr committed
519
	// **************************************************************************************
520

521
	mycomplex_t* tmp = localCCT[omp_get_thread_num()];
522

David Rohr's avatar
David Rohr committed
523
	// **** Multiplying FFTmap with corresponding kernel ****
524
	const mycomplex_t* refCTF = &param.refCTF[iConv * param.FFTMapSize];
525
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
526
	{
527
528
529
		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";
530
531
	}

532
533
534
	//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);

David Rohr's avatar
David Rohr committed
535
	// **** Bringing convoluted Map to real Space ****
David Rohr's avatar
David Rohr committed
536
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, Mapconv);
537

David Rohr's avatar
David Rohr committed
538
	// *** Calculating Cross-correlations of cal-convoluted map with its self *****
539
540
	sumC = 0;
	sumsquareC = 0;
541
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
542
	{
David Rohr's avatar
David Rohr committed
543
544
		sumC += Mapconv[i];
		sumsquareC += Mapconv[i] * Mapconv[i];
545
	}
David Rohr's avatar
David Rohr committed
546
	// *** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***
547
	// Normalizing
548
549
550
551
	myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
	myfloat_t norm4 = norm2 * norm2;
	sumC = sumC / norm2;
	sumsquareC = sumsquareC / norm4;
552
553

	return(0);
554
555
}

556
int bioem::calcross_cor(myfloat_t* localmap, myfloat_t& sum, myfloat_t& sumsquare)
557
{
David Rohr's avatar
David Rohr committed
558
	// *********************** Routine to calculate Cross correlations***********************
559

560
561
	sum = 0.0;
	sumsquare = 0.0;
562
563
564
565
566
	for (int i = 0; i < param.param_device.NumberPixels; i++)
	{
		for (int j = 0; j < param.param_device.NumberPixels; j++)
		{
			// Calculate Sum of pixels
567
			sum += localmap[i * param.param_device.NumberPixels + j];
568
			// Calculate Sum of pixels squared
569
			sumsquare += localmap[i * param.param_device.NumberPixels + j] * localmap[i * param.param_device.NumberPixels + j];
570
571
572
		}
	}
	return(0);
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
}

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

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

int bioem::deviceFinishRun()
{
	return(0);
}
589
590
591
592
593
594
595
596
597
598

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

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