bioem.cpp 47.6 KB
Newer Older
Pilar Cossio's avatar
License  
Pilar Cossio committed
1
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2
3
4
   < BioEM software for Bayesian inference of Electron Microscopy images>
   Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer.
   Max Planck Institute of Biophysics, Frankfurt, Germany.
5

6
   See license statement for terms of distribution.
Pilar Cossio's avatar
License  
Pilar Cossio committed
7
8
9

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

10
11
#include <mpi.h>

12
13
14
15
16
#define MPI_CHK(expr)							\
  if (expr != MPI_SUCCESS)						\
    {									\
      fprintf(stderr, "Error in MPI function %s: %d\n", __FILE__, __LINE__); \
    }
17

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

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

#include "param.h"
#include "bioem.h"
#include "model.h"
#include "map.h"
40
#include "MersenneTwister.h"
41

42
43
44
45
46
47
#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]);

48
49
50
51
52
53
54
55
56
57
58
59
#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);			\
  }
60
61
62
63
64
#define cuda_custom_timeslot_end nvtxRangePop();
#else
#define cuda_custom_timeslot(name,cid)
#define cuda_custom_timeslot_end
#endif
65

66
67
68
69
70
71
72
#include "bioem_algorithm.h"

using namespace boost;
namespace po = boost::program_options;

using namespace std;

73
74
75
76
77
/* For dvl nodes in hydra with problem in boost
namespace std {
     typedef decltype(nullptr) nullptr_t;
}*/

78
79
80
81
// A helper function of Boost
template<class T>
ostream& operator<<(ostream& os, const vector<T>& v)
{
82
83
  copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
  return os;
84
85
86
87
}

bioem::bioem()
{
88
89
90
  FFTAlgo = getenv("FFTALGO") == NULL ? 1 : atoi(getenv("FFTALGO"));
  DebugOutput = getenv("BIOEM_DEBUG_OUTPUT") == NULL ? 2 : atoi(getenv("BIOEM_DEBUG_OUTPUT"));
  nProjectionsAtOnce = getenv("BIOEM_PROJECTIONS_AT_ONCE") == NULL ? 1 : atoi(getenv("BIOEM_PROJECTIONS_AT_ONCE"));
91
92
93
94
95
96
97
98
}

bioem::~bioem()
{
}

