bioem.cpp 59.1 KB
Newer Older
Pilar Cossio's avatar
License    
Pilar Cossio committed
1
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2
   < BioEM software for Bayesian inference of Electron Microscopy images>
3
4
   Copyright (C) 2016 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp, 
	Volker Lindenstruth and Gerhard Hummer.
Pilar Cossio's avatar
Pilar Cossio committed
5

6
   Max Planck Institute of Biophysics, Frankfurt, Germany.
Pilar Cossio's avatar
Pilar Cossio committed
7
8
   Frankfurt Institute for Advanced Studies, Goethe University Frankfurt, Germany.
   Max Planck Computing and Data Facility, Garching, Germany. 
9

Pilar Cossio's avatar
Pilar Cossio committed
10
   Released under the GNU Public License, v3. 
11
   See license statement for terms of distribution.
Pilar Cossio's avatar
License    
Pilar Cossio committed
12
13

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
Pilar Cossio's avatar
Pilar Cossio committed
14
#ifdef WITH_MPI
15
16
#include <mpi.h>

17
18
19
20
21
#define MPI_CHK(expr)							\
  if (expr != MPI_SUCCESS)						\
    {									\
      fprintf(stderr, "Error in MPI function %s: %d\n", __FILE__, __LINE__); \
    }
Pilar Cossio's avatar
Pilar Cossio committed
22
#endif
23

24
25
#include <fstream>
#include <boost/program_options.hpp>
Pilar Cossio's avatar
Pilar Cossio committed
26
27
28
#include <boost/random/normal_distribution.hpp>
#include <boost/random/uniform_int_distribution.hpp>
#include <boost/random/mersenne_twister.hpp>
29
30
31
32
33
34
35
#include <iostream>
#include <algorithm>
#include <iterator>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cmath>
36

37
#ifdef WITH_OPENMP
38
#include <omp.h>
39
#endif
40
41
42

#include <fftw3.h>
#include <math.h>
43
#include "timer.h"
44
#include "autotuner.h"
45
46
47
48
49
50

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

51
52
53
54
55
56
#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]);

57
58
59
60
61
62
63
64
65
66
67
68
#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);			\
  }
69
70
71
72
73
#define cuda_custom_timeslot_end nvtxRangePop();
#else
#define cuda_custom_timeslot(name,cid)
#define cuda_custom_timeslot_end
#endif
74

75
76
77
78
#include "bioem_algorithm.h"

using namespace boost;
namespace po = boost::program_options;
Pilar Cossio's avatar
Pilar Cossio committed
79
namespace bran= boost::random;
80
81
82

using namespace std;

83
/* For dvl nodes in hydra with problem in boost
84
85
86
   namespace std {
   typedef decltype(nullptr) nullptr_t;
   }*/
87

88
89
90
91
// A helper function of Boost
template<class T>
ostream& operator<<(ostream& os, const vector<T>& v)
{
92
93
  copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
  return os;
94
95
96
97
}

bioem::bioem()
{
98
99
100
  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"));
101
  Autotuning = false;
102
103
104
105
106
107
108
109
}

bioem::~bioem()
{
}

