bioem.cpp 60.5 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, COLOR_INIT };
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
  FFTAlgo = getenv("FFTALGO") == NULL ? 1 : atoi(getenv("FFTALGO"));
106
  DebugOutput = getenv("BIOEM_DEBUG_OUTPUT") == NULL ? 0 : atoi(getenv("BIOEM_DEBUG_OUTPUT"));
107
  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  *************************
477
  cuda_custom_timeslot("Initialization", -1, -1, COLOR_INIT);
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
  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);
545 546

  cuda_custom_timeslot_end; //Ending initialization
547 548 549

  HighResTimer timer, timer2;

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

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

560 561


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

566 567 568
  /* Vectors for computing statistic on different parts of the code */
  TimeStat ts((iOrientEnd - iOrientStart), param.nTotCTFs);
  if (DebugOutput >= 1) ts.InitTimeStat(4);
569 570 571

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

572 573 574 575
  for (int iOrientAtOnce = iOrientStart; iOrientAtOnce < iOrientEnd; iOrientAtOnce += nProjectionsAtOnce)
    {
      // ***************************************************************************************
      // ***** Creating Projection for given orientation and transforming to Fourier space *****
576 577 578 579 580
      if (DebugOutput >= 1)
	{
	  timer2.ResetStart();
	  timer.ResetStart();
	}
581
      int iTmpEnd = std::min(iOrientEnd, iOrientAtOnce + nProjectionsAtOnce);
582 583 584

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

585 586 587 588 589
#pragma omp parallel for
      for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
	{
	  createProjection(iOrient, &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize]);
	}
590 591 592 593 594 595
      if (DebugOutput >= 1)
	{
	  ts.time = timer.GetCurrentElapsedTime();
	  ts.Add(TS_PROJECTION);
	  if (DebugOutput >= 2) printf("\tTime Projection %d: %f (rank %d)\n", iOrientAtOnce, ts.time, mpi_rank);
	}
Luka Stanisic's avatar
Luka Stanisic committed
596 597 598 599
      /* Recalibrate if needed */
      if (Autotuning && ((iOrientAtOnce - iOrientStart) % RECALIB_FACTOR == 0) && ((iOrientEnd - iOrientAtOnce) > RECALIB_FACTOR) && (iOrientAtOnce != iOrientStart))
	{
	  aut.Reset();
600
	  rebalanceWrapper(aut.Workload());
Luka Stanisic's avatar
Luka Stanisic committed
601 602
	}

603 604 605
      for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
	{
	  mycomplex_t* proj_mapFFT = &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize];
606

607
	  // ***************************************************************************************
608 609
	  // ***** **** Internal Loop over PSF/CTF convolutions **** *****

610 611 612
	  for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
	    {
	      // *** Calculating convolutions of projection map and crosscorrelations ***
613
	      if (DebugOutput >= 1) timer.ResetStart();
614
	      createConvolutedProjectionMap(iOrient, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
615 616 617 618 619 620
	      if (DebugOutput >= 1)
		{
		  ts.time = timer.GetCurrentElapsedTime();
		  ts.Add(TS_CONVOLUTION);
		  if (DebugOutput >= 2) printf("\t\tTime Convolution %d %d: %f (rank %d)\n", iOrient, iConv, ts.time, mpi_rank);
		}
621

622 623
	      if ((DebugOutput >= 1) || (Autotuning && aut.Needed(iConv))) timer.ResetStart();
	      myfloat_t amp,pha,env;
624

625 626 627
	      amp=param.CtfParam[iConv].pos[0];
	      pha=param.CtfParam[iConv].pos[1];
	      env=param.CtfParam[iConv].pos[2];
628

629 630 631
	      // ******************Internal loop over Reference images CUDA or OpenMP******************
	      // *** Comparing each calculated convoluted map with all experimental maps ***

632
	      compareRefMaps(iOrient, iConv, amp, pha, env, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
633

634 635 636 637 638 639
	      ts.time = 0.;
	      if (DebugOutput >= 1)
		{
		  ts.time = timer.GetCurrentElapsedTime();
		  ts.Add(TS_COMPARISON);
		}
640
	      if (DebugOutput >= 2)
641
		{
642 643
		  const int nShifts = 2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1;
		  const double nFlops = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
644
		    (((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 5. + 25.) / ts.time;
645
		  const double nGBs = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
646 647
		    (((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) / ts.time;
		  const double nGBs2 = (double) RefMap.ntotRefMap * ((double) param.param_device.NumberPixels * (double) param.param_device.NumberPixels + 8.) * (double) sizeof(myfloat_t) / ts.time;
648

649 650
		  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, ts.time, 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, ts.time, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000., mpi_rank);
Luka Stanisic's avatar
Luka Stanisic committed
651 652 653
		}
	      if (Autotuning && aut.Needed(iConv))
		{
654 655 656
		  if (ts.time == 0.) ts.time = timer.GetCurrentElapsedTime();
		  aut.Tune(ts.time);
		  if (aut.Finished() && DebugOutput >= 1) printf("\tOptimal GPU workload %d%% (rank %d)\n", aut.Workload(), mpi_rank);
657
		  rebalanceWrapper(aut.Workload());
658
		}
659 660 661
	    }
	  if (DebugOutput >= 1)
	    {
662 663 664
	      ts.time = timer2.GetCurrentElapsedTime();
	      ts.Add(TS_TPROJECTION);
	      printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient, ts.time, mpi_rank);
665 666
	      timer2.ResetStart();
	    }
667
	}
668
    }
669 670 671 672 673 674 675
  /* Statistical summary on different parts of the code */
  if (DebugOutput >= 1)
    {
      ts.PrintTimeStat(mpi_rank);
      ts.EmptyTimeStat();
    }

676 677 678 679
  //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
680

681
  deviceFinishRun();
682 683


684 685

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

David Rohr's avatar
David Rohr committed
687
#ifdef WITH_MPI
688 689 690 691 692 693 694 695 696 697 698 699 700
  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);
701

702 703 704
	for (int i = 0;i < RefMap.ntotRefMap;i++)
	  {
	    bioem_Probability_map& pProbMap = pProb.getProbMap(i);
Pilar Cossio's avatar
Pilar Cossio committed
705
#ifdef DEBUG
706 707
	    cout << "Reduction " << mpi_rank << " Map " << i << " Prob " << pProbMap.Total << " Const " << pProbMap.Constoadd  << "\n";     
#endif
708
	    tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
709

710 711 712 713 714
	  }
	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
715
	{
716 717 718 719 720 721
	  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;
722
              //temporary array that has the mpirank for the highest pProb.constant
723 724 725 726 727
	    }
	  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
728
		{
729
		  if (mpi_rank == 0) printf("Error: Could not find highest probability\n");
David Rohr's avatar
David Rohr committed
730
		}
731
	      else if (tmpi2[i] != 0) //Skip if rank 0 already has highest probability
David Rohr's avatar
David Rohr committed
732
		{
733 734 735 736 737 738 739 740
		  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
741
		}
742 743 744
	    }
	  delete[] tmpi1;
	  delete[] tmpi2;