int bioem::configure(int ac, char* av[])
{
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
  // **************************************************************************************
  // **** Configuration Routine using boost for extracting parameters, models and maps ****
  // **************************************************************************************
  // ****** And Precalculating necessary grids, map crosscorrelations and kernels  ********
  // *************************************************************************************

  HighResTimer timer;

  if (mpi_rank == 0)
    {
      // *** Inizialzing default variables ***
      std::string infile, modelfile, mapfile,Inputanglefile,Inputbestmap;
      Model.readPDB = false;
      param.param_device.writeAngles = false;
      param.dumpMap = false;
      param.loadMap = false;
      RefMap.readMRC = false;
      RefMap.readMultMRC = false;
      param.notuniformangles=false;
Pilar Cossio's avatar
Pilar Cossio committed
118
      yesoutfilename=false;
119
120
121
122
123
124
125
126
127
128
129
130
131
132

      // *************************************************************************************
      cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
      // *************************************************************************************

      // ********************* Command line reading input with BOOST ************************

      try {
	po::options_description desc("Command line inputs");
	desc.add_options()
	  ("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file")
	  ("Particlesfile", po::value< std::string>(), "if BioEM (Mandatory) Name of paricles file")
	  ("Inputfile", po::value<std::string>(), "if BioEM (Mandatory) Name of input parameter file") 
	  ("PrintBestCalMap", po::value< std::string>(), "(Optional) Only print best calculated map (file nec.). NO BioEM (!)")
133
	  ("ReadOrientation", po::value< std::string>(), "(Optional) Read orientation list instead of uniform grid (file nec.)")
134
135
136
137
138
	  ("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")
Pilar Cossio's avatar
Pilar Cossio committed
139
	  ("OutputFile",  po::value< std::string>(), "(Optional) For changing the outputfile name")
140
141
142
143
144
145
146
147
148
149
150
	  ("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);
151
	p.add("ReadOrientation",-1);
152
153
154
	p.add("PrintBestCalMap",-1);
	p.add("DumpMaps", -1);
	p.add("LoadMapDump", -1);
Pilar Cossio's avatar
Pilar Cossio committed
155
        p.add("OutputFile",-1);
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170

	po::variables_map vm;
	po::store(po::command_line_parser(ac, av).
		  options(desc).positional(p).run(), vm);
	po::notify(vm);

	if((ac < 4)) {
	  std::cout << desc << std::endl;
	  return 1;
	}
	if (vm.count("help")) {
	  cout << "Usage: options_description [options]\n";
	  cout << desc;
	  return 1;
	}
171

172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
	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;
	  }
189
	if (vm.count("ReadOrientation"))
190
	  {
191
192
193
194
195
196
197
198
	    cout << "Reading Orientation from file: "
		 << vm["ReadOrientation"].as<  std::string  >() << "\n";
	    cout << "Important! if using Quaternions, include \n";
	    cout << "QUATERNIONS keyword in INPUT PARAMETER FILE\n";
	    cout << "First row in file should be the total number of orientations (int)\n";
	    cout << "Euler angle format should be alpha (12.6f) beta (12.6f) gamma (12.6f)\n";
	    cout << "Quaternion format q1 (12.6f) q2 (12.6f) q3 (12.6f) q4 (12.6f)\n";
	    Inputanglefile = vm["ReadOrientation"].as<  std::string  >();
199
200
	    param.notuniformangles=true;
	  }
Pilar Cossio's avatar
Pilar Cossio committed
201
202
203
204
205
206
	if (vm.count("OutputFile"))
          {
	OutfileName = vm["OutputFile"].as< std::string >();
 	cout << "Writing OUTPUT to: " <<  vm["OutputFile"].as<  std::string  >() << "\n";
	yesoutfilename=true;
	}
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
	if (vm.count("PrintBestCalMap"))
	  {
	    cout << "Reading Euler Angles from file: "
		 << vm["PrintBestCalMap"].as<  std::string  >() << "\n";
	    Inputbestmap = vm["PrintBestCalMap"].as<  std::string  >();
	    param.printModel=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 >();
	  }
      }
      catch(std::exception& e)
	{
	  cout << e.what() << "\n";
	  return 1;
	}
David Rohr's avatar
David Rohr committed
251

252
253
      //check for consitency in multiple MRCs
      if(RefMap.readMultMRC && not(RefMap.readMRC))
David Rohr's avatar
David Rohr committed
254
	{
255
256
257
	  cout << "For Multiple MRCs command --ReadMRC is necesary too";
	  exit(1);
	}
David Rohr's avatar
David Rohr committed
258

259
260
261
262
      if(!Model.readPDB){
	cout << "Note: Reading model in simple text format (not PDB)\n";
	cout << "----  x   y   z  radius  density ------- \n";
      } 
David Rohr's avatar
David Rohr committed
263

264
265
266
267
268
      if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
      // ********************* Reading Parameter Input ***************************
      if(!param.printModel){
	// Standard definition for BioEM
	param.readParameters(infile.c_str(),Inputanglefile.c_str());
David Rohr's avatar
David Rohr committed
269

270
271
	// ********************* Reading Particle Maps Input **********************
	RefMap.readRefMaps(param, mapfile.c_str());
David Rohr's avatar
David Rohr committed
272

273

274
275
      } else{
	// Reading parameters for only writting down Best projection
276

277
278
	param.forprintBest(Inputbestmap.c_str());
      }	
279

280
      // ********************* Reading Model Input ******************************
281
      Model.readModel(param, modelfile.c_str());
David Rohr's avatar
David Rohr committed
282

283
      cout << "Remark: look at file COORDREAD to confirm that the Model coordinates are correct\n";
David Rohr's avatar
David Rohr committed
284

285
      if (DebugOutput >= 2 && mpi_rank == 0) printf("Reading Input Data %f\n", timer.GetCurrentElapsedTime());
286
287
288
289
290
291
    
	if(param.param_device.writeCC && mpi_size>1){
	cout << "Exiting::: WRITE CROSS-CORRELATION ONLY VAILD FOR 1 MPI PROCESS\n";
        exit(1);
	}
}
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
#ifdef WITH_MPI
  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);
      if (DebugOutput >= 2 && mpi_rank == 0) printf("MPI Broadcast of Input Data %f\n", timer.GetCurrentElapsedTime());
    }
#endif
309

310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
  // ****************** Precalculating Necessary Stuff *********************

  if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
  param.PrepareFFTs();
  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Prepare FFTs %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
  precalculate();

  if (getenv("BIOEM_DEBUG_BREAK"))
    {
      const int cut = atoi(getenv("BIOEM_DEBUG_BREAK"));
      if (param.nTotGridAngles > cut) param.nTotGridAngles = cut;
      if (param.nTotCTFs > cut) param.nTotCTFs = cut;
    }

  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Precalculate %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
  if(!param.printModel)pProb.init(RefMap.ntotRefMap, param.nTotGridAngles, param.nTotCC, *this);

  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Init Probabilities %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
  deviceInit();
  if (DebugOutput >= 2 && mpi_rank == 0) printf("Time Device Init %f\n", timer.GetCurrentElapsedTime());

  return(0);