int bioem::configure(int ac, char* av[])
{
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;

Pilar Cossio's avatar
Pilar Cossio committed
118
119


120
  std::string infile, modelfile, mapfile,Inputanglefile,Inputbestmap; 
121
122
123
124
125
126
127
128
129
130
131
  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
132
      yesoutfilename=false;
133
134
135
136
137
138
139
140
141
142
143

      // *************************************************************************************
      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")
Pilar Cossio's avatar
Pilar Cossio committed
144
	  ("Particlesfile", po::value< std::string>(), "if BioEM (Mandatory) Name of particle-image file")
145
	  ("Inputfile", po::value<std::string>(), "if BioEM (Mandatory) Name of input parameter file") 
Pilar Cossio's avatar
Pilar Cossio committed
146
147
	  ("PrintBestCalMap", po::value< std::string>(), "(Optional) Only print best calculated map. NO BioEM (!)")
	  ("ReadOrientation", po::value< std::string>(), "(Optional) Read file name containing orientations")
148
149
150
	  ("ReadPDB", "(Optional) If reading model file in PDB format")
	  ("ReadMRC", "(Optional) If reading particle file in MRC format")
	  ("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
Pilar Cossio's avatar
Pilar Cossio committed
151
152
	  ("DumpMaps", "(Optional) Dump maps after they were read from particle-image file")
	  ("LoadMapDump", "(Optional) Read Maps from dump option")
Pilar Cossio's avatar
Pilar Cossio committed
153
	  ("OutputFile",  po::value< std::string>(), "(Optional) For changing the outputfile name")
154
155
156
157
158
159
160
161
162
163
164
	  ("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);
165
	p.add("ReadOrientation",-1);
166
167
168
	p.add("PrintBestCalMap",-1);
	p.add("DumpMaps", -1);
	p.add("LoadMapDump", -1);
Pilar Cossio's avatar
Pilar Cossio committed
169
        p.add("OutputFile",-1);
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184

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

186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
	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;
	  }
203
	if (vm.count("ReadOrientation"))
204
	  {
205
206
207
208
209
210
211
212
	    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  >();
213
214
	    param.notuniformangles=true;
	  }
Pilar Cossio's avatar
Pilar Cossio committed
215
216
	if (vm.count("OutputFile"))
          {
217
218
219
220
	    OutfileName = vm["OutputFile"].as< std::string >();
	    cout << "Writing OUTPUT to: " <<  vm["OutputFile"].as<  std::string  >() << "\n";
	    yesoutfilename=true;
	  }
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
251
252
253
254
255
256
257
258
259
260
261
262
263
264
	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
265

266
267
      //check for consitency in multiple MRCs
      if(RefMap.readMultMRC && not(RefMap.readMRC))
David Rohr's avatar
David Rohr committed
268
	{
269
270
271
	  cout << "For Multiple MRCs command --ReadMRC is necesary too";
	  exit(1);
	}
David Rohr's avatar
David Rohr committed
272

273
274
275
276
      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
277

278
279
280
281
      if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
      // ********************* Reading Parameter Input ***************************
      if(!param.printModel){
	// Standard definition for BioEM
Pilar Cossio's avatar
Pilar Cossio committed
282
	param.readParameters(infile.c_str());
David Rohr's avatar
David Rohr committed
283

284
285
	// ********************* Reading Particle Maps Input **********************
	RefMap.readRefMaps(param, mapfile.c_str());
David Rohr's avatar
David Rohr committed
286

287

288
289
      } else{
	// Reading parameters for only writting down Best projection
290

291
292
	param.forprintBest(Inputbestmap.c_str());
      }	
293

294
      // ********************* Reading Model Input ******************************
295
      Model.readModel(param, modelfile.c_str());
David Rohr's avatar
David Rohr committed
296

297
      cout << "**NOTE:: look at file COORDREAD to confirm that the Model coordinates are correct\n";
David Rohr's avatar
David Rohr committed
298

299
      if (DebugOutput >= 2 && mpi_rank == 0) printf("Reading Input Data Time: %f\n", timer.GetCurrentElapsedTime());
300
    
301
      if(param.param_device.writeCC && mpi_size>1){
302
303
	cout << "Exiting::: WRITE CROSS-CORRELATION ONLY VAILD FOR 1 MPI PROCESS\n";
        exit(1);
304
      }
Pilar Cossio's avatar
Pilar Cossio committed
305

306
307
      // Generating Grids of orientations 
      if(!param.printModel)param.CalculateGridsParam(Inputanglefile.c_str());
308
    }
309

310
#ifdef WITH_MPI
311

Pilar Cossio's avatar
Pilar Cossio committed
312

313
314

  // ********************* MPI inizialization/ Transfer of parameters******************
315
316
317
318
  if (mpi_size > 1)
    {
      if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
      MPI_Bcast(&param, sizeof(param), MPI_BYTE, 0, MPI_COMM_WORLD);
319
      //We have to reinitialize all pointers !!!!!!!!!!!!
Pilar Cossio's avatar
Pilar Cossio committed
320
      if (mpi_rank != 0) param.angprior = NULL;
Pilar Cossio's avatar
Pilar Cossio committed
321

322
323
      if (mpi_rank != 0)param.angles =  (myfloat3_t*) mallocchk(param.nTotGridAngles  * sizeof (myfloat3_t));
      MPI_Bcast(param.angles, param.nTotGridAngles  * sizeof (myfloat3_t),MPI_BYTE, 0, MPI_COMM_WORLD);
324

Pilar Cossio's avatar
Pilar Cossio committed
325
#ifdef DEBUG
326
327
      for(int n=0;n<param.nTotGridAngles;n++){
	cout << "CHECK: Angle orient " << mpi_rank << " "<< n << " " <<  param.angles[n].pos[0] << " " <<  param.angles[n].pos[1] << " " << param.angles[n].pos[2] << " " << param.angles[n].quat4  << " " << "\n";} 
Pilar Cossio's avatar
Pilar Cossio committed
328
329

#endif
330
      //****refCtf, CtfParam, angles automatically filled by precalculate function bellow
331
332
333
334
335
336
337
338
339

      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());
Pilar Cossio's avatar
Pilar Cossio committed
340

341
342
    }
