bioem.cpp 59.8 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"
Luka Stanisic's avatar
Luka Stanisic committed
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
#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]);
56
enum myColor { COLOR_PROJECTION, COLOR_CONVOLUTION, COLOR_COMPARISON, COLOR_WORKLOAD };
57

58
59
60
61
// Projection number is stored in category attribute
// Convolution number is stored in payload attribute

#define cuda_custom_timeslot(name,iMap,iConv,cid) {	\
62
63
64
65
66
67
68
    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];		\
69
70
71
    eventAttrib.category = iMap;			\
    eventAttrib.payloadType = NVTX_PAYLOAD_TYPE_UNSIGNED_INT64; \
    eventAttrib.payload.llValue = iConv;		\
72
73
74
75
    eventAttrib.messageType = NVTX_MESSAGE_TYPE_ASCII;	\
    eventAttrib.message.ascii = name;			\
    nvtxRangePushEx(&eventAttrib);			\
  }
76
77
#define cuda_custom_timeslot_end nvtxRangePop();
#else
78
#define cuda_custom_timeslot(name,iMap,iConv,cid)
79
80
#define cuda_custom_timeslot_end
#endif
81

82
83
84
85
#include "bioem_algorithm.h"

using namespace boost;
namespace po = boost::program_options;
Pilar Cossio's avatar
Pilar Cossio committed
86
namespace bran= boost::random;
87
88
89

using namespace std;

90
/* For dvl nodes in hydra with problem in boost
91
92
93
   namespace std {
   typedef decltype(nullptr) nullptr_t;
   }*/
94

95
96
97
98
// A helper function of Boost
template<class T>
ostream& operator<<(ostream& os, const vector<T>& v)
{
99
100
  copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
  return os;
101
102
103
104
}

bioem::bioem()
{
105
106
107
  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"));
Luka Stanisic's avatar
Luka Stanisic committed
108
  Autotuning = false;
109
110
111
112
113
114
115
116
}

bioem::~bioem()
{
}