344
345
}

346
347
void bioem::cleanup()
{
348
349
350
  //Deleting allocated pointers
  free_device_host(pProb.ptr);
  RefMap.freePointers();
351
352
}

353
354
int bioem::precalculate()
{
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
  // **************************************************************************************
  // **Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels**
  // **************************************************************************************
  HighResTimer timer;

  // Generating Grids of orientations 
  param.CalculateGridsParam();                

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

  if (DebugOutput >= 3)
    {
      printf("\tTime Precalculate CTFs: %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
  //Precalculate Maps
  if(!param.printModel) RefMap.precalculate(param, *this);
  if (DebugOutput >= 3) printf("\tTime Precalculate Maps: %f\n", timer.GetCurrentElapsedTime());

  return(0);
381
382
383
384
}

int bioem::run()
{
David Rohr's avatar
David Rohr committed
385

386
387
388
  // **************************************************************************************
  // ********** Secondary routine for printing out the only best projection ***************
  // **************************************************************************************
David Rohr's avatar
David Rohr committed
389

390
  if(mpi_rank == 0 && param.printModel){ //Only works for 1 MPI process (not parallelized)
391

392
393
394
395
396
    cout << "\nAnalysis for printing best projection::: \n \n" ; 
    mycomplex_t* proj_mapsFFT;
    myfloat_t* conv_map = NULL;
    mycomplex_t* conv_mapFFT;
    myfloat_t sumCONV, sumsquareCONV;
397

398
399
400
    proj_mapsFFT = (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);
    conv_map = (myfloat_t*) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
401

402
    cout << "...... Calculating Projection .......................\n " ;
403

404
    createProjection(0, proj_mapsFFT);
405

406
    cout << "...... Calculating Convolution .......................\n " ;
David Rohr's avatar
David Rohr committed
407

408
    createConvolutedProjectionMap(0, 0, proj_mapsFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
409

410
  }
411

412
413
414
  // **************************************************************************************
  // **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
  // **************************************************************************************
415

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

419
  if (mpi_rank == 0) printf("\tInitializing Probabilities\n");
420

421
422
423
424
  // Inizialzing Probabilites to zero and constant to -Infinity
  for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
    {
      bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);
David Rohr's avatar
David Rohr committed
425

426
      pProbMap.Total = 0.0;
427
428
      pProbMap.Constoadd = -FLT_MAX; //Problem if using double presicion

429
      if (param.param_device.writeAngles)
430
	{
431
432
433
434
435
	  for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
	    {
	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);

	      pProbAngle.forAngles = 0.0;
436
	      pProbAngle.ConstAngle = -FLT_MAX;
437
438
	    }
	}
439
440

    if (param.param_device.writeCC)
441
	{      int  cc=0;
442
	  for (int cent_x = 0; cent_x < param.param_device.NumberPixels; cent_x = cent_x + param.param_device.CCdisplace)
443
	    {
444
	      for (int cent_y = 0; cent_y < param.param_device.NumberPixels; cent_y = cent_y + param.param_device.CCdisplace)
445
		{
446
		  bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
447
448
449
450
451
452
453
454
455
			//Debuggin:: cout << iRefMap << " " << cc << " " << cent_x << " " << cent_y << "\n";

			if(!param.param_device.CCwithBayes) {
					pProbCC.forCC=-FLT_MAX;
                        	}else {
                                        pProbCC.forCC = 0.0;
					pProbCC.ConstCC=-FLT_MAX;
                                }
		  	cc++;
456
		}
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
	    }
	}                 
    }
  // **************************************************************************************
  deviceStartRun();

  // ******************************** MAIN CYCLE ******************************************

  mycomplex_t* proj_mapsFFT;
  myfloat_t* conv_map = NULL;
  mycomplex_t* conv_mapFFT;
  myfloat_t sumCONV, sumsquareCONV;

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

  HighResTimer timer, timer2;

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

  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;

  for (int iOrientAtOnce = iOrientStart; iOrientAtOnce < iOrientEnd; iOrientAtOnce += nProjectionsAtOnce)
    {
      // ***************************************************************************************
      // ***** Creating Projection for given orientation and transforming to Fourier space *****
      if (DebugOutput >= 1) timer2.ResetStart();
      if (DebugOutput >= 2) timer.ResetStart();
      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);
498

499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
      for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
	{
	  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 ***

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

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

	      if (DebugOutput >= 2)
518
		{
519
520
521
522
523
524
525
526
527
		  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);
528
		}
529
530
531
532
533
534
	    }
	  if (DebugOutput >= 1)
	    {
	      printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient, timer2.GetCurrentElapsedTime(), mpi_rank);
	      timer2.ResetStart();
	    }