#endif
343

344
345
346
  // ****************** Precalculating Necessary Stuff *********************
  if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
  param.PrepareFFTs();
347

348
349
350
351
352
353
354
  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Prepare FFTs %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
  precalculate();

355
  // ****************** For debugging *********************
356
357
358
359
360
361
362
363
364
365
366
367
  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();
    }
Luka Stanisic's avatar
Luka Stanisic committed
368

369
370
371
372
373
374
375
376
377
378
379
380
  // ****************** For autotuning **********************
  if ((getenv("GPU") && atoi(getenv("GPU"))) && ((!getenv("GPUWORKLOAD") || (atoi(getenv("GPUWORKLOAD")) == -1))) && (!getenv("BIOEM_DEBUG_BREAK") || (atoi(getenv("BIOEM_DEBUG_BREAK")) > FIRST_STABLE)))
    {
      Autotuning = true;
      if (mpi_rank == 0) printf("Autotuning of GPUWorkload enabled:\n\tAlgorithm %d\n\tRecalibration at every %d projections\n\tComparisons are considered stable after first %d comparisons\n", AUTOTUNING_ALGORITHM, RECALIB_FACTOR, FIRST_STABLE);
    }
  else
    {
      Autotuning = false;
      if (mpi_rank == 0) printf("Autotuning of GPUWorkload disabled\n");
    }

Luka Stanisic's avatar
Luka Stanisic committed
381
382
383
  // ****************** Initializing pointers *********************

  deviceInit();
384
385
386

  if (DebugOutput >= 2 && mpi_rank == 0)
    {
Luka Stanisic's avatar
Luka Stanisic committed
387
      printf("Time Device Init %f\n", timer.GetCurrentElapsedTime());
388
389
      timer.ResetStart();
    }
390

Luka Stanisic's avatar
Luka Stanisic committed
391
  if(!param.printModel)pProb.init(RefMap.ntotRefMap, param.nTotGridAngles, param.nTotCC, *this);
392

Luka Stanisic's avatar
Luka Stanisic committed
393
394
395
396
397
  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Init Probabilities %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
398
399

  return(0);
400
401
}

402
403
void bioem::cleanup()
{
404
  //Deleting allocated pointers
405
406
  free_device_host(pProb.ptr);
  RefMap.freePointers();
407
408
}

409
410
int bioem::precalculate()
{
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
  // **************************************************************************************
  // **Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels**
  // **************************************************************************************
  HighResTimer timer;
  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);
433
434
435
436
}

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

438
439
440
  // **************************************************************************************
  // ********** Secondary routine for printing out the only best projection ***************
  // **************************************************************************************
David Rohr's avatar
David Rohr committed
441

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

444
445
446
447
448
    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;
449

450
451
452
    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);
453

454
    cout << "...... Calculating Projection .......................\n " ;
455

456
    createProjection(0, proj_mapsFFT);
457

458
    cout << "...... Calculating Convolution .......................\n " ;
David Rohr's avatar
David Rohr committed
459

460
    createConvolutedProjectionMap(0, 0, proj_mapsFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
461

462
  }
463

464
465
466
  // **************************************************************************************
  // **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
  // **************************************************************************************
467

468
469
  // **** 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
470

471
  if (mpi_rank == 0) printf("\tInitializing Probabilities\n");
472

Pilar Cossio's avatar
Pilar Cossio committed
473
474
475
476
  // Contros for MPI
  if(mpi_size > param.nTotGridAngles){
    cout << "EXIT: Wrong MPI setup More MPI processes than orientations\n"; exit(1);
  }
477

478
479
480
481
  // 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
482

483
      pProbMap.Total = 0.0;
484
485
      pProbMap.Constoadd = -FLT_MAX; //Problem if using double presicion

