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

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

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

10 11
#include <mpi.h>

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

18 19 20 21 22 23 24 25 26
#include <fstream>
#include <boost/program_options.hpp>
#include <iostream>
#include <algorithm>
#include <iterator>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cmath>
27

28
#ifdef WITH_OPENMP
29
#include <omp.h>
30
#endif
31 32 33 34 35 36 37 38 39

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

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

42 43 44 45 46 47
#ifdef BIOEM_USE_NVTX
#include "nvToolsExt.h"

const uint32_t colors[] = { 0x0000ff00, 0x000000ff, 0x00ffff00, 0x00ff00ff, 0x0000ffff, 0x00ff0000, 0x00ffffff };
const int num_colors = sizeof(colors)/sizeof(colors[0]);

48 49 50 51 52 53 54 55 56 57 58 59
#define cuda_custom_timeslot(name,cid) {		\
    int color_id = cid;					\
    color_id = color_id%num_colors;			\
    nvtxEventAttributes_t eventAttrib = {0};		\
    eventAttrib.version = NVTX_VERSION;			\
    eventAttrib.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;	\
    eventAttrib.colorType = NVTX_COLOR_ARGB;		\
    eventAttrib.color = colors[color_id];		\
    eventAttrib.messageType = NVTX_MESSAGE_TYPE_ASCII;	\
    eventAttrib.message.ascii = name;			\
    nvtxRangePushEx(&eventAttrib);			\
  }
60 61 62 63 64
#define cuda_custom_timeslot_end nvtxRangePop();
#else
#define cuda_custom_timeslot(name,cid)
#define cuda_custom_timeslot_end
#endif
65

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

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

using namespace std;

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

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

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

bioem::~bioem()
{
}

int bioem::configure(int ac, char* av[])
{
99 100 101 102 103 104 105 106
  // **************************************************************************************
  // **** Configuration Routine using boost for extracting parameters, models and maps ****
  // **************************************************************************************
  // ****** And Precalculating necessary grids, map crosscorrelations and kernels  ********
  // *************************************************************************************

  HighResTimer timer;

107
  std::string infile, modelfile, mapfile,Inputanglefile,Inputbestmap; 
108 109 110 111 112 113 114 115 116 117 118
  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
119
      yesoutfilename=false;
120 121 122 123 124 125 126 127 128 129 130 131 132 133

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

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

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

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

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

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

260 261 262 263
      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
264

265 266 267 268
      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
269
	param.readParameters(infile.c_str());
David Rohr's avatar
David Rohr committed
270

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

274

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

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

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

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

286
      if (DebugOutput >= 2 && mpi_rank == 0) printf("Reading Input Data Time: %f\n", timer.GetCurrentElapsedTime());
287
    
288
      if(param.param_device.writeCC && mpi_size>1){
289 290
	cout << "Exiting::: WRITE CROSS-CORRELATION ONLY VAILD FOR 1 MPI PROCESS\n";
        exit(1);
291
      }
Pilar Cossio's avatar
Pilar Cossio committed
292

293 294
      // Generating Grids of orientations 
      if(!param.printModel)param.CalculateGridsParam(Inputanglefile.c_str());
295
    }
296

297
#ifdef WITH_MPI
298 299 300 301 302 303 304

  // Contros for MPI
  if(mpi_size > param.nTotGridAngles){
    cout << "EXIT: Wrong MPI setup More MPI processes than orientations\n"; exit(1);
  }

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

312 313
      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);
314

Pilar Cossio's avatar
Pilar Cossio committed
315
#ifdef PILAR_DEBUG
316 317
      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
318 319

#endif
320
      //****refCtf, CtfParam, angles automatically filled by precalculate function bellow
321 322 323 324 325 326 327 328 329

      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
330

331 332
    }
#endif
333

334 335 336
  // ****************** Precalculating Necessary Stuff *********************
  if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
  param.PrepareFFTs();
337