745
	}
746

David Rohr's avatar
David Rohr committed
747
	if (mpi_rank == 0)
748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764
	  {
	    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)
765
	{
766 767 768 769 770 771
	  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++)
	    {
772
	      for (int j = 0;j < param.nTotGridAngles;j++)
773
                {
774 775 776 777
		  //	      tmp1[i] = pProb.getProbMap(i).Constoadd;
		  //	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
		  tmp1[i * param.nTotGridAngles + j]= pProb.getProbAngle(i, j).ConstAngle;
		}
778
	    }
779

780 781 782 783
	  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
784
		{
785 786
		  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
787
		}
788 789 790 791 792
	    }
	  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++)
793
		{
794 795 796 797 798 799
		  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];
		    }
800
		}
801 802 803 804 805 806 807
	    }
	  delete[] tmp1;
	  delete[] tmp2;
	  delete[] tmp3;
	}
    }
#endif
David Rohr's avatar
David Rohr committed
808

809 810

  // ************* Writing Out Probabilities ***************
811 812
  if (mpi_rank == 0)
    {
813
 
814 815
      // Output for Angle Probability File
      ofstream angProbfile;
816 817 818
      if(param.param_device.writeAngles)
	{
	  angProbfile.open ("ANG_PROB");
819
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
820 821
          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" ;};
822
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
823 824
	  //          angProbfile <<"Model Used: " << modelfile.c_str() << "\n";
	  //          angProbfile <<"Input Used: " << infile.c_str() << "\n";
825
	}
826
      // Output for Cross Correlation File
827 828 829 830
      ofstream ccProbfile;
      if(param.param_device.writeCC)
	{
	  ccProbfile.open ("CROSS_CORRELATION");
831
	  ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
832 833 834
          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";
835
          ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
836 837
	}

838 839 840 841 842 843 844 845 846
      // 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{
847
	  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";}
848 849 850
      }else { 
	if(param.usepsf){
	  //     if( localcc[rx * param.param_device.NumberPixels + ry] <
Pilar Cossio's avatar
Pilar Cossio committed
851
	  outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - q4 -PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";
852
	}else{
Pilar Cossio's avatar
Pilar Cossio committed
853
          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";
854
        }}
855
      if(param.writeCTF) outputProbFile << " RefMap:  MapNumber ; CTFMaxParm: defocus - b-Env (B ref. Penzeck 2010)\n";
Pilar Cossio's avatar
Pilar Cossio committed
856
      if(param.yespriorAngles) outputProbFile << "**** Remark: Using Prior Proability in Angles ****\n";
857 858
      outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";

Pilar Cossio's avatar
Pilar Cossio committed
859
       
860
      // Loop over reference maps
861
      // ************* Over all maps ***************
862

863 864 865 866 867
      for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
	  // **** Total Probability ***
	  bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

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

872 873
	    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: ";
874 875 876
            // *** 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)) << " ";

877

878
	  }else{ 
879 880 881 882 883 884 885 886 887
	    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
888 889 890 891 892 893

	  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] << " [] ";
894 895 896 897 898 899
	  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] " ;
900 901
	  if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_norm << " [] " ;}else{outputProbFile << "N.A." << " [] ";}
	  if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_mu << " [] ";}else{outputProbFile << "N.A." << " [] ";} 
902
	  outputProbFile << "\n";
903 904 905 906 907 908 909 910 911

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

915
	  //*************** Writing Individual Angle probabilities
916 917 918
	  if(param.param_device.writeAngles)
	    {
	      for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
David Rohr's avatar
David Rohr committed
919
		{
920
		  bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
921

922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944
		  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";
		      }
		  }
945 946
		}
	    }
947
	
948
	  //************* Writing Cross-Correlations if requiered
949
          //************* This is currently not in the manual *****
950
	  if(param.param_device.writeCC){
951

952 953 954 955 956
	    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
957
            int used[(param.param_device.NumberPixels+1) * (param.param_device.NumberPixels+1)];
958

959 960
	    halfPix = param.param_device.NumberPixels / 2 ;
	    // Ordering the centers of the Cross Correlation
961

962 963 964 965 966
	    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
967
			used[ rx * param.param_device.NumberPixels + ry ]= 0;
968 969 970 971 972 973 974 975 976 977 978 979 980 981