486
      if (param.param_device.writeAngles)
487
	{
488
489
490
491
492
	  for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
	    {
	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);

	      pProbAngle.forAngles = 0.0;
493
	      pProbAngle.ConstAngle = -FLT_MAX;
494
495
	    }
	}
496

497
      if (param.param_device.writeCC)
498
	{      int  cc=0;
499
	  for (int cent_x = 0; cent_x < param.param_device.NumberPixels; cent_x = cent_x + param.param_device.CCdisplace)
500
	    {
501
	      for (int cent_y = 0; cent_y < param.param_device.NumberPixels; cent_y = cent_y + param.param_device.CCdisplace)
502
		{
503
		  bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
504
505
506
507
508
509
510
511
512
		  //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++;
513
		}
514
	    }
515
	  if(!FFTAlgo){cout << "Cross correlation calculation must be with enviormental variable FFTALGO=1\n"; exit(1);}
516
517
	}                 
    }
518
519
520

  if(!FFTAlgo){cout << "Remark: Not using FFT algorithm. Not using Prior in B-Env.";}

521
  // **************************************************************************************
522

523
524
525
526
527
528
529
530
531
532
533
  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
534
  //******** Alocating Vectors *************
535
536
537
538
539
540
541
  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;

542
543
544
  /* Autotuning */
  Autotuner aut;
  if (Autotuning)
Luka Stanisic's avatar
Luka Stanisic committed
545
    {
546
547
      aut.Initialize(AUTOTUNING_ALGORITHM, FIRST_STABLE);
      rebalance(aut.Workload());      
Luka Stanisic's avatar
Luka Stanisic committed
548
    }
549

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

552
553


554
555
556
557
  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;

558
559
560

  // **************************Loop Over orientations***************************************

561
562
563
564
565
566
567
  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);
568
569
570

      // **************************Parallel orientations for projections at once***************

571
572
573
574
575
576
#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);
577

Luka Stanisic's avatar
Luka Stanisic committed
578
579
      /* Recalibrate if needed */
      if (Autotuning && ((iOrientAtOnce - iOrientStart) % RECALIB_FACTOR == 0) && ((iOrientEnd - iOrientAtOnce) > RECALIB_FACTOR) && (iOrientAtOnce != iOrientStart))
580
	{
Luka Stanisic's avatar
Luka Stanisic committed
581
582
583
	  aut.Reset();
	  rebalance(aut.Workload());
	}
584

Luka Stanisic's avatar
Luka Stanisic committed
585
586
      for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
	{
587
	  mycomplex_t* proj_mapFFT = &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize];
588

589
	  // ***************************************************************************************
590
591
	  // ***** **** Internal Loop over PSF/CTF convolutions **** *****

592
593
594
595
596
597
598
599
	  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);

600
      	      if ((DebugOutput >= 2) || (Autotuning && aut.Needed(iConv))) timer.ResetStart();
601
602
603
604
605
606
     	      myfloat_t amp,pha,env;

              amp=param.CtfParam[iConv].pos[0];
              pha=param.CtfParam[iConv].pos[1];
              env=param.CtfParam[iConv].pos[2];

607
608
609
	      // ******************Internal loop over Reference images CUDA or OpenMP******************
	      // *** Comparing each calculated convoluted map with all experimental maps ***

610
	      compareRefMaps(iOrient, iConv, amp, pha, env, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
611

612
	      double compTime=0.;
613
	      if (DebugOutput >= 2)
614
		{
615
		  compTime = timer.GetCurrentElapsedTime();
616
617
618
619
620
621
622
		  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;

623
624
		  if (Autotuning) printf("\t\tTime Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s, with GPU workload %d%%) (rank %d)\n", iOrient, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000., aut.Workload(), mpi_rank);
		  else 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);
625
		}
626
	      if (Autotuning && aut.Needed(iConv))
627
628
		{
		  if (compTime == 0.) compTime = timer.GetCurrentElapsedTime();
629
630
631
		  aut.Tune(compTime);
		  if (aut.Finished() && DebugOutput >= 2) printf("\t\tOptimal GPU workload %d%% (rank %d)\n", aut.Workload(), mpi_rank);
		  rebalance(aut.Workload());
632
		}
633
634
635
636
	    }
	  if (DebugOutput >= 1)
	    {
	      printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient, timer2.GetCurrentElapsedTime(), mpi_rank);