338 339 340 341 342 343 344
  if (DebugOutput >= 2 && mpi_rank == 0)
    {
      printf("Time Prepare FFTs %f\n", timer.GetCurrentElapsedTime());
      timer.ResetStart();
    }
  precalculate();

345
  // ****************** For debugging *********************
346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364
  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();
    }
365 366 367

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

368
  deviceInit();
369

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

  return(0);
373 374
}

375 376
void bioem::cleanup()
{
377
  //Deleting allocated pointers
378 379
  free_device_host(pProb.ptr);
  RefMap.freePointers();
380 381
}

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

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

411 412 413
  // **************************************************************************************
  // ********** Secondary routine for printing out the only best projection ***************
  // **************************************************************************************
David Rohr's avatar
David Rohr committed
414

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

417 418 419 420 421
    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;
422

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

427
    cout << "...... Calculating Projection .......................\n " ;
428

429
    createProjection(0, proj_mapsFFT);
430

431
    cout << "...... Calculating Convolution .......................\n " ;
David Rohr's avatar
David Rohr committed
432

433
    createConvolutedProjectionMap(0, 0, proj_mapsFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
434

435
  }
436

437 438 439
  // **************************************************************************************
  // **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
  // **************************************************************************************
440

441 442
  // **** 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
443

444
  if (mpi_rank == 0) printf("\tInitializing Probabilities\n");
445

446

447 448 449 450
  // 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
451

452
      pProbMap.Total = 0.0;
453 454
      pProbMap.Constoadd = -FLT_MAX; //Problem if using double presicion

455
      if (param.param_device.writeAngles)
456
	{
457 458 459 460 461
	  for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
	    {
	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);

	      pProbAngle.forAngles = 0.0;
462
	      pProbAngle.ConstAngle = -FLT_MAX;
463 464
	    }
	}
465

466
      if (param.param_device.writeCC)
467
	{      int  cc=0;
468
	  for (int cent_x = 0; cent_x < param.param_device.NumberPixels; cent_x = cent_x + param.param_device.CCdisplace)
469
	    {
470
	      for (int cent_y = 0; cent_y < param.param_device.NumberPixels; cent_y = cent_y + param.param_device.CCdisplace)
471
		{
472
		  bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
473 474 475 476 477 478 479 480 481
		  //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++;
482
		}
483
	    }
484
	  if(!FFTAlgo){cout << "Cross correlation calculation must be with enviormental variable FFTALGO=1\n"; exit(1);}
485 486
	}                 
    }
487 488 489

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

490
  // **************************************************************************************
491

492 493 494 495 496 497 498 499 500 501 502
  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
503
  //******** Alocating Vectors *************
504 505 506 507 508 509 510 511 512
  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);

513 514


515 516 517 518
  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;

519 520 521

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

522 523 524 525 526 527 528
  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);
529 530 531

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

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

539 540 541
      for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
	{
	  mycomplex_t* proj_mapFFT = &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize];
542

543
	  // ***************************************************************************************
544 545
	  // ***** **** Internal Loop over PSF/CTF convolutions **** *****

546 547 548 549 550 551 552 553
	  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);

554
	  
555
	      if (DebugOutput >= 2) timer.ResetStart();
556 557 558 559 560 561
     	      myfloat_t amp,pha,env;

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

562 563 564
	      // ******************Internal loop over Reference images CUDA or OpenMP******************
	      // *** Comparing each calculated convoluted map with all experimental maps ***

565
	      compareRefMaps(iOrient, iConv, amp, pha, env, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
566 567

	      if (DebugOutput >= 2)
568
		{
569 570 571 572 573 574 575 576 577
		  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);
578
		}
579 580 581 582 583 584
	    }
	  if (DebugOutput >= 1)
	    {
	      printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient, timer2.GetCurrentElapsedTime(), mpi_rank);
	      timer2.ResetStart();
	    }
585
	}
586 587 588 589 590
    }
  //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
591

592
  deviceFinishRun();
593 594


595 596

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

David Rohr's avatar
David Rohr committed
598
#ifdef WITH_MPI
599 600 601 602 603 604 605 606 607 608 609 610 611
  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);
