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);
}