637
638
	      timer2.ResetStart();
	    }
639
	}
640
641
642
643
644
    }
  //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
645

646
  deviceFinishRun();
647
648


649
650

  // ************* Collecing all the probabilities from MPI replicas ***************
David Rohr's avatar
David Rohr committed
651

David Rohr's avatar
David Rohr committed
652
#ifdef WITH_MPI
653
654
655
656
657
658
659
660
661
662
663
664
665
  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);
666

667
668
669
	for (int i = 0;i < RefMap.ntotRefMap;i++)
	  {
	    bioem_Probability_map& pProbMap = pProb.getProbMap(i);
Pilar Cossio's avatar
Pilar Cossio committed
670
#ifdef DEBUG
671
672
	    cout << "Reduction " << mpi_rank << " Map " << i << " Prob " << pProbMap.Total << " Const " << pProbMap.Constoadd  << "\n";     
#endif
673
	    tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
674

675
676
677
678
679
	  }
	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
680
	{
681
682
683
684
685
686
	  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;
687
              //temporary array that has the mpirank for the highest pProb.constant
688
689
690
691
692
	    }
	  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
693
		{
694
		  if (mpi_rank == 0) printf("Error: Could not find highest probability\n");
David Rohr's avatar
David Rohr committed
695
		}
696
	      else if (tmpi2[i] != 0) //Skip if rank 0 already has highest probability
David Rohr's avatar
David Rohr committed
697
		{
698
699
700
701
702
703
704
705
		  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
706
		}
707
708
709
	    }
	  delete[] tmpi1;
	  delete[] tmpi2;
710
	}
711

David Rohr's avatar
David Rohr committed
712
	if (mpi_rank == 0)
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
	  {
	    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)
730
	{
731
732
733
734
735
736
	  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++)
	    {
737
	      for (int j = 0;j < param.nTotGridAngles;j++)
738
                {
739
740
741
742
		  //	      tmp1[i] = pProb.getProbMap(i).Constoadd;
		  //	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
		  tmp1[i * param.nTotGridAngles + j]= pProb.getProbAngle(i, j).ConstAngle;
		}
743
	    }
744

745
746
747
748
	  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
749
		{
750
751
		  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
752
		}
753
754
755
756
757
	    }
	  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++)
758
		{
759
760
761
762
763
764
		  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];
		    }
765
		}
766
767
768
769
770
771
772
	    }
	  delete[] tmp1;
	  delete[] tmp2;
	  delete[] tmp3;
	}
    }
#endif
David Rohr's avatar
David Rohr committed
773

774
775

  // ************* Writing Out Probabilities ***************