535
	}
536
537
538
539
540
    }
  //deallocating fftw_complex vector
  myfftw_free(proj_mapsFFT);
  myfftw_free(conv_mapFFT);
  if (!FFTAlgo) myfftw_free(conv_map);
David Rohr's avatar
David Rohr committed
541

542
  deviceFinishRun();
543

544
  // ************* Writing Out Probabilities ***************
545

546
  // *** Angular Probability ***
David Rohr's avatar
David Rohr committed
547

David Rohr's avatar
David Rohr committed
548
#ifdef WITH_MPI
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
  if (mpi_size > 1)
    {
      if (DebugOutput >= 1 && mpi_rank == 0) timer.ResetStart();
      //Reduce Constant and summarize probabilities
      {
	myfloat_t* tmp1 = new myfloat_t[RefMap.ntotRefMap];
	myfloat_t* tmp2 = new myfloat_t[RefMap.ntotRefMap];
	myfloat_t* tmp3 = new myfloat_t[RefMap.ntotRefMap];
	for (int i = 0;i < RefMap.ntotRefMap;i++)
	  {
	    tmp1[i] = pProb.getProbMap(i).Constoadd;
	  }
	MPI_Allreduce(tmp1, tmp2, RefMap.ntotRefMap, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
	for (int i = 0;i < RefMap.ntotRefMap;i++)
	  {
	    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;
David Rohr's avatar
David Rohr committed
571
	{
572
573
574
575
576
577
578
579
580
581
582
	  int* tmpi1 = new int[RefMap.ntotRefMap];
	  int* tmpi2 = new int[RefMap.ntotRefMap];
	  for (int i = 0;i < RefMap.ntotRefMap;i++)
	    {
	      bioem_Probability_map& pProbMap = pProb.getProbMap(i);
	      tmpi1[i] = tmp2[i] <= pProbMap.Constoadd ? mpi_rank : -1;
	    }
	  MPI_Allreduce(tmpi1, tmpi2, RefMap.ntotRefMap, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
	  for (int i = 0;i < RefMap.ntotRefMap;i++)
	    {
	      if (tmpi2[i] == -1)
David Rohr's avatar
David Rohr committed
583
		{
584
		  if (mpi_rank == 0) printf("Error: Could not find highest probability\n");
David Rohr's avatar
David Rohr committed
585
		}
586
	      else if (tmpi2[i] != 0) //Skip if rank 0 already has highest probability
David Rohr's avatar
David Rohr committed
587
		{
588
589
590
591
592
593
594
595
		  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
596
		}
597
598
599
	    }
	  delete[] tmpi1;
	  delete[] tmpi2;
600
	}
601

David Rohr's avatar
David Rohr committed
602
	if (mpi_rank == 0)
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
	  {
	    for (int i = 0;i < RefMap.ntotRefMap;i++)
	      {
		bioem_Probability_map& pProbMap = pProb.getProbMap(i);
		pProbMap.Total = tmp3[i];
		pProbMap.Constoadd = tmp2[i];
	      }
	  }

	delete[] tmp1;
	delete[] tmp2;
	delete[] tmp3;
	if (DebugOutput >= 1 && mpi_rank == 0 && mpi_size > 1) printf("Time MPI Reduction: %f\n", timer.GetCurrentElapsedTime());
      }

      //Angle Reduction and Probability summation for individual angles
      if (param.param_device.writeAngles)
620
	{
621
622
623
624
625
626
627
628
629
630
631
632
	  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++)
	    {
	      tmp1[i] = pProb.getProbMap(i).Constoadd;
	    }
	  MPI_Allreduce(tmp1, tmp2, count, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
	  for (int i = 0;i < RefMap.ntotRefMap;i++)
	    {
	      for (int j = 0;j < param.nTotGridAngles;j++)
David Rohr's avatar
David Rohr committed
633
		{
634
635
		  bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
		  tmp1[i * param.nTotGridAngles + j] = pProbAngle.forAngles * exp(pProbAngle.ConstAngle - tmp2[i * param.nTotGridAngles + j]);
David Rohr's avatar
David Rohr committed
636
		}
637
638
639
640
641
	    }
	  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++)
642
		{
643
644
645
646
647
648
		  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];
		    }
649
		}
650
651
652
653
654
655
656
	    }
	  delete[] tmp1;
	  delete[] tmp2;
	  delete[] tmp3;
	}
    }