int bioem::configure(int ac, char* av[])
{
117
118
119
120
121
122
123
124
  // **************************************************************************************
  // **** 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
125
126


127
  std::string infile, modelfile, mapfile,Inputanglefile,Inputbestmap; 
128
129
130
131
132
133
134
135
136
137
138
  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
139
      yesoutfilename=false;
140
141
142
143
144
145
146
147
148
149
150

      // *************************************************************************************
      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
151
	  ("Particlesfile", po::value< std::string>(), "if BioEM (Mandatory) Name of particle-image file")
152
	  ("Inputfile", po::value<std::string>(), "if BioEM (Mandatory) Name of input parameter file") 
Pilar Cossio's avatar
Pilar Cossio committed
153
154
	  ("PrintBestCalMap", po::value< std::string>(), "(Optional) Only print best calculated map. NO BioEM (!)")
	  ("ReadOrientation", po::value< std::string>(), "(Optional) Read file name containing orientations")
155
156
157
	  ("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
158
159
	  ("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
160
	  ("OutputFile",  po::value< std::string>(), "(Optional) For changing the outputfile name")
161
162
163
164
165
166
167
168
169
170
171
	  ("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);
172
	p.add("ReadOrientation",-1);
173
174
175
	p.add("PrintBestCalMap",-1);
	p.add("DumpMaps", -1);
	p.add("LoadMapDump", -1);
Pilar Cossio's avatar
Pilar Cossio committed
176
        p.add("OutputFile",-1);
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191

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

193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
	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;
	  }
210
	if (vm.count("ReadOrientation"))
211
	  {
212
213
214
215
216
217
218
219
	    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  >();
220
221
	    param.notuniformangles=true;
	  }
Pilar Cossio's avatar
Pilar Cossio committed
222
223
	if (vm.count("OutputFile"))
          {
224
225
226
227
	    OutfileName = vm["OutputFile"].as< std::string >();
	    cout << "Writing OUTPUT to: " <<  vm["OutputFile"].as<  std::string  >() << "\n";
	    yesoutfilename=true;
	  }
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
265
266
267
268
269
270
271
	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
272

273
274
      //check for consitency in multiple MRCs
      if(RefMap.readMultMRC && not(RefMap.readMRC))
David Rohr's avatar
David Rohr committed
275
	{
276
277
278
	  cout << "For Multiple MRCs command --ReadMRC is necesary too";
	  exit(1);
	}
David Rohr's avatar
David Rohr committed
279

280
281
282
283
      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
284

285
286
287
288
      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
289
	param.readParameters(infile.c_str());
David Rohr's avatar
David Rohr committed
290

291
292
	// ********************* Reading Particle Maps Input **********************
	RefMap.readRefMaps(param, mapfile.c_str());
David Rohr's avatar
David Rohr committed
293

294

295
296
      } else{
	// Reading parameters for only writting down Best projection
297

298
299
	param.forprintBest(Inputbestmap.c_str());
      }	
300

301
      // ********************* Reading Model Input ******************************
302
      Model.readModel(param, modelfile.c_str());
David Rohr's avatar
David Rohr committed
303

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

306
      if (DebugOutput >= 2 && mpi_rank == 0) printf("Reading Input Data Time: %f\n", timer.GetCurrentElapsedTime());
307
    
308
      if(param.param_device.writeCC && mpi_size>1){
309
310
	cout << "Exiting::: WRITE CROSS-CORRELATION ONLY VAILD FOR 1 MPI PROCESS\n";
        exit(1);
311
      }
Pilar Cossio's avatar
Pilar Cossio committed
312

313
314
      // Generating Grids of orientations 
      if(!param.printModel)param.CalculateGridsParam(Inputanglefile.c_str());
315
    }
316

317
#ifdef WITH_MPI
318

Pilar Cossio's avatar
Pilar Cossio committed
319

320
321

  // ********************* MPI inizialization/ Transfer of parameters******************
322
323
324
325
  if (mpi_size > 1)
    {
      if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
      MPI_Bcast(&param, sizeof(param), MPI_BYTE, 0, MPI_COMM_WORLD);
326
      //We have to reinitialize all pointers !!!!!!!!!!!!
Pilar Cossio's avatar
Pilar Cossio committed
327
      if (mpi_rank != 0) param.angprior = NULL;
Pilar Cossio's avatar
Pilar Cossio committed
328

329
330
      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);
331

Pilar Cossio's avatar
Pilar Cossio committed
332
#ifdef DEBUG
333
334
      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
335
336

#endif
337
      //****refCtf, CtfParam, angles automatically filled by precalculate function bellow
338
339
340
341
342
343
344
345
346

      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
347

348
349
    }
#endif
350

351
352
353
  // ****************** Precalculating Necessary Stuff *********************
  if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
  param.PrepareFFTs();
354

355
356
357
358
359
360
361
  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Prepare FFTs %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
  precalculate();

362
  // ****************** For debugging *********************
363
364
365
366
367
368
369
370
371
372
373
374
  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
375
376
377
378
379
380
381
382
383
384
385
386
387

  // ****************** 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");
    }

388
389
390
  // ****************** Initializing pointers *********************

  deviceInit();
391
392
393

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

398
  if(!param.printModel)pProb.init(RefMap.ntotRefMap, param.nTotGridAngles, param.nTotCC, *this);
399

400
401
402
403
404
  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Init Probabilities %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
405
406

  return(0);
407
408
}

409
410
void bioem::cleanup()
{
411
  //Deleting allocated pointers
412
413
  free_device_host(pProb.ptr);
  RefMap.freePointers();
414
415
}

416
417
int bioem::precalculate()
{
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
  // **************************************************************************************
  // **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);
440
441
442
443
}

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

445
446
447
  // **************************************************************************************
  // ********** Secondary routine for printing out the only best projection ***************
  // **************************************************************************************
David Rohr's avatar
David Rohr committed
448

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

451
452
453
454
455
    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;
456

457
458
459
    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);
460

461
    cout << "...... Calculating Projection .......................\n " ;
462

463
    createProjection(0, proj_mapsFFT);
464

465
    cout << "...... Calculating Convolution .......................\n " ;
David Rohr's avatar
David Rohr committed
466

467
    createConvolutedProjectionMap(0, 0, proj_mapsFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
468

469
  }
470

471
472
473
  // **************************************************************************************
  // **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
  // **************************************************************************************
474

475
476
  // **** 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
477