776
777
  if (mpi_rank == 0)
    {
778
 
779
780
      // Output for Angle Probability File
      ofstream angProbfile;
781
782
783
      if(param.param_device.writeAngles)
	{
	  angProbfile.open ("ANG_PROB");
784
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
785
786
          if(!param.doquater){ angProbfile <<" RefMap:  MapNumber ; alpha[rad] - beta[rad] - gamma[rad] - logP - cal log Probability + Constant: Numerical Const.+ log (volume) + prior ang\n" ;}
	  else { angProbfile <<" RefMap:  MapNumber ; q1 - q2 -q3 - logP- cal log Probability + Constant: Numerical Const. + log (volume) + prior ang\n" ;};
787
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
788
789
	  //          angProbfile <<"Model Used: " << modelfile.c_str() << "\n";
	  //          angProbfile <<"Input Used: " << infile.c_str() << "\n";
790
	}
791
      // Output for Cross Correlation File
792
793
794
795
      ofstream ccProbfile;
      if(param.param_device.writeCC)
	{
	  ccProbfile.open ("CROSS_CORRELATION");
796
	  ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
797
798
799
          ccProbfile <<" RefMap:  MapNumber ; Pixel x - Pixel y - Cross-Correlation \n";
          ccProbfile <<"Note that the highest Cross-correlation is the best.\n";
          ccProbfile <<"If the particles are flipped, include the keyward FLIPPED in the Param file.\n";
800
          ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
801
802
	}

803
804
805
806
807
808
809
810
811
      // Output for Standard Probability
      ofstream outputProbFile;
      if(!yesoutfilename)OutfileName="Output_Probabilities";
      outputProbFile.open (OutfileName.c_str());
      outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n";   
      outputProbFile << "Notation= RefMap:  MapNumber ; LogProb natural logarithm of posterior Probability ; Constant: Numerical Const. for adding Probabilities \n";
      if(!param.doquater){
	if(param.usepsf){
	  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{
812
	  outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - alpha[rad] - beta[rad] - gamma[rad] - CTF amp - CTF defocus - CTF B-Env - center x - center y - normalization - offsett \n";}
813
814
815
      }else { 
	if(param.usepsf){
	  //     if( localcc[rx * param.param_device.NumberPixels + ry] <
Pilar Cossio's avatar
Pilar Cossio committed
816
	  outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - q4 -PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";
817
	}else{
Pilar Cossio's avatar
Pilar Cossio committed
818
          outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - q4 - CTF amp - CTF defocus - CTF B-Env - center x - center y - normalization - offsett \n";
819
        }}
820
      if(param.writeCTF) outputProbFile << " RefMap:  MapNumber ; CTFMaxParm: defocus - b-Env (B ref. Penzeck 2010)\n";
Pilar Cossio's avatar
Pilar Cossio committed
821
      if(param.yespriorAngles) outputProbFile << "**** Remark: Using Prior Proability in Angles ****\n";
822
823
      outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";

Pilar Cossio's avatar
Pilar Cossio committed
824
       
825
      // Loop over reference maps
826
      // ************* Over all maps ***************
827

828
829
830
831
832
      for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
	  // **** Total Probability ***
	  bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

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

837
838
	    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";
	    outputProbFile << "RefMap: " << iRefMap << " Maximizing Param: ";
839
840
841
            // *** Param that maximize probability****
            outputProbFile << (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)) << " ";

842

843
	  }else{ 
844
845
846
847
848
849
850
851
852
	    outputProbFile << "Warining! with Map " << iRefMap << "Numerical Integrated Probability without constant = 0.0;\n";
	    outputProbFile << "Warining RefMap: " << iRefMap << "Check that constant is finite: " << pProbMap.Constoadd  << "\n"; 
	    outputProbFile << "Warining RefMap: i) check model, ii) check refmap , iii) check GPU on/off command inconsitency\n";
	    //	    outputProbFile << "Warning! " << iRefMap << " LogProb:  "  << 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";
	  }
	  //	    outputProbFile << "RefMap: " << iRefMap << " Maximizing Param: ";

	  // *** 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
853
854
855
856
857
858

	  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] << " [] ";
	  if(param.doquater)outputProbFile << param.angles[pProbMap.max.max_prob_orient].quat4 << " [] "; 
	  outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[0] << " [] ";
859
860
861
862
863
864
	  if(!param.usepsf){outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/ 2.f /M_PI / param.elecwavel * 0.0001 << " [micro-m] ";
	  }else{outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << " [1/A²] ";}
	  if(!param.usepsf){outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " [A²] ";}
	  else{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] " ;
865
866
	  if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_norm << " [] " ;}else{outputProbFile << "N.A." << " [] ";}
	  if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_mu << " [] ";}else{outputProbFile << "N.A." << " [] ";} 
867
	  outputProbFile << "\n";
868
869
870
871
872
873
874
875
876

	  // Writing out CTF parameters if requiered
	  if(param.writeCTF && param.usepsf){

	    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*0.0001 << " [micro-m] "; 
877
	    outputProbFile << 4*M_PI*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi << " [A²] \n";
878
	  }
879

880
	  //*************** Writing Individual Angle probabilities
881
882
883
	  if(param.param_device.writeAngles)
	    {
	      for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
David Rohr's avatar
David Rohr committed
884
		{
885
		  bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
886

887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
		  myfloat_t logp=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);
		  if(!param.doquater){
		    // For Euler Angles
		    if(param.yespriorAngles){
		      logp+=param.angprior[iOrient];
		      angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << logp << " Separated: "
				  << 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) << " " << param.angprior[iOrient] << "\n";
		    } else
		      {
			angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " <<  logp << " Separated: "<<
			  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";
		      }
		  }else {
		    // Samething but for Quaternions
		    if(param.yespriorAngles){
		      logp+=param.angprior[iOrient];
		      angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << param.angles[iOrient].quat4 << " " << logp << " Separated: " << 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) << " " << param.angprior[iOrient] << "\n";
		    } else
		      {
			angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << param.angles[iOrient].quat4 << " " << logp << " Separated: "<<
			  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";
		      }
		  }