#endif
David Rohr's avatar
David Rohr committed
657

658
659
  if (mpi_rank == 0)
    {
660
661
662
 
// Output for Angle Probability File
	ofstream angProbfile;
663
664
665
      if(param.param_device.writeAngles)
	{
	  angProbfile.open ("ANG_PROB");
666
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
Pilar Cossio's avatar
Pilar Cossio committed
667
668
          if(!param.doquater){ angProbfile <<" RefMap:  MapNumber ; alpha[rad] - beta[rad] - gamma[rad] - log Probability\n" ;}
	  else { angProbfile <<" RefMap:  MapNumber ; q1 - q2 -q3 - log Probability\n" ;};
669
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
670
	}
671
// Output for Cross Correlation File
672
673
674
675
      ofstream ccProbfile;
      if(param.param_device.writeCC)
	{
	  ccProbfile.open ("CROSS_CORRELATION");
676
677
678
	  ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
          ccProbfile <<" RefMap:  MapNumber ; Pixel x - Pixel y - Cross-Correlation \n"; 
          ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
679
680
	}

681
682
// Output for Standard Probability
	ofstream outputProbFile;
Pilar Cossio's avatar
Pilar Cossio committed
683
684
        if(!yesoutfilename)OutfileName="Output_Probabilities";
      	outputProbFile.open (OutfileName.c_str());
685
 	outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n";   
Pilar Cossio's avatar
Pilar Cossio committed
686
687
688
689
690
	outputProbFile << "Notation= RefMap:  MapNumber ; LogProb: natural logarithm of posterior Probability ; Constant: Numerical Const. for adding Probabilities \n";
        if(!param.doquater){
	outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - alpha[rad] - beta[rad] - gamma[rad] - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";}
	else { 
	outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";}
691
692
693
694
	if(param.writeCTF) outputProbFile << " RefMap:  MapNumber ; CTFMaxParm: defocus - b-factor (ref. Penzeck 2010)\n";
	outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";

// Loop over reference maps
695

696
697
698
699
700
      for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
	  // **** Total Probability ***
	  bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

701
	 //Controll for Value of Total Probability
Pilar Cossio's avatar
Pilar Cossio committed
702
          // cout << pProbMap.Total << " " <<  pProbMap.Constoadd << " " << FLT_MAX <<" " << log(FLT_MAX) << "\n";
703
704
          if(pProbMap.Total>1.e-38){

705
	  outputProbFile << "RefMap: " << iRefMap << " LogProb:  "  << 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";
706

707
	  outputProbFile << "RefMap: " << iRefMap << " Maximizing Param: ";
708
709
710

	  // *** 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)) << " ";
Pilar Cossio's avatar
Pilar Cossio committed
711
712
713
	  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] << " [ ] ";
714
715
716
717
718
719
720
	  outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[0] << " [ ] ";
	  outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << " [1./A²] ";
	  outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " [1./A²] ";
	  outputProbFile << pProbMap.max.max_prob_cent_x << " [pix] ";
	  outputProbFile << pProbMap.max.max_prob_cent_y << " [pix] " ;
	  outputProbFile << pProbMap.max.max_prob_norm << " [ ] " ; 
	  outputProbFile << pProbMap.max.max_prob_mu << " [ ] ";
721
	  outputProbFile << "\n";
722
	  }else{ 
Pilar Cossio's avatar
Pilar Cossio committed
723
	 outputProbFile << "PROBLEM!! with Map " << iRefMap << "Numerical Integrated Probability = 0.0; \n - Try re-normalizing experimental map \n-looking at your model density. Maybe RADIUS IS LESS THAN PIXEL SIZE \n If so try KEYWORD NO_PROJECT_RADIUS in parameter inputfile";
724
	 outputProbFile << "RefMap: " << iRefMap << "Most probabily Numerical Constant Out of Float-point range: " << pProbMap.Constoadd  << "\n"; }
725

726
	 // Writing out CTF parameters if requiered
727
	  if(param.writeCTF){
728

729
730
731
732
733
734
735
736
737
	          myfloat_t denomi;
        	  denomi = param.CtfParam[pProbMap.max.max_prob_conv].pos[1] * param.CtfParam[pProbMap.max.max_prob_conv].pos[1] + 
				param.CtfParam[pProbMap.max.max_prob_conv].pos[2] * param.CtfParam[pProbMap.max.max_prob_conv].pos[2];
		  outputProbFile << "RefMap: " << iRefMap << " CTFMaxParam: ";
		  outputProbFile <<  2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel << " [micro-m]; ";
          	outputProbFile << "2*(" << sqrt(4*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi) << ")² [1./A²] \n";
		}

	// Writing Individual Angle probabilities
738
739
740
	  if(param.param_device.writeAngles)
	    {
	      for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
David Rohr's avatar
David Rohr committed
741
		{
742
		  bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
743

744
		  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) << "\n";
745
746
		}
	    }
747
748
749
750
	
	 // Writing Cross-Correlations if requiered
              if(param.param_device.writeCC){

751
752
753
754
755
756
757
758
759
	      int  cc=0;
	      int halfPix;
	      int rx=0;
	      int ry=0;
	      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

760
	      for (int rx = 0; rx < param.param_device.NumberPixels ; rx++)
David Rohr's avatar
David Rohr committed
761
		{
762
		  for (int ry = 0; ry < param.param_device.NumberPixels ; ry++)
763
764
765
		    {
		      localcc[ rx * param.param_device.NumberPixels + ry ] = 0.0;
		    }
David Rohr's avatar
David Rohr committed
766
		}
767

768
	      for (int cent_x = 0; cent_x < param.param_device.NumberPixels ; cent_x = cent_x + param.param_device.CCdisplace)
769
		{
770
		  for (int cent_y = 0; cent_y < param.param_device.NumberPixels ; cent_y = cent_y + param.param_device.CCdisplace)
771
772
773
		    {
		      //localcc[ rx * param.param_device.NumberPixels + ry ] = 0.0;
		      bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
774
775

			// Applying Periodic boundary conditions to the CC
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
		      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";
793
794
795
796
797
		      if(!param.param_device.CCwithBayes){
				localcc[ rx * param.param_device.NumberPixels + ry ] = pProbCC.forCC;
	                      }else{ 
				localcc[ rx * param.param_device.NumberPixels + ry ] = log(pProbCC.forCC)+pProbCC.ConstCC;
 			}
798
799
800
		      cc++;
		    }
		  //              ccProbfile << "\n";
801
		}
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
	  	if(!param.ignoreCCoff){
		      for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels ; rx = rx + param.param_device.CCdisplace)
			{
			  for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry = ry + param.param_device.CCdisplace)
			    {
					ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
//				 cout << " cc " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] <<"\n" ;
		   	 }
				  ccProbfile << "\n";
			}			
  		}else{
                      for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels ; rx++)
                        {
                          for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry++)
                            {
		                 if(localcc[ rx * param.param_device.NumberPixels + ry ]!=0.0) ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
                         }
                                  ccProbfile << "\n";
                        }

		}