478
  if (mpi_rank == 0) printf("\tInitializing Probabilities\n");
479

Pilar Cossio's avatar
Pilar Cossio committed
480
481
482
483
  // Contros for MPI
  if(mpi_size > param.nTotGridAngles){
    cout << "EXIT: Wrong MPI setup More MPI processes than orientations\n"; exit(1);
  }
484

485
486
487
488
  // 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
489

490
      pProbMap.Total = 0.0;
491
492
      pProbMap.Constoadd = -FLT_MAX; //Problem if using double presicion

493
      if (param.param_device.writeAngles)
494
	{
495
496
497
498
499
	  for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
	    {
	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);

	      pProbAngle.forAngles = 0.0;
500
	      pProbAngle.ConstAngle = -FLT_MAX;
501
502
	    }
	}
503

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

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

528
  // **************************************************************************************
529

530
531
532
533
534
535
536
537
538
539
540
  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
541
  //******** Alocating Vectors *************
542
543
544
545
546
547
548
  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;

Luka Stanisic's avatar
Luka Stanisic committed
549
550
551
552
553
  /* Autotuning */
  Autotuner aut;
  if (Autotuning)
    {
      aut.Initialize(AUTOTUNING_ALGORITHM, FIRST_STABLE);
554
      rebalanceWrapper(aut.Workload());
Luka Stanisic's avatar
Luka Stanisic committed
555
556
    }

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

559
560


561
562
563
564
  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;

565
566
567

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

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

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

578
579
580
581
582
583
#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);
584

Luka Stanisic's avatar
Luka Stanisic committed
585
586
587
588
      /* Recalibrate if needed */
      if (Autotuning && ((iOrientAtOnce - iOrientStart) % RECALIB_FACTOR == 0) && ((iOrientEnd - iOrientAtOnce) > RECALIB_FACTOR) && (iOrientAtOnce != iOrientStart))
	{
	  aut.Reset();
589
	  rebalanceWrapper(aut.Workload());
Luka Stanisic's avatar
Luka Stanisic committed
590
591
	}

592
593
594
      for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
	{
	  mycomplex_t* proj_mapFFT = &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize];
595

596
	  // ***************************************************************************************
597
598
	  // ***** **** Internal Loop over PSF/CTF convolutions **** *****

599
600
601
602
603
604
605
606
	  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);

Luka Stanisic's avatar
Luka Stanisic committed
607
      	      if ((DebugOutput >= 2) || (Autotuning && aut.Needed(iConv))) timer.ResetStart();
608
609
610
611
612
613
     	      myfloat_t amp,pha,env;

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

614
615
616
	      // ******************Internal loop over Reference images CUDA or OpenMP******************
	      // *** Comparing each calculated convoluted map with all experimental maps ***

617
	      compareRefMaps(iOrient, iConv, amp, pha, env, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
618

Luka Stanisic's avatar
Luka Stanisic committed
619
	      double compTime=0.;
620
	      if (DebugOutput >= 2)
621
		{
Luka Stanisic's avatar
Luka Stanisic committed
622
		  compTime = timer.GetCurrentElapsedTime();
623
624
625
626
627
628
629
		  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;

Luka Stanisic's avatar
Luka Stanisic committed
630
631
632
633
634
635
636
637
		  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);
		}
	      if (Autotuning && aut.Needed(iConv))
		{
		  if (compTime == 0.) compTime = timer.GetCurrentElapsedTime();
		  aut.Tune(compTime);
		  if (aut.Finished() && DebugOutput >= 2) printf("\t\tOptimal GPU workload %d%% (rank %d)\n", aut.Workload(), mpi_rank);
638
		  rebalanceWrapper(aut.Workload());
639
		}
640
641
642
643
644
645
	    }
	  if (DebugOutput >= 1)
	    {
	      printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient, timer2.GetCurrentElapsedTime(), mpi_rank);
	      timer2.ResetStart();
	    }
646
	}
647
648
649
650
651
    }
  //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
652

653
  deviceFinishRun();
654
655


656
657

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

David Rohr's avatar
David Rohr committed
659
#ifdef WITH_MPI
660
661
662
663
664
665
666
667
668
669
670
671
672
  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);
673

