bioem.cpp 57.2 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.
5
   Max Planck Institute of Biophysics, Frankfurt, Germany.
6

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

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
Pilar Cossio's avatar
Pilar Cossio committed
10
#ifdef WITH_MPI
11
12
#include <mpi.h>

13
14
15
16
17
#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
18
#endif
19

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

33
#ifdef WITH_OPENMP
34
#include <omp.h>
35
#endif
36
37
38

#include <fftw3.h>
#include <math.h>
39
#include "timer.h"
40
41
42
43
44

#include "param.h"
#include "bioem.h"
#include "model.h"
#include "map.h"
45
#include "MersenneTwister.h"
46

47
48
49
50
51
52
#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]);

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

71
72
73
74
#include "bioem_algorithm.h"

using namespace boost;
namespace po = boost::program_options;
Pilar Cossio's avatar
Pilar Cossio committed
75
namespace bran= boost::random;
76
77
78

using namespace std;

79
/* For dvl nodes in hydra with problem in boost
80
81
82
   namespace std {
   typedef decltype(nullptr) nullptr_t;
   }*/
83

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

bioem::bioem()
{
94
95
96
  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"));
97
98
99
100
101
102
103
104
}

bioem::~bioem()
{
}

int bioem::configure(int ac, char* av[])
{
105
106
107
108
109
110
111
112
  // **************************************************************************************
  // **** 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
113
114


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

      // *************************************************************************************
      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
139
	  ("Particlesfile", po::value< std::string>(), "if BioEM (Mandatory) Name of particle-image file")
140
	  ("Inputfile", po::value<std::string>(), "if BioEM (Mandatory) Name of input parameter file") 
Pilar Cossio's avatar
Pilar Cossio committed
141
142
	  ("PrintBestCalMap", po::value< std::string>(), "(Optional) Only print best calculated map. NO BioEM (!)")
	  ("ReadOrientation", po::value< std::string>(), "(Optional) Read file name containing orientations")
143
144
145
	  ("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
146
147
	  ("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
148
	  ("OutputFile",  po::value< std::string>(), "(Optional) For changing the outputfile name")
149
150
151
152
153
154
155
156
157
158
159
	  ("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);
160
	p.add("ReadOrientation",-1);
161
162
163
	p.add("PrintBestCalMap",-1);
	p.add("DumpMaps", -1);
	p.add("LoadMapDump", -1);
Pilar Cossio's avatar
Pilar Cossio committed
164
        p.add("OutputFile",-1);
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179

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

181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
	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;
	  }
198
	if (vm.count("ReadOrientation"))
199
	  {
200
201
202
203
204
205
206
207
	    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  >();
208
209
	    param.notuniformangles=true;
	  }
Pilar Cossio's avatar
Pilar Cossio committed
210
211
	if (vm.count("OutputFile"))
          {
212
213
214
215
	    OutfileName = vm["OutputFile"].as< std::string >();
	    cout << "Writing OUTPUT to: " <<  vm["OutputFile"].as<  std::string  >() << "\n";
	    yesoutfilename=true;
	  }
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
	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
260

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

268
269
270
271
      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
272

273
274
275
276
      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
277
	param.readParameters(infile.c_str());
David Rohr's avatar
David Rohr committed
278

279
280
	// ********************* Reading Particle Maps Input **********************
	RefMap.readRefMaps(param, mapfile.c_str());
David Rohr's avatar
David Rohr committed
281

282

283
284
      } else{
	// Reading parameters for only writting down Best projection
285

286
287
	param.forprintBest(Inputbestmap.c_str());
      }	
288

289
      // ********************* Reading Model Input ******************************
290
      Model.readModel(param, modelfile.c_str());
David Rohr's avatar
David Rohr committed
291

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

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

301
302
      // Generating Grids of orientations 
      if(!param.printModel)param.CalculateGridsParam(Inputanglefile.c_str());
303
    }
304

305
#ifdef WITH_MPI
306

Pilar Cossio's avatar
Pilar Cossio committed
307

308
309

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

317
318
      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);
319

Pilar Cossio's avatar
Pilar Cossio committed
320
#ifdef DEBUG
321
322
      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
323
324

#endif
325
      //****refCtf, CtfParam, angles automatically filled by precalculate function bellow
326
327
328
329
330
331
332
333
334

      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
335

336
337
    }
#endif
338

339
340
341
  // ****************** Precalculating Necessary Stuff *********************
  if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
  param.PrepareFFTs();
342

343
344
345
346
347
348
349
  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Prepare FFTs %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
  precalculate();