823
824
	    }
	}
825

826
827
828
829
830
831
832
833
      if(param.param_device.writeAngles)
	{
	  angProbfile.close();
	}

      if(param.param_device.writeCC)
	{
	  ccProbfile.close();
834
	}
835

836
837
838
839
      outputProbFile.close();
    }

  return(0);
840
841
}

842
int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
843
{
844
845
846
847
848
849
850
  //***************************************************************************************
  //***** BioEM routine for comparing reference maps to convoluted maps *****
  if (FFTAlgo)
    {
      //With FFT Algorithm
#pragma omp parallel for schedule(dynamic, 1)
      for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
851
	{
852
853
	  const int num = omp_get_thread_num();
	  calculateCCFFT(iRefMap, iOrient, iConv, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
854
	}
855
856
857
858
859
860
    }
  else
    {
      //Without FFT Algorithm
#pragma omp parallel for schedule(dynamic, 1)
      for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
861
	{
862
	  compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, conv_map, pProb, param.param_device, RefMap);
863
	}
864
865
    }
  return(0);
866
867
}

868
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)
869
{
870
871
  //***************************************************************************************
  //***** Calculating cross correlation in FFTALGOrithm *****
Pilar Cossio's avatar
Pilar Cossio committed
872

873
874
875
876
877
878
  const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize];
  for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
    {
      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];
    }
879

880
  myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
881

882
  doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
Pilar Cossio's avatar
Pilar Cossio committed
883
884

#ifdef PILAR_DEBUG
885
886
  if (param.param_device.writeCC)
    {      int  cc=0;
887
      for (int cent_x = 0; cent_x < param.param_device.NumberPixels ; cent_x = cent_x + param.param_device.CCdisplace)
888
	{
889
	  for (int cent_y = 0; cent_y < param.param_device.NumberPixels ; cent_y = cent_y + param.param_device.CCdisplace)
890
891
892
893
894
895
	    {
	      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++;
	    }
	}
    }
