bioem.cpp 37.5 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
			outputProbFile << pProbMap.max.max_prob_cent_y << " " ;
Pilar Cossio's avatar
Pilar Cossio committed
604
	                outputProbFile << pProbMap.max.max_prob_norm << " " ; // param.param_device.Ntotpi / param.param_device.Ntotpi << " " ;
Pilar Cossio's avatar
Pilar Cossio committed
605 606
        	        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
						localproj[ii * param.param_device.NumberPixels + jj] += param.pixelSize * 2 * sqrt( rad2 - dist ) * Model.points[n].density
Pilar Cossio's avatar
Pilar Cossio committed
827 828
//							/ Model.NormDen * 3 / (4 * M_PI * rad2 *  Model.points[n].radius); 	
							* 3 / (4 * M_PI * rad2 *  Model.points[n].radius); 
Pilar Cossio's avatar
Pilar Cossio committed
829
					}
830 831
				}
			}
Pilar Cossio's avatar
Pilar Cossio committed
832 833 834
		}
	}
	else 
835
	{
836
		// ************ Simple Projection over the Z axis********************
Pilar Cossio's avatar
Pilar Cossio committed
837 838
		for(int n = 0; n < Model.nPointsModel; n++)
		{
839 840 841
			//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);
842

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

Pilar Cossio's avatar
Pilar Cossio committed
849
			localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density;//Model.NormDen;
850
		}
851 852
	}

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

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

877 878
	cuda_custom_timeslot_end;

879 880 881
	return(0);
}

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

890 891
	cuda_custom_timeslot("Convolution", 1);

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

Pilar Cossio's avatar
Pilar Cossio committed
894 895
 	myfloat_t norm = (myfloat_t) (param.param_device.NumberPixels *  param.param_device.NumberPixels);

David Rohr's avatar
David Rohr committed
896
	// **** Multiplying FFTmap with corresponding kernel ****
897
	const mycomplex_t* refCTF = &param.refCTF[iConv * param.FFTMapSize];
898
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
899
	{
Pilar Cossio's avatar
Pilar Cossio committed
900 901
		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] ) ;
902
		// cout << "GG " << i << " " << j << " " << refCTF[i][0] << " " << refCTF[i][1] <<" " <<lproj[i][0] <<" " <<lproj[i][1] << "\n";
903 904
	}

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

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