350
  // ****************** For debugging *********************
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
  if (getenv("BIOEM_DEBUG_BREAK"))
    {
      const int cut = atoi(getenv("BIOEM_DEBUG_BREAK"));
      if (param.nTotGridAngles > cut) param.nTotGridAngles = cut;
      if (param.nTotCTFs > cut) param.nTotCTFs = cut;
    }

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

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

  // ****************** Initializng pointers *********************

373
  deviceInit();
374

375
376
377
  if (DebugOutput >= 2 && mpi_rank == 0) printf("Time Device Init %f\n", timer.GetCurrentElapsedTime());

  return(0);
378
379
}

380
381
void bioem::cleanup()
{
382
  //Deleting allocated pointers
383
384
  free_device_host(pProb.ptr);
  RefMap.freePointers();
385
386
}

387
388
int bioem::precalculate()
{
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
  // **************************************************************************************
  // **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);
411
412
413
414
}

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

416
417
418
  // **************************************************************************************
  // ********** Secondary routine for printing out the only best projection ***************
  // **************************************************************************************
David Rohr's avatar
David Rohr committed
419

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

422
423
424
425
426
    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;
427

428
429
430
    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);
431

432
    cout << "...... Calculating Projection .......................\n " ;
433

434
    createProjection(0, proj_mapsFFT);
435

436
    cout << "...... Calculating Convolution .......................\n " ;
David Rohr's avatar
David Rohr committed
437

438
    createConvolutedProjectionMap(0, 0, proj_mapsFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
439

440
  }
441

442
443
444
  // **************************************************************************************
  // **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
  // **************************************************************************************
445

446
447
  // **** 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
448

449
  if (mpi_rank == 0) printf("\tInitializing Probabilities\n");
450

Pilar Cossio's avatar
Pilar Cossio committed
451
452
453
454
  // Contros for MPI
  if(mpi_size > param.nTotGridAngles){
    cout << "EXIT: Wrong MPI setup More MPI processes than orientations\n"; exit(1);
  }
455

456
457
458
459
  // 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
460

461
      pProbMap.Total = 0.0;
462
463
      pProbMap.Constoadd = -FLT_MAX; //Problem if using double presicion

464
      if (param.param_device.writeAngles)
465
	{
466
467
468
469
470
	  for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
	    {
	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);

	      pProbAngle.forAngles = 0.0;
471
	      pProbAngle.ConstAngle = -FLT_MAX;
472
473
	    }
	}
474

475
      if (param.param_device.writeCC)
476
	{      int  cc=0;
477
	  for (int cent_x = 0; cent_x < param.param_device.NumberPixels; cent_x = cent_x + param.param_device.CCdisplace)
478
	    {
479
	      for (int cent_y = 0; cent_y < param.param_device.NumberPixels; cent_y = cent_y + param.param_device.CCdisplace)
480
		{
481
		  bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
482
483
484
485
486
487
488
489
490
		  //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++;
491
		}
492
	    }
493
	  if(!FFTAlgo){cout << "Cross correlation calculation must be with enviormental variable FFTALGO=1\n"; exit(1);}
494
495
	}                 
    }
496
497
498

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

499
  // **************************************************************************************
500

501
502
503
504
505
506
507
508
509
510
511
  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
512
  //******** Alocating Vectors *************
513
514
515
516
517
518
519
520
521
  proj_mapsFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * ProjMapSize * nProjectionsAtOnce);
  conv_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
  if (!FFTAlgo) conv_map = (myfloat_t*) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
             

  HighResTimer timer, timer2;

  if (DebugOutput >= 1 && mpi_rank == 0) printf("\tMain Loop GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)², Pixels %d², OMP Threads %d, MPI Ranks %d\n", param.nTotGridAngles, param.nTotCTFs, RefMap.ntotRefMap, 2 * param.param_device.maxDisplaceCenter + param.param_device.GridSpaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, omp_get_max_threads(), mpi_size);

522
523


524
525
526
527
  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;

528
529
530

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

531
532
533
534
535
536
537
  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);
538
539
540

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

541
542
543
544
545
546
#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);
547

548
549
550
      for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
	{
	  mycomplex_t* proj_mapFFT = &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize];
551

552
	  // ***************************************************************************************
553
554
	  // ***** **** Internal Loop over PSF/CTF convolutions **** *****

555
556
557
558
559
560
561
562
	  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);

563
	  
564
	      if (DebugOutput >= 2) timer.ResetStart();
565
566
567
568
569
570
     	      myfloat_t amp,pha,env;

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

571
572
573
	      // ******************Internal loop over Reference images CUDA or OpenMP******************
	      // *** Comparing each calculated convoluted map with all experimental maps ***