612

613 614 615
	for (int i = 0;i < RefMap.ntotRefMap;i++)
	  {
	    bioem_Probability_map& pProbMap = pProb.getProbMap(i);
616 617 618
#ifdef PILAR_DEBUG
	    cout << "Reduction " << mpi_rank << " Map " << i << " Prob " << pProbMap.Total << " Const " << pProbMap.Constoadd  << "\n";     
#endif
619
	    tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
620

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

David Rohr's avatar
David Rohr committed
658
	if (mpi_rank == 0)
659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675
	  {
	    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)
676
	{
677 678 679 680 681 682
	  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++)
	    {
683
	      for (int j = 0;j < param.nTotGridAngles;j++)
684
                {
685 686 687 688
		  //	      tmp1[i] = pProb.getProbMap(i).Constoadd;
		  //	      bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j);
		  tmp1[i * param.nTotGridAngles + j]= pProb.getProbAngle(i, j).ConstAngle;
		}
689
	    }
690

691 692 693 694
	  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
695
		{
696 697
		  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
698
		}
699 700 701 702 703
	    }
	  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++)
704
		{
705 706 707 708 709 710
		  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];
		    }
711
		}
712 713 714 715 716 717 718
	    }
	  delete[] tmp1;
	  delete[] tmp2;
	  delete[] tmp3;
	}
    }
#endif
David Rohr's avatar
David Rohr committed
719

720 721

  // ************* Writing Out Probabilities ***************
722 723
  if (mpi_rank == 0)
    {
724
 
725 726
      // Output for Angle Probability File
      ofstream angProbfile;
727 728 729
      if(param.param_device.writeAngles)
	{
	  angProbfile.open ("ANG_PROB");
730
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
731 732
          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" ;};
733
	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
734 735
	  //          angProbfile <<"Model Used: " << modelfile.c_str() << "\n";
	  //          angProbfile <<"Input Used: " << infile.c_str() << "\n";
736
	}
737
      // Output for Cross Correlation File
738 739 740 741
      ofstream ccProbfile;
      if(param.param_device.writeCC)
	{
	  ccProbfile.open ("CROSS_CORRELATION");
742
	  ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
743 744 745
          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";
746
          ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
747 748
	}

749 750 751 752 753 754 755 756 757
      // 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{
758
	  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";}
759 760 761 762 763
      }else { 
	if(param.usepsf){
	  //     if( localcc[rx * param.param_device.NumberPixels + ry] <
	  outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";
	}else{
764
          outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - CTF amp - CTF defocus - CTF B-Env - center x - center y - normalization - offsett \n";
765
        }}
766
      if(param.writeCTF) outputProbFile << " RefMap:  MapNumber ; CTFMaxParm: defocus - b-Env (B ref. Penzeck 2010)\n";
767 768 769
      outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";

      // Loop over reference maps
770
      // ************* Over all maps ***************
771

772 773 774 775 776
      for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
	  // **** Total Probability ***
	  bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

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

781 782
	    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: ";
783 784 785
            // *** 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)) << " ";

786

787
	  }else{ 
788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809
	    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)) << " ";
	  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] << " [ ] ";
	  outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[0] << " [ ] ";
	  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] " ;
	  outputProbFile << pProbMap.max.max_prob_norm << " [ ] " ; 
	  outputProbFile << pProbMap.max.max_prob_mu << " [ ] ";
	  outputProbFile << "\n";
810 811 812 813 814 815 816 817 818 819 820

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

824
	  //*************** Writing Individual Angle probabilities
825 826 827
	  if(param.param_device.writeAngles)
	    {
	      for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
David Rohr's avatar
David Rohr committed
828
		{
829
		  bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
830

831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853
		  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";
		      }
		  }
854 855
		}
	    }
856
	
857
	  //************* Writing Cross-Correlations if requiered