Pilar Cossio's avatar
Pilar Cossio committed
896
897
#endif

898
}
899

900
int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
901
{
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
  // **************************************************************************************
  // ****  BioEM Create Projection routine in Euler angle predefined grid******************
  // ********************* and turns projection into Fourier space ************************
  // **************************************************************************************

  cuda_custom_timeslot("Projection", 0);

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

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

917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
  if(param.doquater){

  myfloat_t quater[4],cq[4],temp[4];
	myfloat_t a[4];
    //quaternion
	    quater[1]=param.angles[iMap].pos[1];
	    quater[2]=param.angles[iMap].pos[2];
	    quater[3]=param.angles[iMap].pos[0];  
	    quater[0]=param.angles[iMap].quat4;

	// its conjugate
	    cq[0]=quater[0];
	    cq[1]=-quater[1];
	    cq[2]=-quater[2];
	    cq[3]=-quater[3];

for(int n = 0; n < Model.nPointsModel; n++)
    {
	//Rotation with quaterion

	// Initial coords	
	a[0]=0.;
	a[1]=Model.points[n].point.pos[0];
	a[2]=Model.points[n].point.pos[1];
	a[3]=Model.points[n].point.pos[2];

	temp[0]=a[0]*cq[0]-a[1]*cq[1]-a[2]*cq[2]-a[3]*cq[3];
	temp[1]=a[0]*cq[1]+a[1]*cq[0]+a[2]*cq[3]-a[3]*cq[2];
	temp[2]=a[0]*cq[2]-a[1]*cq[3]+a[2]*cq[0]+a[3]*cq[1];
	temp[3]=a[0]*cq[3]+a[1]*cq[2]-a[2]*cq[1]+a[3]*cq[0];

	// First coordinate is zero Going back to real space
        RotatedPointsModel[n].pos[0] = quater[0]*temp[1]+quater[1]*temp[0]+quater[2]*temp[3]-quater[3]*temp[2];
	RotatedPointsModel[n].pos[1] = quater[0]*temp[2]-quater[1]*temp[3]+quater[2]*temp[0]+quater[3]*temp[1];
        RotatedPointsModel[n].pos[2] = quater[0]*temp[3]+quater[1]*temp[2]-quater[2]*temp[1]+quater[3]*temp[0]; 
    
    
	}


	} else{
 // Doing Euler angles instead of Quaternions
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
  alpha = param.angles[iMap].pos[0];
  beta = param.angles[iMap].pos[1];
  gam = param.angles[iMap].pos[2];

  // **** To see how things are going: cout << "Id " << omp_get_thread_num() <<  " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***

  // ********** Creat Rotation with pre-defiend grid of orientations**********
  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++)
    {
      RotatedPointsModel[n].pos[0] = 0.0;
      RotatedPointsModel[n].pos[1] = 0.0;
      RotatedPointsModel[n].pos[2] = 0.0;
    }
  for(int n = 0; n < Model.nPointsModel; n++)
    {
      for(int k = 0; k < 3; k++)
985
	{
986
987
988
989
	  for(int j = 0; j < 3; j++)
	    {
	      RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.points[n].point.pos[j];
	    }
990
	}
991
    }
992

993
994
995
996
997
998
}

if(param.printrotmod) {
 for(int n = 0; n < Model.nPointsModel; n++) cout << "ROTATED " << iMap << " " << n <<" "<< RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " <<  RotatedPointsModel[n].pos[2] << "\n";

}
999
  int i, j;
1000

1001
1002
1003
1004
1005
  //********** Projection with radius ***************
  if(param.doaaradius)
    {
      int irad;
      myfloat_t dist, rad2;
Pilar Cossio's avatar
Pilar Cossio committed
1006

1007
1008
1009
1010
1011
1012
1013
1014
1015
      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);
	  //Getting the radius
	  irad=int( Model.points[n].radius / param.pixelSize );
	  rad2= Model.points[n].radius * Model.points[n].radius;

Pilar Cossio's avatar
Pilar Cossio committed
1016
1017
1018
1019
1020
1021
1022
	 if (i < 0 || j < 0 || i >= param.param_device.NumberPixels || j >= param.param_device.NumberPixels)
            {
              if (DebugOutput >= 0) cout << "WARNING::: Model Point out of Projection map: " << i << ", " << j << "\n";
              continue;
            }