574
	      compareRefMaps(iOrient, iConv, amp, pha, env, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
575
576

	      if (DebugOutput >= 2)
577
		{
578
579
580
581
582
583
584
585
586
		  const double compTime = timer.GetCurrentElapsedTime();
		  const int nShifts = 2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1;
		  const double nFlops = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
		    (((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 5. + 25.) / compTime;
		  const double nGBs = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
		    (((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 2. + 8.) * (double) sizeof(myfloat_t) / compTime;
		  const double nGBs2 = (double) RefMap.ntotRefMap * ((double) param.param_device.NumberPixels * (double) param.param_device.NumberPixels + 8.) * (double) sizeof(myfloat_t) / compTime;

		  printf("\t\tTime Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s) (rank %d)\n", iOrient, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000., mpi_rank);
587
		}
588
589
590
591
592
593
	    }
	  if (DebugOutput >= 1)
	    {
	      printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient, timer2.GetCurrentElapsedTime(), mpi_rank);
	      timer2.ResetStart();
	    }
594
	}
595
596
597
598
599
    }
  //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
600

601
  deviceFinishRun();
602
603


604
605

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

David Rohr's avatar
David Rohr committed
607
#ifdef WITH_MPI
608
609
610
611
612
613
614
615
616
617
618
619
620
  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);
621

622
623
624
	for (int i = 0;i < RefMap.ntotRefMap;i++)
	  {
	    bioem_Probability_map& pProbMap = pProb.getProbMap(i);
Pilar Cossio's avatar
Pilar Cossio committed
625
#ifdef DEBUG
626
627
	    cout << "Reduction " << mpi_rank << " Map " << i << " Prob " << pProbMap.Total << " Const " << pProbMap.Constoadd  << "\n";     
#endif
628
	    tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
629

630
631
632
633
634
	  }
	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
635
	{
636
637
638
639
640
641
	  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;
642
              //temporary array that has the mpirank for the highest pProb.constant
643
644
645
646
647
	    }
	  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
648
		{
649
		  if (mpi_rank == 0) printf("Error: Could not find highest probability\n");
David Rohr's avatar
David Rohr committed
650
		}
651
	      else if (tmpi2[i] != 0) //Skip if rank 0 already has highest probability
David Rohr's avatar
David Rohr committed
652
		{
653
654
655
656
657
658
659
660
		  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
661
		}
662
663
664
	    }
	  delete[] tmpi1;
	  delete[] tmpi2;
665
	}
666

David Rohr's avatar
David Rohr committed
667
	if (mpi_rank == 0)
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
	  {
	    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)
685
	{
686
687
688
689
690
691
	  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++)
	    {
692
	      for (int j = 0;j < param.nTotGridAngles;j++)
693
                {
694
695
696
697
		  //	      tmp1[i] = pProb.getProbMap(i).Constoadd;
		  //	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
		  tmp1[i * param.nTotGridAngles + j]= pProb.getProbAngle(i, j).ConstAngle;
		}
698
	    }
699

700
701
702
703
	  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
704
		{
705
706
		  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
707
		}
708
709
710
711
712
	    }
	  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++)
713
		{
714
715
716
717
718
719
		  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];
		    }
720
		}
721
722
723
724
725
726
727
	    }
	  delete[] tmp1;
	  delete[] tmp2;
	  delete[] tmp3;
	}
    }
#endif
David Rohr's avatar
David Rohr committed
728

729
730

  // ************* Writing Out Probabilities ***************
731
732
  if (mpi_rank == 0)
    {
733
 
734
735
      // Output for Angle Probability File
      ofstream angProbfile;
736
737
738
      if(param.param_device.writeAngles)
	{
	  angProbfile.open ("ANG_PROB");
739
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
740
741
          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" ;};
742
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
743
744
	  //          angProbfile <<"Model Used: " << modelfile.c_str() << "\n";
	  //          angProbfile <<"Input Used: " << infile.c_str() << "\n";
745
	}
746
      // Output for Cross Correlation File
747
748
749
750
      ofstream ccProbfile;
      if(param.param_device.writeCC)
	{
	  ccProbfile.open ("CROSS_CORRELATION");
751
	  ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
752
753
754
          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";
755
          ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
756
757
	}

758
759
760
761
762
763
764
765
766
      // 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{
767
	  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";}
768
769
770
      }else { 
	if(param.usepsf){
	  //     if( localcc[rx * param.param_device.NumberPixels + ry] <
Pilar Cossio's avatar
Pilar Cossio committed
771
	  outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - q4 -PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";
772
	}else{
Pilar Cossio's avatar
Pilar Cossio committed
773
          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";
774
        }}