858
	  if(param.param_device.writeCC){
859

860 861 862 863 864
	    int  cc=0;
	    int halfPix;
	    int rx=0;
	    int ry=0;
	    myfloat_t localcc[ (param.param_device.NumberPixels+1) * (param.param_device.NumberPixels+1) ];
865

866 867
	    halfPix = param.param_device.NumberPixels / 2 ;
	    // Ordering the centers of the Cross Correlation
868

869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 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
	    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;
		  }
	      }

	    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;
		    }
		    cc++;
		  }
		//              ccProbfile << "\n";
	      }
	    if(!param.ignoreCCoff){
	      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
913
		{
914
		  for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry = ry + param.param_device.CCdisplace)
915
		    {
916
		      if(localcc[ rx * param.param_device.NumberPixels + ry ] > -FLT_MAX){
917
			ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
918
		      }else{
919
			ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << -FLT_MAX << "\n" ;
920
		      }
921
		      //				 cout << " cc " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] <<"\n" ;
922
		    }
923 924 925 926
		  ccProbfile << "\n";
		}			
	    }else{
	      for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels ; rx++)
927
		{
928
		  for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry++)
929
		    {
930
		      if(localcc[ rx * param.param_device.NumberPixels + ry ] > -FLT_MAX){
931
                        ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
932
		      }else{
933
                        ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << -FLT_MAX << "\n" ;
934
		      }
935
		    }
936
		  ccProbfile << "\n";
937
		}
938

939
	    }
940
	  }
941
	}
942

943 944 945 946 947 948 949 950
      if(param.param_device.writeAngles)
	{
	  angProbfile.close();
	}

      if(param.param_device.writeCC)
	{
	  ccProbfile.close();
951
	}
952

953 954 955 956
      outputProbFile.close();
    }

  return(0);
957 958
}

959
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)
960
{
961

962 963 964 965 966 967 968
  //***************************************************************************************
  //***** BioEM routine for comparing reference maps to convoluted maps *****
  if (FFTAlgo)
    {
      //With FFT Algorithm
#pragma omp parallel for schedule(dynamic, 1)
      for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
969
	{
970
	  const int num = omp_get_thread_num();
971
	  calculateCCFFT(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
972
	}
973 974 975 976 977 978
    }
  else
    {
      //Without FFT Algorithm
#pragma omp parallel for schedule(dynamic, 1)
      for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
979
	{
980
	  compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, conv_map, pProb, param.param_device, RefMap);
981
	}
982 983
    }
  return(0);
984 985
}

986
inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv,  myfloat_t amp, myfloat_t pha, myfloat_t env, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC)
987
{
988 989
  //***************************************************************************************
  //***** Calculating cross correlation in FFTALGOrithm *****
Pilar Cossio's avatar
Pilar Cossio committed
990

991 992 993 994 995 996
  for(int i = 0; i < param.param_device.NumberPixels; i++)
    {
      for(int j = 0; j < param.param_device.NumberPixels; j++) lCC[i * param.param_device.NumberPixels + j] = 0.f; 
    }


997 998 999 1000 1001 1002
  const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize];
  for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
    {
      localCCT[i][0] = localConvFFT[i][0] * RefMapFFT[i][0] + localConvFFT[i][1] * RefMapFFT[i][1];
      localCCT[i][1] = localConvFFT[i][1] * RefMapFFT[i][0] - localConvFFT[i][0] * RefMapFFT[i][1];
    }
1003

1004
  myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
1005 1006
 
  // printf("HereCC %p %f %d %d %d %d \n", &lCC[139 * param.param_device.NumberPixels + 139],lCC[139 * param.param_device.NumberPixels + 139],mpi_rank,iConv,iOrient,iRefMap);
1007

1008
  doRefMapFFT(iRefMap, iOrient, iConv, amp, pha, env, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
Pilar Cossio's avatar
Pilar Cossio committed
1009

1010

Pilar Cossio's avatar
Pilar Cossio committed
1011
#ifdef PILAR_DEBUG
1012 1013
  if (param.param_device.writeCC)
    {      int  cc=0;