bioem.cpp 74.3 KB
Newer Older
Pilar Cossio's avatar
License    
Pilar Cossio committed
1
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2
   < BioEM software for Bayesian inference of Electron Microscopy images>
Luka Stanisic's avatar
Luka Stanisic committed
3
4
   Copyright (C) 2017 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
        Luka Stanisic, Volker Lindenstruth and Gerhard Hummer.
5
   Max Planck Institute of Biophysics, Frankfurt, Germany.
Luka Stanisic's avatar
Luka Stanisic committed
6
7
8
   Frankfurt Institute for Advanced Studies, Goethe University Frankfurt,
   Germany.
   Max Planck Computing and Data Facility, Garching, Germany.
9

Luka Stanisic's avatar
Luka Stanisic 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

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
Luka Stanisic's avatar
Luka Stanisic committed
14

Pilar Cossio's avatar
Pilar Cossio committed
15
#ifdef WITH_MPI
16
17
#include <mpi.h>

Luka Stanisic's avatar
Luka Stanisic committed
18
19
20
21
22
#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
23
#endif
24

Luka Stanisic's avatar
Luka Stanisic committed
25
26
27
#include "MersenneTwister.h"
#include <algorithm>
#include <cmath>
28
#include <fstream>
Luka Stanisic's avatar
Luka Stanisic committed
29
#include <getopt.h>
30
31
#include <iostream>
#include <iterator>
Luka Stanisic's avatar
Luka Stanisic committed
32
#include <queue>
33
34
#include <stdio.h>
#include <stdlib.h>
Luka Stanisic's avatar
Luka Stanisic committed
35
#include <string.h>
36
#include <string>
Luka Stanisic's avatar
Luka Stanisic committed
37
#include <vector>
38

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

Luka Stanisic's avatar
Luka Stanisic committed
43
44
#include "autotuner.h"
#include "timer.h"
45
46
47
48
49
#include <fftw3.h>
#include <math.h>

#include "bioem.h"
#include "map.h"
Luka Stanisic's avatar
Luka Stanisic committed
50
51
#include "model.h"
#include "param.h"
52

53
54
55
#ifdef BIOEM_USE_NVTX
#include "nvToolsExt.h"

Luka Stanisic's avatar
Luka Stanisic committed
56
57
58
59
60
61
62
63
64
65
66
const uint32_t colors[] = {0x0000ff00, 0x000000ff, 0x00ffff00, 0x00ff00ff,
                           0x0000ffff, 0x00ff0000, 0x00ffffff};
const int num_colors = sizeof(colors) / sizeof(colors[0]);
enum myColor
{
  COLOR_PROJECTION,
  COLOR_CONVOLUTION,
  COLOR_COMPARISON,
  COLOR_WORKLOAD,
  COLOR_INIT
};
67

68
69
70
// Projection number is stored in category attribute
// Convolution number is stored in payload attribute

Luka Stanisic's avatar
Luka Stanisic committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#define cuda_custom_timeslot(name, iMap, iConv, 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.category = iMap;                                               \
    eventAttrib.payloadType = NVTX_PAYLOAD_TYPE_UNSIGNED_INT64;                \
    eventAttrib.payload.llValue = iConv;                                       \
    eventAttrib.messageType = NVTX_MESSAGE_TYPE_ASCII;                         \
    eventAttrib.message.ascii = name;                                          \
    nvtxRangePushEx(&eventAttrib);                                             \
86
  }
87
88
#define cuda_custom_timeslot_end nvtxRangePop();
#else
Luka Stanisic's avatar
Luka Stanisic committed
89
#define cuda_custom_timeslot(name, iMap, iConv, cid)
90
91
#define cuda_custom_timeslot_end
#endif
92

93
94
95
96
97
98
#include "bioem_algorithm.h"

using namespace std;

bioem::bioem()
{
Luka Stanisic's avatar
Luka Stanisic committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
  BioEMAlgo = getenv("BIOEM_ALGO") == NULL ? 1 : atoi(getenv("BIOEM_ALGO"));

  DebugOutput = getenv("BIOEM_DEBUG_OUTPUT") == NULL ?
                    0 :
                    atoi(getenv("BIOEM_DEBUG_OUTPUT"));

  if (getenv("BIOEM_PROJ_CONV_AT_ONCE") != NULL)
  {
    nProjectionsAtOnce = atoi(getenv("BIOEM_PROJ_CONV_AT_ONCE"));
    if (BioEMAlgo == 1 && getenv("GPU") && atoi(getenv("GPU")) &&
        nProjectionsAtOnce > 1)
    {
      printf("Warning: using parallel convolutions with GPUs can create race "
             "condition and lead to inaccurate results. "
             "BIOEM_PROJ_CONV_AT_ONCE is going to be set 1.\n");
      nProjectionsAtOnce = 1;
    }
  }
  else if (BioEMAlgo == 1)
    nProjectionsAtOnce = 1;
  else
    nProjectionsAtOnce =
        getenv("OMP_NUM_THREADS") == NULL ? 1 : atoi(getenv("OMP_NUM_THREADS"));

  if (getenv("BIOEM_CUDA_THREAD_COUNT") != NULL)
    CudaThreadCount = atoi(getenv("BIOEM_CUDA_THREAD_COUNT"));
  else if (BioEMAlgo == 1)
    CudaThreadCount = CUDA_THREAD_COUNT_ALGO1;
  else
    CudaThreadCount = CUDA_THREAD_COUNT_ALGO2;

Luka Stanisic's avatar
Luka Stanisic committed
130
  Autotuning = false;
131
132
}

Luka Stanisic's avatar
Luka Stanisic committed
133
134
135
bioem::~bioem() {}

void bioem::printOptions(myoption_t *myoptions, int myoptions_length)
136
{
Luka Stanisic's avatar
Luka Stanisic committed
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
  printf("\nCommand line inputs:\n");

  // Find longest column width
  int maxlen = 0;
  for (int i = 0; i < myoptions_length; i++)
  {
    if (myoptions[i].hidden)
      continue;
    if (maxlen < strlen(myoptions[i].name))
      maxlen = strlen(myoptions[i].name);
  }

  for (int i = 0; i < myoptions_length; i++)
  {
    if (myoptions[i].hidden)
      continue;
    printf("  --%-*s", maxlen, myoptions[i].name);
    if (myoptions[i].arg == required_argument)
      printf(" arg");
    else
      printf("    ");
    printf(" %s\n", myoptions[i].desc);
  }
  printf("\n");
161
162
}

Luka Stanisic's avatar
Luka Stanisic committed
163
int bioem::readOptions(int ac, char *av[])
164
{
165
166
  HighResTimer timer;

Luka Stanisic's avatar
Luka Stanisic committed
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
  // *** Inizialzing default variables ***
  std::string infile, modelfile, mapfile, Inputanglefile, Inputbestmap;
  Model.readPDB = false;
  param.param_device.writeAngles = 0;
  param.dumpMap = false;
  param.loadMap = false;
  param.printModel = false;
  RefMap.readMRC = false;
  RefMap.readMultMRC = false;
  param.notuniformangles = false;
  OutfileName = "Output_Probabilities";

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

  // Write your options here
  myoption_t myoptions[] = {
      {"Modelfile", required_argument, "(Mandatory) Name of model file", false},
      {"Particlesfile", required_argument,
       "(Mandatory) Name of particle-image file", false},
      {"Inputfile", required_argument,
       "(Mandatory) Name of input parameter file", false},
      {"PrintBestCalMap", required_argument,
       "(Optional) Only print best calculated map. NO BioEM!", true},
      {"ReadOrientation", required_argument,
       "(Optional) Read file name containing orientations", false},
      {"ReadPDB", no_argument, "(Optional) If reading model file in PDB format",
       false},
      {"ReadMRC", no_argument,
       "(Optional) If reading particle file in MRC format", false},
      {"ReadMultipleMRC", no_argument, "(Optional) If reading Multiple MRCs",
       false},
      {"DumpMaps", no_argument,
       "(Optional) Dump maps after they were read from particle-image file",
       false},
      {"LoadMapDump", no_argument, "(Optional) Read Maps from dump option",
       false},
      {"OutputFile", required_argument,
       "(Optional) For changing the outputfile name", false},
      {"help", no_argument, "(Optional) Produce help message", false}};
  int myoptions_length = sizeof(myoptions) / sizeof(myoption_t);

  // If not all Mandatory parameters are defined
  if ((ac < 2))
  {
    printf("Error: Need to specify all mandatory options\n");
    printOptions(myoptions, myoptions_length);
    return 1;
  }
Pilar Cossio's avatar
Pilar Cossio committed
215

Luka Stanisic's avatar
Luka Stanisic committed
216
217
218
219
220
221
222
223
  // Creating options structure for getopt_long()
  struct option *long_options =
      (option *) calloc((myoptions_length + 1), sizeof(option));
  for (int i = 0; i < myoptions_length; i++)
  {
    long_options[i].name = myoptions[i].name;
    long_options[i].has_arg = myoptions[i].arg;
  }
Pilar Cossio's avatar
Pilar Cossio committed
224

Luka Stanisic's avatar
Luka Stanisic committed
225
226
227
228
229
230
231
232
233
234
  int myopt;
  while (1)
  {
    /* getopt_long stores the option index here. */
    int option_index = 0;
    myopt = getopt_long(ac, av, "", long_options, &option_index);

    /* Detect the end of the options. */
    if (myopt == -1)
      break;
Pilar Cossio's avatar
Pilar Cossio committed
235

Luka Stanisic's avatar
Luka Stanisic committed
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
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
    switch (myopt)
    {
      case 0:
#ifdef DEBUG
        printf("option %s", long_options[option_index].name);
        if (optarg)
          printf(" with arg %s", optarg);
        printf("\n");
#endif
        // Here write actions for each option
        if (!strcmp(long_options[option_index].name, "help"))
        {
          cout << "Usage: options_description [options]\n";
          printOptions(myoptions, myoptions_length);
          return 1;
        }
        if (!strcmp(long_options[option_index].name, "Inputfile"))
        {
          cout << "Input file is: " << optarg << "\n";
          infile = optarg;
        }
        if (!strcmp(long_options[option_index].name, "Modelfile"))
        {
          cout << "Model file is: " << optarg << "\n";
          modelfile = optarg;
        }
        if (!strcmp(long_options[option_index].name, "ReadPDB"))
        {
          cout << "Reading model file in PDB format.\n";
          Model.readPDB = true;
        }
        if (!strcmp(long_options[option_index].name, "ReadOrientation"))
        {
          cout << "Reading Orientation from file: " << optarg << "\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 = optarg;
          param.notuniformangles = true;
        }
        if (!strcmp(long_options[option_index].name, "OutputFile"))
        {
          cout << "Writing OUTPUT to: " << optarg << "\n";
          OutfileName = optarg;
        }
        if (!strcmp(long_options[option_index].name, "PrintBestCalMap"))
        {
          cout << "Reading Best Parameters from file: " << optarg << "\n";
          Inputbestmap = optarg;
          param.printModel = true;
        }
        if (!strcmp(long_options[option_index].name, "ReadMRC"))
        {
          cout << "Reading particle file in MRC format.\n";
          RefMap.readMRC = true;
        }
        if (!strcmp(long_options[option_index].name, "ReadMultipleMRC"))
        {
          cout << "Reading Multiple MRCs.\n";
          RefMap.readMultMRC = true;
        }
        if (!strcmp(long_options[option_index].name, "DumpMaps"))
        {
          cout << "Dumping Maps after reading from file.\n";
          param.dumpMap = true;
        }
        if (!strcmp(long_options[option_index].name, "LoadMapDump"))
        {
          cout << "Loading Map dump.\n";
          param.loadMap = true;
        }
        if (!strcmp(long_options[option_index].name, "Particlesfile"))
        {
          cout << "Particle file is: " << optarg << "\n";
          mapfile = optarg;
        }
        break;
      case '?':
        /* getopt_long already printed an error message. */
        printOptions(myoptions, myoptions_length);
        return 1;
      default:
        abort();
325
    }
Luka Stanisic's avatar
Luka Stanisic committed
326
327
328
329
330
331
332
333
334
335
336
  }
  /* Print any remaining command line arguments (not options) and exit */
  if (optind < ac)
  {
    printf("Error: non-option ARGV-elements: ");
    while (optind < ac)
      printf("%s ", av[optind++]);
    putchar('\n');
    printOptions(myoptions, myoptions_length);
    return 1;
  }
337

Luka Stanisic's avatar
Luka Stanisic committed
338
339
340
341
342
343
  // check for consitency in multiple MRCs
  if (RefMap.readMultMRC && not(RefMap.readMRC))
  {
    cout << "For Multiple MRCs command --ReadMRC is necesary too";
    exit(1);
  }
344

Luka Stanisic's avatar
Luka Stanisic committed
345
346
347
348
349
  if (!Model.readPDB)
  {
    cout << "Note: Reading model in simple text format (not PDB)\n";
    cout << "----  x   y   z  radius  density ------- \n";
  }
Pilar Cossio's avatar
Pilar Cossio committed
350

Luka Stanisic's avatar
Luka Stanisic committed
351
352
  if (DebugOutput >= 2 && mpi_rank == 0)
    timer.ResetStart();
353

Luka Stanisic's avatar
Luka Stanisic committed
354
355
356
357
358
359
360
361
362
363
364
365
366
  // *** Reading Parameter Input ***
  if (!param.printModel)
  {
    // Standard definition for BioEM
    param.readParameters(infile.c_str());
    // *** Reading Particle Maps Input ***
    RefMap.readRefMaps(param, mapfile.c_str());
  }
  else
  {
    // Reading parameters for only writting down Best projection
    param.forprintBest(Inputbestmap.c_str());
  }
Pilar Cossio's avatar
Pilar Cossio committed
367

Luka Stanisic's avatar
Luka Stanisic committed
368
369
  // *** Reading Model Input ***
  Model.readModel(param, modelfile.c_str());
370

Luka Stanisic's avatar
Luka Stanisic committed
371
372
  cout << "**NOTE:: look at file COORDREAD to confirm that the Model "
          "coordinates are correct\n";
Pilar Cossio's avatar
Pilar Cossio committed
373

Luka Stanisic's avatar
Luka Stanisic committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
  if (DebugOutput >= 2 && mpi_rank == 0)
    printf("Reading Input Data Time: %f\n", timer.GetCurrentElapsedTime());

  // Generating Grids of orientations
  if (!param.printModel)
    param.CalculateGridsParam(Inputanglefile.c_str());

  return (0);
}

int bioem::configure(int ac, char *av[])
{
  // **************************************************************************************
  // **** Configuration Routine using getopts for extracting parameters, models
  // and maps ****
  // **************************************************************************************
  // ****** And Precalculating necessary grids, map crosscorrelations and
  // kernels  ********
  // *************************************************************************************

  HighResTimer timer;

  if (mpi_rank == 0 && readOptions(ac, av))
    return 1;

#ifdef WITH_MPI
400

Luka Stanisic's avatar
Luka Stanisic committed
401
402
403
404
405
406
407
408
409
410
  // ********************* MPI inizialization/ Transfer of
  // parameters******************
  if (mpi_size > 1)
  {
    if (DebugOutput >= 2 && mpi_rank == 0)
      timer.ResetStart();
    MPI_Bcast(&param, sizeof(param), MPI_BYTE, 0, MPI_COMM_WORLD);
    // We have to reinitialize all pointers !!!!!!!!!!!!
    if (mpi_rank != 0)
      param.angprior = NULL;
411

Luka Stanisic's avatar
Luka Stanisic committed
412
413
414
415
416
    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);
Pilar Cossio's avatar
Pilar Cossio committed
417

Luka Stanisic's avatar
Luka Stanisic committed
418
419
420
421
422
423
424
#ifdef DEBUG
    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";
425
    }
Luka Stanisic's avatar
Luka Stanisic committed
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448

#endif
    //****refCtf, CtfParam, angles automatically filled by precalculate function
    // bellow

    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());
  }
449
#endif
450

451
  // ****************** Precalculating Necessary Stuff *********************
Luka Stanisic's avatar
Luka Stanisic committed
452
453
  if (DebugOutput >= 2 && mpi_rank == 0)
    timer.ResetStart();
454
  param.PrepareFFTs();
455

456
  if (DebugOutput >= 2 && mpi_rank == 0)
Luka Stanisic's avatar
Luka Stanisic committed
457
458
459
460
  {
    printf("Time Prepare FFTs %f\n", timer.GetCurrentElapsedTime());
    timer.ResetStart();
  }
461
462
  precalculate();

463
  // ****************** For debugging *********************
464
  if (getenv("BIOEM_DEBUG_BREAK"))
Luka Stanisic's avatar
Luka Stanisic committed
465
466
467
468
469
470
471
  {
    const int cut = atoi(getenv("BIOEM_DEBUG_BREAK"));
    if (param.nTotGridAngles > cut)
      param.nTotGridAngles = cut;
    if (param.nTotCTFs > cut)
      param.nTotCTFs = cut;
  }
472
473

  if (DebugOutput >= 2 && mpi_rank == 0)
Luka Stanisic's avatar
Luka Stanisic committed
474
475
476
477
478
479
480
  {
    printf("Time Precalculate %f\n", timer.GetCurrentElapsedTime());
    timer.ResetStart();
  }

  // Number of parallel Convolutions and Comparisons
  param.nTotParallelConv = min(param.nTotCTFs, nProjectionsAtOnce);
Luka Stanisic's avatar
Luka Stanisic committed
481
482

  // ****************** For autotuning **********************
Luka Stanisic's avatar
Luka Stanisic committed
483
484
485
486
487
488
489
490
491
492
493
494
  if ((getenv("GPU") && atoi(getenv("GPU"))) && (BioEMAlgo == 1) &&
      ((!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);
  }
Luka Stanisic's avatar
Luka Stanisic committed
495
  else
Luka Stanisic's avatar
Luka Stanisic committed
496
497
498
  {
    Autotuning = false;
    if (mpi_rank == 0)
Luka Stanisic's avatar
Luka Stanisic committed
499
    {
Luka Stanisic's avatar
Luka Stanisic committed
500
501
502
503
504
505
506
507
      printf("Autotuning of GPUWorkload disabled");
      if (getenv("GPU") && atoi(getenv("GPU")))
        printf(", using GPUWorkload: %d%%\n",
               (getenv("GPUWORKLOAD") && (atoi(getenv("GPUWORKLOAD")) != -1)) ?
                   atoi(getenv("GPUWORKLOAD")) :
                   100);
      else
        printf(", please enable GPUs\n");
Luka Stanisic's avatar
Luka Stanisic committed
508
    }
Luka Stanisic's avatar
Luka Stanisic committed
509
  }
Luka Stanisic's avatar
Luka Stanisic committed
510

511
512
513
  // ****************** Initializing pointers *********************

  deviceInit();
514
515

  if (DebugOutput >= 2 && mpi_rank == 0)
Luka Stanisic's avatar
Luka Stanisic committed
516
517
518
519
  {
    printf("Time Device Init %f\n", timer.GetCurrentElapsedTime());
    timer.ResetStart();
  }
520

Luka Stanisic's avatar
Luka Stanisic committed
521
522
  if (!param.printModel)
    pProb.init(RefMap.ntotRefMap, param.nTotGridAngles, *this);
523

524
  if (DebugOutput >= 2 && mpi_rank == 0)
Luka Stanisic's avatar
Luka Stanisic committed
525
526
527
528
  {
    printf("Time Init Probabilities %f\n", timer.GetCurrentElapsedTime());
    timer.ResetStart();
  }
529

Luka Stanisic's avatar
Luka Stanisic committed
530
  return (0);
531
532
}

533
534
void bioem::cleanup()
{
Luka Stanisic's avatar
Luka Stanisic committed
535
  // Deleting allocated pointers
536
537
  free_device_host(pProb.ptr);
  RefMap.freePointers();
538
539
}

540
541
int bioem::precalculate()
{
542
  // **************************************************************************************
Luka Stanisic's avatar
Luka Stanisic committed
543
544
  // **Precalculating Routine of Orientation grids, Map crosscorrelations and
  // CTF Kernels**
545
546
  // **************************************************************************************
  HighResTimer timer;
Luka Stanisic's avatar
Luka Stanisic committed
547
548
549
550
551
552
  if (DebugOutput >= 2)
  {
    printf("\tTime Precalculate Grids Param: %f\n",
           timer.GetCurrentElapsedTime());
    timer.ResetStart();
  }
553
554
555
  // Precalculating CTF Kernels stored in class Param
  param.CalculateRefCTF();

Luka Stanisic's avatar
Luka Stanisic committed
556
557
558
559
560
561
562
563
564
565
  if (DebugOutput >= 2)
  {
    printf("\tTime Precalculate CTFs: %f\n", timer.GetCurrentElapsedTime());
    timer.ResetStart();
  }
  // Precalculate Maps
  if (!param.printModel)
    RefMap.precalculate(param, *this);
  if (DebugOutput >= 2)
    printf("\tTime Precalculate Maps: %f\n", timer.GetCurrentElapsedTime());
566

Luka Stanisic's avatar
Luka Stanisic committed
567
  return (0);
568
569
}

Luka Stanisic's avatar
Luka Stanisic committed
570
int bioem::printModel()
571
{
572
  // **************************************************************************************
Luka Stanisic's avatar
Luka Stanisic committed
573
574
  // ********** Secondary routine for printing out the only best projection
  // ***************
575
  // **************************************************************************************
David Rohr's avatar
David Rohr committed
576

Luka Stanisic's avatar
Luka Stanisic committed
577
578
579
580
581
  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;
582

Luka Stanisic's avatar
Luka Stanisic committed
583
584
585
586
587
588
589
590
591
  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);
592

Luka Stanisic's avatar
Luka Stanisic committed
593
  cout << "...... Calculating Projection .......................\n ";
594

Luka Stanisic's avatar
Luka Stanisic committed
595
  createProjection(0, proj_mapsFFT);
596

Luka Stanisic's avatar
Luka Stanisic committed
597
  cout << "...... Calculating Convolution .......................\n ";
David Rohr's avatar
David Rohr committed
598

Luka Stanisic's avatar
Luka Stanisic committed
599
600
  createConvolutedProjectionMap_noFFT(proj_mapsFFT, conv_map, conv_mapFFT,
                                      sumCONV, sumsquareCONV);
601

Luka Stanisic's avatar
Luka Stanisic committed
602
603
  return (0);
}
604

Luka Stanisic's avatar
Luka Stanisic committed
605
606
int bioem::run()
{
607
  // **************************************************************************************
Luka Stanisic's avatar
Luka Stanisic committed
608
609
  // **** Main BioEM routine, projects, convolutes and compares with Map using
  // OpenMP ****
610
  // **************************************************************************************
611

Luka Stanisic's avatar
Luka Stanisic committed
612
613
614
615
  // **** If we want to control the number of threads ->
  // omp_set_num_threads(XX); ******
  // ****************** Declarying class of Probability Pointer
  // *************************
616
  cuda_custom_timeslot("Initialization", -1, -1, COLOR_INIT);
Luka Stanisic's avatar
Luka Stanisic committed
617
618
  if (mpi_rank == 0)
    printf("\tInitializing Probabilities\n");
619

Pilar Cossio's avatar
Pilar Cossio committed
620
  // Contros for MPI
Luka Stanisic's avatar
Luka Stanisic committed
621
622
623
624
  if (mpi_size > param.nTotGridAngles)
  {
    cout << "EXIT: Wrong MPI setup More MPI processes than orientations\n";
    exit(1);
Pilar Cossio's avatar
Pilar Cossio committed
625
  }
626

627
  // Inizialzing Probabilites to zero and constant to -Infinity
Luka Stanisic's avatar
Luka Stanisic committed
628
629
630
  for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap++)
  {
    bioem_Probability_map &pProbMap = pProb.getProbMap(iRefMap);
David Rohr's avatar
David Rohr committed
631

Luka Stanisic's avatar
Luka Stanisic committed
632
633
    pProbMap.Total = 0.0;
    pProbMap.Constoadd = MIN_PROB;
634

Luka Stanisic's avatar
Luka Stanisic committed
635
636
637
638
639
640
    if (param.param_device.writeAngles)
    {
      for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
      {
        bioem_Probability_angle &pProbAngle =
            pProb.getProbAngle(iRefMap, iOrient);
641

Luka Stanisic's avatar
Luka Stanisic committed
642
643
644
645
646
        pProbAngle.forAngles = 0.0;
        pProbAngle.ConstAngle = MIN_PROB;
      }
    }
  }
647

648
  // **************************************************************************************
649

650
651
  deviceStartRun();

Luka Stanisic's avatar
Luka Stanisic committed
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
  // ******************************** MAIN CYCLE
  // ******************************************

  mycomplex_t *proj_mapsFFT;
  mycomplex_t *conv_mapsFFT;
  myparam5_t *comp_params =
      new myparam5_t[param.nTotParallelConv * PIPELINE_LVL];
  int iPipeline = 0;

  // 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
  //******** Allocating Vectors *************
  proj_mapsFFT = (mycomplex_t *) myfftw_malloc(
      sizeof(mycomplex_t) * ProjMapSize * nProjectionsAtOnce);
  conv_mapsFFT =
      (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.FFTMapSize *
                                    param.nTotParallelConv * PIPELINE_LVL);

  cuda_custom_timeslot_end; // Ending initialization
675
676
677

  HighResTimer timer, timer2;

Luka Stanisic's avatar
Luka Stanisic committed
678
679
680
  /* Autotuning */
  Autotuner aut;
  if (Autotuning)
Luka Stanisic's avatar
Luka Stanisic committed
681
682
683
684
  {
    aut.Initialize(AUTOTUNING_ALGORITHM, FIRST_STABLE);
    rebalanceWrapper(aut.Workload());
  }
685

Luka Stanisic's avatar
Luka Stanisic committed
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
  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);

  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;
701

702
703
  /* Vectors for computing statistic on different parts of the code */
  TimeStat ts((iOrientEnd - iOrientStart), param.nTotCTFs);
Luka Stanisic's avatar
Luka Stanisic committed
704
705
  if (DebugOutput >= 1)
    ts.InitTimeStat(4);
706

Luka Stanisic's avatar
Luka Stanisic committed
707
708
  // **************************Loop Over
  // orientations***************************************
709

Luka Stanisic's avatar
Luka Stanisic committed
710
711
712
713
714
715
716
  for (int iOrientAtOnce = iOrientStart; iOrientAtOnce < iOrientEnd;
       iOrientAtOnce += nProjectionsAtOnce)
  {
    // ***************************************************************************************
    // ***** Creating Projection for given orientation and transforming to
    // Fourier space *****
    if (DebugOutput >= 1)
717
    {
Luka Stanisic's avatar
Luka Stanisic committed
718
719
720
721
722
      timer2.ResetStart();
      timer.ResetStart();
    }
    int iOrientEndAtOnce =
        std::min(iOrientEnd, iOrientAtOnce + nProjectionsAtOnce);
723

Luka Stanisic's avatar
Luka Stanisic committed
724
725
// **************************Parallel orientations for projections at
// once***************
726

727
#pragma omp parallel for
Luka Stanisic's avatar
Luka Stanisic committed
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
    for (int iOrient = iOrientAtOnce; iOrient < iOrientEndAtOnce; iOrient++)
    {
      createProjection(iOrient,
                       &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize]);
    }
    if (DebugOutput >= 1)
    {
      ts.time = timer.GetCurrentElapsedTime();
      ts.Add(TS_PROJECTION);
      if (DebugOutput >= 2)
        printf("\tTime Projection %d-%d: %f (rank %d)\n", iOrientAtOnce,
               iOrientEndAtOnce - 1, ts.time, mpi_rank);
    }
    /* Recalibrate if needed */
    if (Autotuning && ((iOrientAtOnce - iOrientStart) % RECALIB_FACTOR == 0) &&
        ((iOrientEnd - iOrientAtOnce) > RECALIB_FACTOR) &&
        (iOrientAtOnce != iOrientStart))
    {
      aut.Reset();
      rebalanceWrapper(aut.Workload());
    }

    for (int iOrient = iOrientAtOnce; iOrient < iOrientEndAtOnce; iOrient++)
    {
      mycomplex_t *proj_mapFFT =
          &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize];

      // ***************************************************************************************
      // ***** **** Internal Loop over PSF/CTF convolutions **** *****
      for (int iConvAtOnce = 0; iConvAtOnce < param.nTotCTFs;
           iConvAtOnce += param.nTotParallelConv)
      {
        if (DebugOutput >= 1)
          timer.ResetStart();
        int iConvEndAtOnce =
            std::min(param.nTotCTFs, iConvAtOnce + param.nTotParallelConv);
        // Total number of convolutions that can be treated in this iteration in
        // parallel
        int maxParallelConv = iConvEndAtOnce - iConvAtOnce;
#pragma omp parallel for
        for (int iConv = iConvAtOnce; iConv < iConvEndAtOnce; iConv++)
        {
          // *** Calculating convolutions of projection map and
          // crosscorrelations ***
          int i =
              (iPipeline & 1) * param.nTotParallelConv + (iConv - iConvAtOnce);
          mycomplex_t *localmultFFT = &conv_mapsFFT[i * param.FFTMapSize];

          createConvolutedProjectionMap(iOrient, iConv, proj_mapFFT,
                                        localmultFFT, comp_params[i].sumC,
                                        comp_params[i].sumsquareC);

          comp_params[i].amp = param.CtfParam[iConv].pos[0];
          comp_params[i].pha = param.CtfParam[iConv].pos[1];
          comp_params[i].env = param.CtfParam[iConv].pos[2];
        }
        if (DebugOutput >= 1)
        {
          ts.time = timer.GetCurrentElapsedTime();
          ts.Add(TS_CONVOLUTION);
          if (DebugOutput >= 2)
            printf("\t\tTime Convolution %d %d-%d: %f (rank %d)\n", iOrient,
                   iConvAtOnce, iConvEndAtOnce - 1, ts.time, mpi_rank);
        }

        // ******************Internal loop over Reference images CUDA or
        // OpenMP******************
        // *** Comparing each calculated convoluted map with all experimental
        // maps ***
        ts.time = 0.;
        if ((DebugOutput >= 1) || (Autotuning && aut.Needed(iConvAtOnce)))
          timer.ResetStart();
        compareRefMaps(iPipeline++, iOrient, iConvAtOnce, maxParallelConv,
                       conv_mapsFFT, comp_params);
        if (DebugOutput >= 1)
        {
          ts.time = timer.GetCurrentElapsedTime();
          ts.Add(TS_COMPARISON);
        }
        if (DebugOutput >= 2)
        {
          if (Autotuning)
            printf("\t\tTime Comparison %d %d-%d: %f sec with GPU workload "
                   "%d%% (rank %d)\n",
                   iOrient, iConvAtOnce, iConvEndAtOnce - 1, ts.time,
                   aut.Workload(), mpi_rank);
          else
            printf("\t\tTime Comparison %d %d-%d: %f sec (rank %d)\n", iOrient,
                   iConvAtOnce, iConvEndAtOnce - 1, ts.time, mpi_rank);
        }
        if (Autotuning && aut.Needed(iConvAtOnce))
        {
          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);
          rebalanceWrapper(aut.Workload());
        }
      }
829
      if (DebugOutput >= 1)
Luka Stanisic's avatar
Luka Stanisic committed
830
831
832
833
834
835
836
      {
        ts.time = timer2.GetCurrentElapsedTime();
        ts.Add(TS_TPROJECTION);
        printf("\tTotal time for projection %d: %f (rank %d)\n", iOrient,
               ts.time, mpi_rank);
        timer2.ResetStart();
      }
837
    }
Luka Stanisic's avatar
Luka Stanisic committed
838
  }
839
840
  /* Statistical summary on different parts of the code */
  if (DebugOutput >= 1)
Luka Stanisic's avatar
Luka Stanisic committed
841
842
843
844
  {
    ts.PrintTimeStat(mpi_rank);
    ts.EmptyTimeStat();
  }
845

Luka Stanisic's avatar
Luka Stanisic committed
846
  // deallocating fftw_complex vector
847
  myfftw_free(proj_mapsFFT);
Luka Stanisic's avatar
Luka Stanisic committed
848
  myfftw_free(conv_mapsFFT);
David Rohr's avatar
David Rohr committed
849

850
  deviceFinishRun();
851

Luka Stanisic's avatar
Luka Stanisic committed
852
853
854
// *******************************************************************************
// ************* Collecing all the probabilities from MPI replicas
// ***************
David Rohr's avatar
David Rohr committed
855

David Rohr's avatar
David Rohr committed
856
#ifdef WITH_MPI
857
  if (mpi_size > 1)
Luka Stanisic's avatar
Luka Stanisic committed
858
859
860
861
  {
    if (DebugOutput >= 1 && mpi_rank == 0)
      timer.ResetStart();
    // Reduce Constant and summarize probabilities
862
    {
Luka Stanisic's avatar
Luka Stanisic committed
863
864
865
866
867
868
869
870
871
872
873
      myprob_t *tmp1 = new myprob_t[RefMap.ntotRefMap];
      myprob_t *tmp2 = new myprob_t[RefMap.ntotRefMap];
      myprob_t *tmp3 = new myprob_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);

      for (int i = 0; i < RefMap.ntotRefMap; i++)
874
      {
Luka Stanisic's avatar
Luka Stanisic committed
875
        bioem_Probability_map &pProbMap = pProb.getProbMap(i);
Pilar Cossio's avatar
Pilar Cossio committed
876
#ifdef DEBUG
Luka Stanisic's avatar
Luka Stanisic committed
877
878
        cout << "Reduction " << mpi_rank << " Map " << i << " Prob "
             << pProbMap.Total << " Const " << pProbMap.Constoadd << "\n";
879
#endif
Luka Stanisic's avatar
Luka Stanisic committed
880
        tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
881
      }
Luka Stanisic's avatar
Luka Stanisic committed
882
883
      MPI_Reduce(tmp1, tmp3, RefMap.ntotRefMap, MY_MPI_FLOAT, MPI_SUM, 0,
                 MPI_COMM_WORLD);
884

Luka Stanisic's avatar
Luka Stanisic committed
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
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
      // Find MaxProb
      MPI_Status mpistatus;
      {
        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;
          // temporary array that has the mpirank for the highest pProb.constant
        }
        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)
          {
            if (mpi_rank == 0)
              printf("Error: Could not find highest probability\n");
          }
          else if (tmpi2[i] !=
                   0) // Skip if rank 0 already has highest probability
          {
            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);
            }
          }
        }
        delete[] tmpi1;
        delete[] tmpi2;
      }

      if (mpi_rank == 0)
      {
        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());
941
    }
David Rohr's avatar
David Rohr committed
942

Luka Stanisic's avatar
Luka Stanisic committed
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
    // Angle Reduction and Probability summation for individual angles
    if (param.param_device.writeAngles)
    {
      const int count = RefMap.ntotRefMap * param.nTotGridAngles;
      myprob_t *tmp1 = new myprob_t[count];
      myprob_t *tmp2 = new myprob_t[count];
      myprob_t *tmp3 = new myprob_t[count];
      for (int i = 0; i < RefMap.ntotRefMap; i++)
      {
        for (int j = 0; j < param.nTotGridAngles; j++)
        {
          //	      tmp1[i] = pProb.getProbMap(i).Constoadd;
          //	      bioem_Probability_angle& pProbAngle =
          // pProb.getProbAngle(i, j);
          tmp1[i * param.nTotGridAngles + j] =
              pProb.getProbAngle(i, j).ConstAngle;
        }
      }

      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++)
        {
          bioem_Probability_angle &pProbAngle = pProb.getProbAngle(i, j);
          tmp1[i * param.nTotGridAngles + j] =
              pProbAngle.forAngles *
              exp(pProbAngle.ConstAngle - tmp2[i * param.nTotGridAngles + j]);
        }
      }
      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++)
        {
          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];
          }
        }
      }
      delete[] tmp1;
      delete[] tmp2;
      delete[] tmp3;
    }
  }
#endif
992
993

  // ************* Writing Out Probabilities ***************
994
  if (mpi_rank == 0)
Luka Stanisic's avatar
Luka Stanisic committed
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
  {

    // Output for Angle Probability File
    ofstream angProbfile;
    angProbfile.precision(OUTPUT_PRECISION);
    angProbfile.setf(ios::fixed);
    if (param.param_device.writeAngles)
    {
      angProbfile.open("ANG_PROB");
      angProbfile << "************************* HEADER:: NOTATION "
                     "*******************************************\n";
      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";
      };
      angProbfile << "************************* HEADER:: NOTATION "
                     "*******************************************\n";
      //          angProbfile <<"Model Used: " << modelfile.c_str() << "\n";
      //          angProbfile <<"Input Used: " << infile.c_str() << "\n";
    }

    // Output for Standard Probability
    ofstream outputProbFile;
    outputProbFile.precision(OUTPUT_PRECISION);
    outputProbFile.setf(ios::fixed);
    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
      {
        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";
      }
    }
    else
1052
    {
Luka Stanisic's avatar
Luka Stanisic committed
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
      if (param.usepsf)
      {
        //     if( localcc[rx * param.param_device.NumberPixels + ry] <
        outputProbFile << "Notation= RefMap:  MapNumber ; Maximizing Param: "
                          "MaxLogProb - q1 - q2 - q3 - q4 -PSF amp - PSF phase "
                          "- PSF envelope - center x - center y - "
                          "normalization - offsett \n";
      }
      else
      {
        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";
      }
1068
    }
Luka Stanisic's avatar
Luka Stanisic committed
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
    if (param.writeCTF)
      outputProbFile << " RefMap:  MapNumber ; CTFMaxParm: defocus - b-Env (B "
                        "ref. Penzeck 2010)\n";
    if (param.yespriorAngles)
      outputProbFile << "**** Remark: Using Prior Proability in Angles ****\n";
    outputProbFile << "************************* HEADER:: NOTATION "
                      "*******************************************\n\n";

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

    for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap++)
    {
      // **** Total Probability ***
      bioem_Probability_map &pProbMap = pProb.getProbMap(iRefMap);
1084

Luka Stanisic's avatar
Luka Stanisic committed
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
      // Controll for Value of Total Probability
      // cout << pProbMap.Total << " " <<  pProbMap.Constoadd << " " << FLT_MAX
      // <<" " << log(FLT_MAX) << "\n";
      if (pProbMap.Total > 1.e-38)
      {

        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: ";
        // *** 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))
                       << " ";
      }
      else
      {
        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]
                     << " [] ";
      if (param.doquater)
        outputProbFile << param.angles[pProbMap.max.max_prob_orient].quat4
                       << " [] ";
      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";

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

      //*************** Writing Individual Angle probabilities
      if (param.param_device.writeAngles)
      {
        // Finding the best param.param_device.writeAngles probabilities
        // This implementation is clean, but not the most optimal one
        // and it supposes param.param_device.writeAngles <<
        // param.nTotGridAngles
        unsigned K =
            param.param_device.writeAngles; // number of best probabilities
                                            // clang-format off
        std::priority_queue<std::pair<double, int>,
                            std::vector<std::pair<double, int> >,
                            std::greater<std::pair<double, int> > >
            q;
        // clang-format on
        for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
        {
          bioem_Probability_angle &pProbAngle =
              pProb.getProbAngle(iRefMap, iOrient);

          myprob_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 (q.size() < K)
            q.push(std::pair<double, int>(logp, iOrient));
          else if (q.top().first < logp)
          {
            q.pop();
            q.push(std::pair<double, int>(logp, iOrient));
          }
        }
        K = q.size();
        int *rev_iOrient = (int *) malloc(K * sizeof(int));
        myprob_t *rev_logp = (myprob_t *) malloc(K * sizeof(myprob_t));
        for (int i = K - 1; i >= 0; i--)
        {
          rev_iOrient[i] = q.top().second;
          rev_logp[i] = q.top().first;
          q.pop();
        }
        for (unsigned i = 0; i < K; i++)
        {
          int iOrient = rev_iOrient[i];
          bioem_Probability_angle &pProbAngle =
              pProb.getProbAngle(iRefMap, iOrient);
          myprob_t logp = rev_logp[i];

          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";
            }
          }
        }
        free(rev_iOrient);
        free(rev_logp);
      }
    }

    if (param.param_device.writeAngles)
    {
      angProbfile.close();
    }

    outputProbFile.close();
  }

  return (0);
1322
1323
}

Luka Stanisic's avatar
Luka Stanisic committed
1324
1325
1326
int bioem::compareRefMaps(int iPipeline, int iOrient, int iConvStart,
                          int maxParallelConv, mycomplex_t *localmultFFT,
                          myparam5_t *comp_params, const int startMap)
1327
{
1328
1329
  //***************************************************************************************
  //***** BioEM routine for comparing reference maps to convoluted maps *****
1330
  //***************************************************************************************
Luka Stanisic's avatar
Luka Stanisic committed
1331
1332
1333
1334
1335
1336
  cuda_custom_timeslot("Comparison", iOrient, iConvStart, COLOR_COMPARISON);

  int k = (iPipeline & 1) * param.nTotParallelConv;

  if (BioEMAlgo == 1)
  {
1337
#pragma omp parallel for schedule(dynamic, 1)
Luka Stanisic's avatar
Luka Stanisic committed
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
    for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap++)
    {
      const int num = omp_get_thread_num();
      for (int iConv = 0; iConv < maxParallelConv; iConv++)
      {
        calculateCCFFT(iRefMap, &localmultFFT[(k + iConv) * param.FFTMapSize],
                       param.fft_scratch_complex[num],
                       param.fft_scratch_real[num]);
        doRefMapFFT(
            iRefMap, iOrient, iConvStart + iConv, comp_params[k + iConv].amp,
            comp_params[k + iConv].pha, comp_params[k + iConv].env,
            comp_params[k + iConv].sumC, comp_params[k + iConv].sumsquareC,
            param.fft_scratch_real[num], pProb, param.param_device, RefMap);
      }
1352
    }
Luka Stanisic's avatar
Luka Stanisic committed
1353
  }
1354
  else
Luka Stanisic's avatar
Luka Stanisic committed
1355
1356
1357
  {
    myblockCPU_t *comp_blocks = new myblockCPU_t[maxParallelConv];
    for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap++)
1358
1359
    {
#pragma omp parallel for schedule(dynamic, 1)
Luka Stanisic's avatar
Luka Stanisic committed
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
      for (int iConv = 0; iConv < maxParallelConv; iConv++)
      {
        const int num = omp_get_thread_num();
        calculateCCFFT(iRefMap, &localmultFFT[(k + iConv) * param.FFTMapSize],
                       param.fft_scratch_complex[num],
                       param.fft_scratch_real[num]);
        doRefMap_CPU_Parallel(iRefMap, iOrient, iConv,
                              param.fft_scratch_real[num], &comp_params[k],
                              comp_blocks);
      }
      doRefMap_CPU_Reduce(iRefMap, iOrient, iConvStart, maxParallelConv,
                          &comp_params[k], comp_blocks);
1372
    }
Luka Stanisic's avatar
Luka Stanisic committed
1373
1374
    delete[] comp_blocks;
  }
1375
1376

  cuda_custom_timeslot_end;
Luka Stanisic's avatar
Luka Stanisic committed
1377
  return (0);
1378
1379
}

Luka Stanisic's avatar
Luka Stanisic committed
1380
1381
inline void bioem::calculateCCFFT(int iRefMap, mycomplex_t *localConvFFT,
                                  mycomplex_t *localCCT, myfloat_t *lCC)
1382
{
1383
  //***************************************************************************************
Luka Stanisic's avatar
Luka Stanisic committed
1384
  //***** Calculating cross correlation with FFT algorithm *****
Pilar Cossio's avatar
Pilar Cossio committed
1385

Luka Stanisic's avatar
Luka Stanisic committed
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
  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;
  }

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

  myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
}
1405

Luka Stanisic's avatar
Luka Stanisic committed
1406
1407
1408
1409
1410
1411
1412
inline void bioem::doRefMap_CPU_Parallel(int iRefMap, int iOrient, int iConv,
                                         myfloat_t *lCC,
                                         myparam5_t *comp_params,
                                         myblockCPU_t *comp_block)
{
  //***************************************************************************************
  //***** Computation of log probabilities, done in parallel by OMP
1413

Luka Stanisic's avatar
Luka Stanisic committed
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
  int myGlobalId = iConv * param.param_device.NtotDisp;
  myfloat_t bestLogpro = MIN_PROB;
  int dispC =
      param.param_device.NumberPixels - param.param_device.maxDisplaceCenter;
  int cent_x, cent_y, address, bestId = 0;
  myfloat_t value, bestValue = 0.;
  myprob_t logpro = 0., sumExp = 0.;

  for (int myX = 0; myX < param.param_device.NxDisp; myX++)
  {
    for (int myY = 0; myY < param.param_device.NxDisp; myY++, myGlobalId++)
1425
    {
Luka Stanisic's avatar
Luka Stanisic committed
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
      cent_x = (myX * param.param_device.GridSpaceCenter + dispC) %
               param.param_device.NumberPixels;
      cent_y = (myY * param.param_device.GridSpaceCenter + dispC) %
               param.param_device.NumberPixels;
      address = cent_x * param.param_device.NumberPixels + cent_y;
      value = (myfloat_t) lCC[address] /
              (myfloat_t)(param.param_device.NumberPixels *
                          param.param_device.NumberPixels);

      logpro = calc_logpro(
          param.param_device, comp_params[iConv].amp, comp_params[iConv].pha,
          comp_params[iConv].env, comp_params[iConv].sumC,
          comp_params[iConv].sumsquareC, value, RefMap.sum_RefMap[iRefMap],
          RefMap.sumsquare_RefMap[iRefMap]);
#ifdef DEBUG_PROB
      printf("\t\t\tProb: iRefMap %d, iOrient %d, iConv %d, "
             "cent_x %d, cent_y %d, address %d, value %f, logpro %f\n",
             iRefMap, iOrient, iConv, cent_x, cent_y, address, value, logpro);
#endif
      if (bestLogpro < logpro)
      {
        sumExp *= exp(-logpro + bestLogpro);
        bestLogpro = logpro;
        bestId = myGlobalId;
        bestValue = value;
      }
      sumExp += exp(logpro - bestLogpro);
1453
    }
Luka Stanisic's avatar
Luka Stanisic committed
1454
  }
1455

Luka Stanisic's avatar
Luka Stanisic committed
1456
1457
1458
1459
1460
  comp_block[iConv].logpro = bestLogpro;
  comp_block[iConv].sumExp = sumExp;
  comp_block[iConv].id = bestId;
  comp_block[iConv].value = bestValue;
}
1461

Luka Stanisic's avatar
Luka Stanisic committed
1462
1463
1464
1465
1466
1467
1468
inline void bioem::doRefMap_CPU_Reduce(int iRefMap, int iOrient, int iConvStart,
                                       int maxParallelConv,
                                       myparam5_t *comp_params,
                                       myblockCPU_t *comp_block)
{
  //***************************************************************************************
  //***** Reduction of previously compututed log probabilities
Pilar Cossio's avatar
Pilar Cossio committed
1469

Luka Stanisic's avatar
Luka Stanisic committed
1470
  bioem_Probability_map &pProbMap = pProb.getProbMap(iRefMap);
1471

Luka Stanisic's avatar
Luka Stanisic committed
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
  for (int i = 0; i < maxParallelConv; i++)
  {
    if (pProbMap.Constoadd < comp_block[i].logpro)
    {
      pProbMap.Total *= exp(-comp_block[i].logpro + pProbMap.Constoadd);
      pProbMap.Constoadd = comp_block[i].logpro;

      // ********** Getting parameters that maximize the probability ***********
      int myGlobalId = comp_block[i].id;
      int myConv = myGlobalId / param.param_device.NtotDisp;
      myGlobalId -= myConv * param.param_device.NtotDisp;
      int myX = myGlobalId / param.param_device.NxDisp;
      myGlobalId -= myX * param.param_device.NxDisp;
      int myY = myGlobalId;

      int dispC = param.param_device.NumberPixels -
                  param.param_device.maxDisplaceCenter;
      myfloat_t value = comp_block[i].value;

      pProbMap.max.max_prob_cent_x =
          -((myX * param.param_device.GridSpaceCenter + dispC) -
            param.param_device.NumberPixels);
      pProbMap.max.max_prob_cent_y =
          -((myY * param.param_device.GridSpaceCenter + dispC) -
            param.param_device.NumberPixels);
      pProbMap.max.max_prob_orient = iOrient;
      pProbMap.max.max_prob_conv = iConvStart + myConv;

      pProbMap.max.max_prob_norm =
          -(-comp_params[myConv].sumC * RefMap.sum_RefMap[iRefMap] +
            param.param_device.Ntotpi * value) /
          (comp_params[myConv].sumC * comp_params[myConv].sumC -
           comp_params[myConv].sumsquareC * param.param_device.Ntotpi);
      pProbMap.max.max_prob_mu =
          -(-comp_params[myConv].sumC * value +
            comp_params[myConv].sumsquareC * RefMap.sum_RefMap[iRefMap]) /
          (comp_params[myConv].sumC * comp_params[myConv].sumC -
           comp_params[myConv].sumsquareC * param.param_device.Ntotpi);

#ifdef DEBUG_PROB
      printf("\tProbabilities change: iRefMap %d, iOrient %d, iConv %d, "
             "Total %f, Const %f, bestlogpro %f, sumExp %f, bestId %d\n",
             iRefMap, iOrient, iConvStart + myConv, pProbMap.Total,
             pProbMap.Constoadd, comp_block[i].logpro, comp_block[i].sumExp,
             comp_block[i].id);
      printf("\tParameters: myConv %d, myX %d, myY %d, cent_x -, cent_y -, "
             "probX %d, probY %d\n",
             myConv, myX, myY, pProbMap.max.max_prob_cent_x,
             pProbMap.max.max_prob_cent_y);
#endif
1522
    }
Luka Stanisic's avatar
Luka Stanisic committed
1523
1524
1525
1526
1527
1528
1529
    pProbMap.Total +=
        comp_block[i].sumExp * exp(comp_block[i].logpro - pProbMap.Constoadd);
#ifdef DEBUG_PROB
    printf("\t\tProbabilities after Reduce: iRefMap %d, iOrient %d, iConv "
           "%d, Total %f, Const %f, bestlogpro %f, sumExp %f, bestId %d\n",
           iRefMap, iOrient, iConvStart, pProbMap.Total, pProbMap.Constoadd,
           comp_block[i].logpro, comp_block[i].sumExp, comp_block[i].id);
Pilar Cossio's avatar
Pilar Cossio committed
1530
1531
#endif

Luka Stanisic's avatar
Luka Stanisic committed
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
    // Code for writing angles, not used by default
    if (param.param_device.writeAngles)
    {
      bioem_Probability_angle &pProbAngle =
          pProb.getProbAngle(iRefMap, iOrient);
      if (pProbAngle.ConstAngle < comp_block[i].logpro)
      {
        pProbAngle.forAngles *=
            exp(-comp_block[i].logpro + pProbAngle.ConstAngle);
        pProbAngle.ConstAngle = comp_block[i].logpro;
      }
      pProbAngle.forAngles += comp_block[i].sumExp *
                              exp(comp_block[i].logpro - pProbAngle.ConstAngle);
    }
  }
1547
}
1548

Luka Stanisic's avatar
Luka Stanisic committed
1549
int bioem::createProjection(int iMap, mycomplex_t *mapFFT)
1550
{
1551
  // **************************************************************************************
Luka Stanisic's avatar
Luka Stanisic committed
1552
1553
1554
1555
  // ****  BioEM Create Projection routine in Euler angles / Quaternions
  // ******************
  // ********************* and turns projection into Fourier space
  // ************************
1556
1557
  // **************************************************************************************

1558
  cuda_custom_timeslot("Projection", iMap, 0, COLOR_PROJECTION);
1559
1560
1561
1562

  myfloat3_t RotatedPointsModel[Model.nPointsModel];
  myfloat_t rotmat[3][3];
  myfloat_t alpha, gam, beta;
Luka Stanisic's avatar
Luka Stanisic committed
1563
  myfloat_t *localproj;
1564
1565

  localproj = param.fft_scratch_real[omp_get_thread_num()];
Luka Stanisic's avatar
Luka Stanisic committed
1566
1567
1568
  memset(localproj, 0, param.param_device.NumberPixels *
                           param.param_device.NumberPixels *
                           sizeof(*localproj));
1569

1570
1571
  //*************** Rotating the model ****************************
  //*************** Quaternions ****************************
Luka Stanisic's avatar
Luka Stanisic committed
1572
1573
  if (param.doquater)
  {
Pilar Cossio's avatar