775
      if(param.writeCTF) outputProbFile << " RefMap:  MapNumber ; CTFMaxParm: defocus - b-Env (B ref. Penzeck 2010)\n";
Pilar Cossio's avatar
Pilar Cossio committed
776
      if(param.yespriorAngles) outputProbFile << "**** Remark: Using Prior Proability in Angles ****\n";
777
778
      outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";

Pilar Cossio's avatar
Pilar Cossio committed
779
       
780
      // Loop over reference maps
781
      // ************* Over all maps ***************
782

783
784
785
786
787
      for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
	  // **** Total Probability ***
	  bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

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

792
793
	    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: ";
794
795
796
            // *** 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)) << " ";

797

798
	  }else{ 
799
800
801
802
803
804
805
806
807
	    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
808
809
810
811
812
813

	  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] << " [] ";
814
815
816
817
818
819
	  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] " ;
820
821
	  if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_norm << " [] " ;}else{outputProbFile << "N.A." << " [] ";}
	  if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_mu << " [] ";}else{outputProbFile << "N.A." << " [] ";} 
822
	  outputProbFile << "\n";
823
824
825
826
827
828
829
830
831

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

835
	  //*************** Writing Individual Angle probabilities
836
837
838
	  if(param.param_device.writeAngles)
	    {
	      for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
David Rohr's avatar
David Rohr committed
839
		{
840
		  bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
841

842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
		  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";
		      }
		  }
865
866
		}
	    }
867
	
868
	  //************* Writing Cross-Correlations if requiered
869
          //************* This is currently not in the manual *****
870
	  if(param.param_device.writeCC){
871

872
873
874
875
876
	    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
877
            int used[(param.param_device.NumberPixels+1) * (param.param_device.NumberPixels+1)];
878

879
880
	    halfPix = param.param_device.NumberPixels / 2 ;
	    // Ordering the centers of the Cross Correlation
881

882
883
884
885
886
	    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
887
			used[ rx * param.param_device.NumberPixels + ry ]= 0;
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
		  }
	      }

	    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
921
 			used[ rx * param.param_device.NumberPixels + ry] = 1;
922
923
924
925
926
		    cc++;
		  }
		//              ccProbfile << "\n";
	      }
	    if(!param.ignoreCCoff){
Pilar Cossio's avatar
Pilar Cossio committed
927
/*	      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
928
		{
929
		  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
930
931
932
933
934
935
936
		    {*/
  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){
937
			ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
938
		      }else{
Pilar Cossio's avatar
Pilar Cossio committed
939
			if(localcc[ rx * param.param_device.NumberPixels + ry ] <= -FLT_MAX)ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << -FLT_MAX << "\n" ;
940
		      }
941
		      //				 cout << " cc " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] <<"\n" ;
942
		    }
Pilar Cossio's avatar
Pilar Cossio committed
943
	//	  ccProbfile << "\n";
944
945
946
		}			
	    }else{
	      for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels ; rx++)
947
		{
948
		  for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry++)
949
		    {
Pilar Cossio's avatar
Pilar Cossio committed
950
                         if(used[ rx * param.param_device.NumberPixels + ry ] == 1){
951
                        ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
952
		      }else{
Pilar Cossio's avatar
Pilar Cossio committed
953
                        if(localcc[ rx * param.param_device.NumberPixels + ry ] <= -FLT_MAX)ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << -FLT_MAX << "\n" ;
954
		      }
955
		    }
Pilar Cossio's avatar
Pilar Cossio committed
956
	//	  ccProbfile << "\n";
957
		}
958

959
	    }
960
	  }
961
	}
962

963
964
965
966
967
968
969
970
      if(param.param_device.writeAngles)
	{
	  angProbfile.close();
	}

      if(param.param_device.writeCC)
	{
	  ccProbfile.close();
971
	}
972

973
974
975
976
      outputProbFile.close();
    }

  return(0);
977
978
}

979
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)
980
{
981

982
983
  //***************************************************************************************
  //***** BioEM routine for comparing reference maps to convoluted maps *****
984
  //***************************************************************************************
985
986
987
988
989
  if (FFTAlgo)
    {
      //With FFT Algorithm
#pragma omp parallel for schedule(dynamic, 1)
      for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
990
	{
991
	  const int num = omp_get_thread_num();
992
	  calculateCCFFT(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
993
	}
994
995
996
997
998
999
    }
  else
    {
      //Without FFT Algorithm
#pragma omp parallel for schedule(dynamic, 1)
      for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
1000
	{
1001
	  compareRefMapShifted <