674
675
676
	for (int i = 0;i < RefMap.ntotRefMap;i++)
	  {
	    bioem_Probability_map& pProbMap = pProb.getProbMap(i);
Pilar Cossio's avatar
Pilar Cossio committed
677
#ifdef DEBUG
678
679
	    cout << "Reduction " << mpi_rank << " Map " << i << " Prob " << pProbMap.Total << " Const " << pProbMap.Constoadd  << "\n";     
#endif
680
	    tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
681

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

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

752
753
754
755
	  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
756
		{
757
758
		  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
759
		}
760
761
762
763
764
	    }
	  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++)
765
		{
766
767
768
769
770
771
		  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];
		    }
772
		}
773
774
775
776
777
778
779
	    }
	  delete[] tmp1;
	  delete[] tmp2;
	  delete[] tmp3;
	}
    }
#endif
David Rohr's avatar
David Rohr committed
780

781
782

  // ************* Writing Out Probabilities ***************
783
784
  if (mpi_rank == 0)
    {
785
 
786
787
      // Output for Angle Probability File
      ofstream angProbfile;
788
789
790
      if(param.param_device.writeAngles)
	{
	  angProbfile.open ("ANG_PROB");
791
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
792
793
          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" ;};
794
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
795
796
	  //          angProbfile <<"Model Used: " << modelfile.c_str() << "\n";
	  //          angProbfile <<"Input Used: " << infile.c_str() << "\n";
797
	}
798
      // Output for Cross Correlation File
799
800
801
802
      ofstream ccProbfile;
      if(param.param_device.writeCC)
	{
	  ccProbfile.open ("CROSS_CORRELATION");
803
	  ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
804
805
806
          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";
807
          ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
808
809
	}

810
811
812
813
814
815
816
817
818
      // 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{
819
	  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";}
820
821
822
      }else { 
	if(param.usepsf){
	  //     if( localcc[rx * param.param_device.NumberPixels + ry] <
Pilar Cossio's avatar
Pilar Cossio committed
823
	  outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - q4 -PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";
824
	}else{
Pilar Cossio's avatar
Pilar Cossio committed
825
          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";
826
        }}
827
      if(param.writeCTF) outputProbFile << " RefMap:  MapNumber ; CTFMaxParm: defocus - b-Env (B ref. Penzeck 2010)\n";
Pilar Cossio's avatar
Pilar Cossio committed
828
      if(param.yespriorAngles) outputProbFile << "**** Remark: Using Prior Proability in Angles ****\n";
829
830
      outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";

Pilar Cossio's avatar
Pilar Cossio committed
831
       
832
      // Loop over reference maps
833
      // ************* Over all maps ***************
834

835
836
837
838
839
      for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
	  // **** Total Probability ***
	  bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

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

844
845
	    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: ";
846
847
848
            // *** 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)) << " ";

849

850
	  }else{ 
851
852
853
854
855
856
857
858
859
	    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
860
861
862
863
864
865

	  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] << " [] ";
866
867
868
869
870
871
	  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] " ;
872
873
	  if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_norm << " [] " ;}else{outputProbFile << "N.A." << " [] ";}
	  if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_mu << " [] ";}else{outputProbFile << "N.A." << " [] ";} 
874
	  outputProbFile << "\n";
875
876
877
878
879
880
881
882
883

	  // 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] "; 
884
	    outputProbFile << 4*M_PI*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi << " [A²] \n";
885
	  }
886

887
	  //*************** Writing Individual Angle probabilities
888
889
890
	  if(param.param_device.writeAngles)
	    {
	      for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
David Rohr's avatar
David Rohr committed
891
		{
892
		  bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
893

894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
		  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";
		      }
		  }
917
918
		}
	    }
919
	
920
	  //************* Writing Cross-Correlations if requiered
921
          //************* This is currently not in the manual *****
922
	  if(param.param_device.writeCC){
923

924
925
926
927
928
	    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
929
            int used[(param.param_device.NumberPixels+1) * (param.param_device.NumberPixels+1)];
930

931
932
	    halfPix = param.param_device.NumberPixels / 2 ;
	    // Ordering the centers of the Cross Correlation
933

934
935
936
937
938
	    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
939
			used[ rx * param.param_device.NumberPixels + ry ]= 0;
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
966
967
968
969
970
971
972
		  }
	      }

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