910
911
		}
	    }
912
	
913
	  //************* Writing Cross-Correlations if requiered
914
          //************* This is currently not in the manual *****
915
	  if(param.param_device.writeCC){
916

917
918
919
920
921
	    int  cc=0;
	    int halfPix;
	    int rx=0;
	    int ry=0;
	    myfloat_t localcc[ (param.param_device.NumberPixels+1) * (param.param_device.NumberPixels+1) ];
Pilar Cossio's avatar
Pilar Cossio committed
922
            int used[(param.param_device.NumberPixels+1) * (param.param_device.NumberPixels+1)];
923

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

927
928
929
930
931
	    for (int rx = 0; rx < param.param_device.NumberPixels ; rx++)
	      {
		for (int ry = 0; ry < param.param_device.NumberPixels ; ry++)
		  {
		    localcc[ rx * param.param_device.NumberPixels + ry ] = 0.0;
Pilar Cossio's avatar
Pilar Cossio committed
932
			used[ rx * param.param_device.NumberPixels + ry ]= 0;
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
959
960
961
962
963
964
965
		  }
	      }

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

		    // Applying Periodic boundary conditions to the CC
		    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";
		    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;
		    }
Pilar Cossio's avatar
Pilar Cossio committed
966
 			used[ rx * param.param_device.NumberPixels + ry] = 1;
967
968
969
970
971
		    cc++;
		  }
		//              ccProbfile << "\n";
	      }
	    if(!param.ignoreCCoff){
Pilar Cossio's avatar
Pilar Cossio committed
972
/*	      for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels ; rx = rx + param.param_device.CCdisplace)
David Rohr's avatar
David Rohr committed
973
		{
974
		  for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry = ry + param.param_device.CCdisplace)
Pilar Cossio's avatar
Pilar Cossio committed
975
976
977
978
979
980
981
		    {*/
  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(used[ rx * param.param_device.NumberPixels + ry ] == 1){
982
			ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
983
		      }else{
Pilar Cossio's avatar
Pilar Cossio committed
984
			if(localcc[ rx * param.param_device.NumberPixels + ry ] <= -FLT_MAX)ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << -FLT_MAX << "\n" ;
985
		      }
986
		      //				 cout << " cc " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] <<"\n" ;
987
		    }
Pilar Cossio's avatar
Pilar Cossio committed
988
	//	  ccProbfile << "\n";
989
990
991
		}			
	    }else{
	      for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels ; rx++)
992
		{
993
		  for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry++)
994
		    {
Pilar Cossio's avatar
Pilar Cossio committed
995
                         if(used[ rx * param.param_device.NumberPixels + ry ] == 1){
996
                        ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
997
		      }else{
Pilar Cossio's avatar
Pilar Cossio committed
998
                        if(localcc[ rx * param.param_device.NumberPixels + ry ] <= -FLT_MAX)ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << -FLT_MAX << "\n" ;
999
		      }
1000
		    }
Pilar Cossio's avatar
Pilar Cossio committed
1001
	//	  ccProbfile << "\n";
1002
		}
1003

1004
	    }
1005
	  }
1006
	}
1007

1008
1009
1010
1011
1012
1013
1014
1015
      if(param.param_device.writeAngles)
	{
	  angProbfile.close();
	}

      if(param.param_device.writeCC)
	{
	  ccProbfile.close();
1016
	}
1017

1018
1019
1020
1021
      outputProbFile.close();
    }

  return(0);
1022
1023
}

1024
int bioem::compareRefMaps(int iOrient, int iConv,  myfloat_t amp, myfloat_t pha, myfloat_t env, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
1025
{
1026

1027
1028
  //***************************************************************************************
  //***** BioEM routine for comparing reference maps to convoluted maps *****
1029
  //***************************************************************************************
1030
1031
1032
1033
1034
  if (FFTAlgo)
    {
      //With FFT Algorithm
#pragma omp parallel for schedule(dynamic, 1)
      for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
1035
	{
1036
	  const int num = omp_get_thread_num();
1037
	  calculateCCFFT(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
1038
	}
1039
1040
1041
1042
1043
1044
    }
  else
    {
      //Without FFT Algorithm
#pragma omp parallel for schedule(dynamic, 1)
      for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
1045
	{
1046
	  compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, conv_map, pProb, param.param_device, RefMap);