1023
1024
1025
1026
	  //Projecting over the radius
	  for(int ii= i - irad; ii < i + irad; ii++)
	    {	
	      for(int jj = j - irad; jj < j + irad; jj++)
1027
		{
1028
1029
1030
1031
1032
1033
1034
		  dist= ( (myfloat_t) (ii-i)*(ii-i)+(jj-j)*(jj-j) ) *  param.pixelSize ;
		  if( dist < rad2 )
		    {
		      localproj[ii * param.param_device.NumberPixels + jj] += param.pixelSize * 2 * sqrt( rad2 - dist ) * Model.points[n].density
			//							/ Model.NormDen * 3 / (4 * M_PI * rad2 *  Model.points[n].radius); 	
			* 3 / (4 * M_PI * rad2 *  Model.points[n].radius); 
		    }
Pilar Cossio's avatar
Pilar Cossio committed
1035
		}
1036
	    }
Pilar Cossio's avatar
Pilar Cossio committed
1037
	}
1038
1039
1040
1041
1042
    }
  else 
    {
      // ************ Simple Projection over the Z axis********************
      for(int n = 0; n < Model.nPointsModel; n++)
1043
	{
1044
1045
1046
	  //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);
1047

1048
1049
	  if (i < 0 || j < 0 || i >= param.param_device.NumberPixels || j >= param.param_device.NumberPixels)
	    {
Pilar Cossio's avatar
Pilar Cossio committed
1050
	      if (DebugOutput >= 0) cout << "WARNING:::: Model Point out of Projection map: " << i << ", " << j << "\n";
1051
1052
	      continue;
	    }
1053

1054
	  localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density;//Model.NormDen;
1055
	}
1056
    }
1057

1058
  // **** Output Just to check****
1059
#ifdef PILAR_DEBUG
1060
1061
1062
1063
1064
1065
1066
1067
  if(iMap == 10)
    {
      ofstream myexamplemap;
      ofstream myexampleRot;
      myexamplemap.open ("MAP_i10");
      myexampleRot.open ("Rot_i10");
      myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
      for(int k = 0; k < param.param_device.NumberPixels; k++)
1068
	{
1069
	  for(int j = 0; j < param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j << " " << localproj[k * param.param_device.NumberPixels + j];
1070
	}
1071
1072
1073
1074
1075
      myexamplemap << " \n";
      for(int n = 0; n < Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
      myexamplemap.close();
      myexampleRot.close();
    }
1076
#endif
1077

1078
1079
1080
  // ***** Converting projection to Fourier Space for Convolution later with kernel****
  // ********** Omp Critical is necessary with FFTW*******
  myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localproj, mapFFT);
1081

1082
  cuda_custom_timeslot_end;
1083

1084
  return(0);
1085
1086
}

1087
int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
1088
{
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
  // **************************************************************************************
  // ****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
  // **************** calculated Projection with convoluted precalculated Kernel***********
  // *************** and Backtransforming it to real Space ********************************
  // **************************************************************************************


  cuda_custom_timeslot("Convolution", 1);

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

  // **** Multiplying FFTmap with corresponding kernel ****
  const mycomplex_t* refCTF = &param.refCTF[iConv * param.FFTMapSize];
  for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
    {
      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";
    }

  // *** Calculating Cross-correlations of cal-convoluted map with its self *****
  sumC = localmultFFT[0][0];

  sumsquareC = 0;
  if (FFTAlgo)
    {
      int jloopend = param.param_device.NumberFFTPixels1D;
      if ((param.param_device.NumberPixels & 1) == 0) jloopend--;
      for(int i = 0; i < param.param_device.NumberPixels; i++)
	{
	  for (int j = 1;j < jloopend;j++)
	    {
	      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];
	  if ((param.param_device.NumberPixels & 1) == 0)
	    {
	      k += param.param_device.NumberFFTPixels1D - 1;
	      sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1];
	    }
	}
1132

1133
1134
1135
1136
1137
1138
1139
      myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
      sumsquareC = sumsquareC / norm2;
    }
  else
    {
      //FFTW_C2R will destroy the input array, so we have to work on a copy here
      memcpy(tmp, localmultFFT, sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
1140

1141
1142
      // **** Bringing convoluted Map to real Space ****
      myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, Mapconv);
Pilar Cossio's avatar
Pilar Cossio committed
1143

1144
      for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
1145
	{
1146
	  sumsquareC += Mapconv[i] * Mapconv[i];
1147
1148
	}

1149
1150
1151
1152
      myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
      myfloat_t norm4 = norm2 * norm2;
      sumsquareC = sumsquareC / norm4;
    }