Commit beebec5d authored by David Rohr's avatar David Rohr

fix segfault if rotated model point out of map, improve debug message and options, fix memory leaks

parent 05e7d589
......@@ -68,6 +68,7 @@ ostream& operator<<(ostream& os, const vector<T>& v)
bioem::bioem()
{
FFTAlgo = getenv("FFTALGO") == NULL ? 0 : atoi(getenv("FFTALGO"));
DebugOutput = getenv("BIOEM_DEBUG_OUTPUT") == NULL ? 2 : atoi(getenv("BIOEM_DEBUG_OUTPUT"));
}
bioem::~bioem()
......@@ -219,8 +220,9 @@ int bioem::configure(int ac, char* av[])
if (getenv("BIOEM_DEBUG_BREAK"))
{
param.nTotGridAngles = atoi(getenv("BIOEM_DEBUG_BREAK"));
param.nTotCTFs = atoi(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;
}
pProb.init(RefMap.ntotRefMap, param.nTotGridAngles, *this);
......@@ -308,41 +310,44 @@ int bioem::run()
HighResTimer timer;
printf("\tMain Loop (GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)²), Pixels %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);
if (DebugOutput >= 1) printf("\tMain Loop (GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)²), Pixels %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);
// printf("\tInner Loop Count (%d %d %d) %lld\n", param.param_device.maxDisplaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, (long long int) (param.param_device.NumberPixels * param.param_device.NumberPixels * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1) * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1)));
for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
{
// ***************************************************************************************
// ***** Creating Projection for given orientation and transforming to Fourier space *****
timer.ResetStart();
if (DebugOutput >= 1) timer.ResetStart();
createProjection(iOrient, proj_mapFFT);
printf("Time Projection %d: %f\n", iOrient, timer.GetCurrentElapsedTime());
if (DebugOutput >= 1) printf("Time Projection %d: %f\n", iOrient, timer.GetCurrentElapsedTime());
// ***************************************************************************************
// ***** **** Internal Loop over convolutions **** *****
for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
{
printf("\t\tConvolution %d %d\n", iOrient, iConv);
if (DebugOutput >= 2) printf("\t\tConvolution %d %d\n", iOrient, iConv);
// *** Calculating convolutions of projection map and crosscorrelations ***
timer.ResetStart();
if (DebugOutput >= 2) timer.ResetStart();
createConvolutedProjectionMap(iOrient, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
printf("Time Convolution %d %d: %f\n", iOrient, iConv, timer.GetCurrentElapsedTime());
if (DebugOutput >= 2) printf("Time Convolution %d %d: %f\n", iOrient, iConv, timer.GetCurrentElapsedTime());
// ***************************************************************************************
// *** Comparing each calculated convoluted map with all experimental maps ***
timer.ResetStart();
if (DebugOutput >= 2) timer.ResetStart();
compareRefMaps(iOrient, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
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("Time Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s)\n", iOrient, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000.);
if (DebugOutput >= 2)
{
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("Time Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s)\n", iOrient, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000.);
}
}
}
//deallocating fftw_complex vector
......@@ -516,6 +521,12 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
if (i < 0 || j < 0 || i >= param.param_device.NumberPixels || j >= param.param_device.NumberPixels)
{
if (DebugOutput >= 3) cout << "Model Point out of map: " << i << ", " << j << "\n";
continue;
}
localproj[i * param.param_device.NumberPixels + j] += Model.densityPointsModel[n] / Model.NormDen;
}
......
......@@ -450,6 +450,7 @@ int bioem_cuda::deviceFinishRun()
if (GPUWorkload < 100)
{
pProb.copyFrom(pProb_host, *this);
free_device_host(pProb_host->ptr);
delete[] pProb_host;
}
......
......@@ -55,6 +55,7 @@ protected:
int nProjectionMapsTotal; //Maps in total
int FFTAlgo;
int DebugOutput;
mycomplex_t** localCCT;
myfloat_t** lCC;
......
......@@ -16,6 +16,7 @@ typedef float myfloat_t;
#define myfftw_plan_dft_r2c_2d fftwf_plan_dft_r2c_2d
#define myfftw_plan_dft_c2r_2d fftwf_plan_dft_c2r_2d
#define myfftw_plan fftwf_plan
#define myfftw_cleanup fftwf_cleanup
#define MY_CUFFT_C2R CUFFT_C2R
#define mycufftExecC2R cufftExecC2R
#define mycuComplex_t cuComplex
......@@ -32,6 +33,7 @@ typedef double myfloat_t;
#define myfftw_plan_dft_r2c_2d fftw_plan_dft_r2c_2d
#define myfftw_plan_dft_c2r_2d fftw_plan_dft_c2r_2d
#define myfftw_plan fftw_plan
#define myfftw_cleanup fftw_cleanup
#define mycufftExecC2R cufftExecZ2D
#define mycuComplex_t cuDoubleComplex
#define MY_CUFFT_C2R CUFFT_Z2D
......
......@@ -2,7 +2,7 @@
< BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
......@@ -150,6 +150,7 @@ int bioem_model::readModel()
NormDen += densityPointsModel[numres];
numres++;
}
fclose(file);
nPointsModel = numres;
cout << "Protein structure read from Standard File\n";
}
......
......@@ -2,7 +2,7 @@
< BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
......@@ -66,7 +66,7 @@ int bioem_param::readParameters()
bool yesGSPEnv = false ;
bool yesGSPpha = false ;
bool yesMDC = false ;
bool yesGCen = false ;
bool yesGCen = false ;
ifstream input(fileinput);
if (!input.good())
......@@ -105,7 +105,7 @@ int bioem_param::readParameters()
token = strtok(NULL, " ");
param_device.NumberPixels = int(atoi(token));
if (param_device.NumberPixels < 0 ) { cout << "*** Error: Negative Number of Pixels "; exit(1);}
cout << "Number of Pixels " << param_device.NumberPixels << "\n";
cout << "Number of Pixels " << param_device.NumberPixels << "\n";
yesNumPix= true ;
}
else if (strcmp(token, "GRIDPOINTS_ALPHA") == 0)
......@@ -215,7 +215,7 @@ int bioem_param::readParameters()
}
input.close();
// Checks for ALL INPUT
// Checks for ALL INPUT
if( not ( yesPixSi ) ){ cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ; exit (1);};
if( not ( yesNumPix ) ){ cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ; exit (1);};
if( not ( yesGPal ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ; exit (1);};
......@@ -233,15 +233,15 @@ int bioem_param::readParameters()
if( not ( yesGCen ) ) { cout << "**** INPUT MISSING: Please provide PIXEL_GRID_CENTER \n" ; exit (1);};
//if only one grid point for PSF kernel:
if( (myfloat_t) numberGridPointsCTF_amp == 1 ) gridCTF_amp = startGridCTF_amp;
if( (myfloat_t) numberGridPointsCTF_phase == 1 ) gridCTF_phase = startGridCTF_phase;
if( (myfloat_t) numberGridPointsEnvelop == 1 ) gridEnvelop = startGridEnvelop;
//if only one grid point for PSF kernel:
if( (myfloat_t) numberGridPointsCTF_amp == 1 ) gridCTF_amp = startGridCTF_amp;
if( (myfloat_t) numberGridPointsCTF_phase == 1 ) gridCTF_phase = startGridCTF_phase;
if( (myfloat_t) numberGridPointsEnvelop == 1 ) gridEnvelop = startGridEnvelop;
//More checks with input parameters
// Envelope should not have a standard deviation greater than Npix/2
// Envelope should not have a standard deviation greater than Npix/2
if(sqrt(1./( (myfloat_t) numberGridPointsDisplaceCenter * gridEnvelop + startGridEnvelop))> float(param_device.NumberPixels)/2.0) {
cout << "MAX Standar deviation of envelope is larger than Allowed KERNEL Length";
cout << "MAX Standard deviation of envelope is larger than Allowed KERNEL Length\n";
exit(1);
}
// Envelop param should be positive
......@@ -255,7 +255,7 @@ int bioem_param::readParameters()
exit(1);
}
cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
cout << "Preparing FFTs\n";
......@@ -303,6 +303,7 @@ void bioem_param::releaseFFTPlans()
myfftw_destroy_plan(fft_plan_c2c_backward);
myfftw_destroy_plan(fft_plan_r2c_forward);
myfftw_destroy_plan(fft_plan_c2r_backward);
myfftw_cleanup();
}
fft_plans_created = 0;
}
......@@ -339,7 +340,7 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
}
nTotGridAngles = n;
// ********** Calculating normalized volumen element *********
param_device.volu = grid_alpha * grid_alpha * cos_grid_beta * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment