Commit 9768552f authored by Pilar Cossio's avatar Pilar Cossio
Browse files

Correcting bug accessing ilegal memory in CC mat & tidying up

parent 19bc6d4d
......@@ -6,7 +6,7 @@ option (INCLUDE_CUDA "Build BioEM with CUDA support" ON)
option (INCLUDE_OPENMP "Build BioEM with OpenMP support" ON)
option (INCLUDE_MPI "Build BioEM with MPI support" ON)
option (PRINT_CMAKE_VARIABLES "List all CMAKE Variables" OFF)
option (CUDA_FORCE_GCC "Force GCC as host compiler for CUDA part (If standard host compiler is incompatible with CUDA)" OFF)
option (CUDA_FORCE_GCC "Force GCC as host compiler for CUDA part (If standard host compiler is incompatible with CUDA)" ON)
###Set up general variables
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")
......@@ -18,9 +18,9 @@ set (BIOEM_ICC_FLAGS "-xHost -O3 -fno-alias -fno-fnalias -unroll -g0 -ipo")
set (BIOEM_GCC_FLAGS "-O3 -march=native -fweb -mfpmath=sse -frename-registers -minline-all-stringops -ftracer -funroll-loops -fpeel-loops -fprefetch-loop-arrays -ffast-math -ggdb")
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${BIOEM_ICC_FLAGS}")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${BIOEM_ICC_FLAGS}")
else()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${BIOEM_GCC_FLAGS}")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${BIOEM_GCC_FLAGS}")
endif()
set (BIOEM_SOURCE_FILES "bioem.cpp" "main.cpp" "map.cpp" "model.cpp" "param.cpp" "cmodules/timer.cpp")
......@@ -43,93 +43,94 @@ include_directories(${Boost_INCLUDE_DIRS})
###Find CUDA
set (BIOEM_CUDA_STATUS "Disabled")
if (INCLUDE_CUDA)
set (BIOEM_CUDA_STATUS "Not Found")
find_package(CUDA)
set (BIOEM_CUDA_STATUS "Not Found")
find_package(CUDA)
endif()
if (CUDA_FOUND)
if (CUDA_FORCE_GCC)
cmake_minimum_required(VERSION 2.8.10.1)
#Use GCC as host compiler for CUDA even though host compiler for other files is not GCC
set (CUDA_HOST_COMPILER gcc)
endif()
set(CUDA_NVCC_FLAGS "${CUDA_NVCC_FLAGS};--use_fast_math;-ftz=true;-O4;-Xptxas -O4")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_13,code=sm_13")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_20,code=sm_20")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_20,code=sm_21")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_30,code=sm_30")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_35,code=sm_35")
add_definitions(-DWITH_CUDA)
set (BIOEM_CUDA_STATUS "Found")
if (CUDA_FORCE_GCC)
cmake_minimum_required(VERSION 2.8.10.1)
#Use GCC as host compiler for CUDA even though host compiler for other files is not GCC
set (CUDA_HOST_COMPILER gcc)
endif()
set(CUDA_NVCC_FLAGS "${CUDA_NVCC_FLAGS};--use_fast_math;-ftz=true;-O4;-Xptxas -O4")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_13,code=sm_13")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_20,code=sm_20")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_20,code=sm_21")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_30,code=sm_30")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_35,code=sm_35")
list(APPEND CUDA_NVCC_FLAGS "-gencode=arch=compute_52,code=sm_52")
add_definitions(-DWITH_CUDA)
set (BIOEM_CUDA_STATUS "Found")
endif()
###Find OpenMP
set (BIOEM_OPENMP_STATUS "Disabled")
if (INCLUDE_OPENMP)
set (BIOEM_OPENMP_STATUS "Not Found")
find_package(OpenMP)
set (BIOEM_OPENMP_STATUS "Not Found")
find_package(OpenMP)
endif()
if(OPENMP_FOUND)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
add_definitions(-DWITH_OPENMP)
set (BIOEM_OPENMP_STATUS "Found")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
add_definitions(-DWITH_OPENMP)
set (BIOEM_OPENMP_STATUS "Found")
endif()
###Find MPI
set (BIOEM_MPI_STATUS "Disabled")
if (INCLUDE_MPI)
set (BIOEM_MPI_STATUS "Not Found")
find_package(MPI)
set (BIOEM_MPI_STATUS "Not Found")
find_package(MPI)
endif()
if (MPI_FOUND)
include_directories(${MPI_INCLUDE_PATH})
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${MPI_COMPILE_FLAGS} ${MPI_LINK_FLAGS}")
add_definitions(-DWITH_MPI)
set (BIOEM_MPI_STATUS "Found")
include_directories(${MPI_INCLUDE_PATH})
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${MPI_COMPILE_FLAGS} ${MPI_LINK_FLAGS}")
add_definitions(-DWITH_MPI)
set (BIOEM_MPI_STATUS "Found")
endif()
###Build Executable
if (CUDA_FOUND)
set(BIOEM_TMP_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS})
if (CUDA_FORCE_GCC)
#Hack to use GCC flags for GCC host compiler during NVCC compilation, although host compiler is in fact not GCC for other files
set(CMAKE_CXX_FLAGS ${BIOEM_GCC_FLAGS})
endif()
cuda_add_executable(bioEM ${BIOEM_SOURCE_FILES} bioem_cuda.cu)
set(CMAKE_CXX_FLAGS ${BIOEM_TMP_CMAKE_CXX_FLAGS})
set(BIOEM_TMP_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS})
if (CUDA_FORCE_GCC)
#Hack to use GCC flags for GCC host compiler during NVCC compilation, although host compiler is in fact not GCC for other files
set(CMAKE_CXX_FLAGS ${BIOEM_GCC_FLAGS})
endif()
cuda_add_executable(bioEM ${BIOEM_SOURCE_FILES} bioem_cuda.cu)
set(CMAKE_CXX_FLAGS ${BIOEM_TMP_CMAKE_CXX_FLAGS})
else()
add_executable(bioEM ${BIOEM_SOURCE_FILES})
add_executable(bioEM ${BIOEM_SOURCE_FILES})
endif()
#Additional CXX Flags not used by CUDA compiler
#Additional CXX Flags not used by CUDA compiler
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -Wall -pedantic")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -Wall -pedantic")
else()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -Wall -Wno-vla -Wno-unused-result -Wno-unused-local-typedefs -pedantic")
endif()
if (NOT OPENMP_FOUND)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unknown-pragmas")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-vla -Wno-long-long -Wall -pedantic")
endif()
###Add Libraries
if (CUDA_FOUND)
cuda_add_cufft_to_target(bioEM)
target_link_libraries(bioEM ${CUDA_CUDA_LIBRARY})
cuda_add_cufft_to_target(bioEM)
#set(CUDA_LIBRARIES "/afs/ipp/.cs/cuda/6.5-gtx9/amd64_sles11/lib64/stubs/libcuda.so")
#cuda_add_library()
target_link_libraries(bioEM ${CUDA_CUDA_LIBRARY})
# target_link_libraries(bioEM ${CUDA_LIBRARIES})
endif()
if (FFTWF_LIBRARIES)
target_link_libraries(bioEM ${FFTWF_LIBRARIES})
target_link_libraries(bioEM ${FFTWF_LIBRARIES})
else()
target_link_libraries(bioEM -L${FFTW_LIBDIR} -lfftw3 -lfftw3f)
target_link_libraries(bioEM -L${FFTW_LIBDIR} -lfftw3 -lfftw3f)
endif()
target_link_libraries(bioEM -L${Boost_LIBRARY_DIRS} -lboost_program_options)
if (MPI_FOUND)
target_link_libraries(bioEM ${MPI_LIBRARIES})
target_link_libraries(bioEM ${MPI_LIBRARIES})
endif()
###Show Status
......@@ -137,14 +138,16 @@ message(STATUS "Build Status")
message(STATUS "FFTW library: ${FFTW_LIBDIR}")
message(STATUS "Boost directory: ${Boost_LIBRARY_DIRS}")
message(STATUS "FFTW includedir: ${FFTW_INCLUDEDIR}")
message(STATUS "CUDA libraries: ${CUDA_CUDA_LIBRARY}")
message(STATUS "CUDA libraries: ${CUDA_LIBRARIES}")
message(STATUS "CUDA: ${BIOEM_CUDA_STATUS}")
message(STATUS "OpenMP: ${BIOEM_OPENMP_STATUS}")
message(STATUS "MPI: ${BIOEM_MPI_STATUS}")
if (PRINT_CMAKE_VARIABLES)
get_cmake_property(_variableNames VARIABLES)
foreach (_variableName ${_variableNames})
message(STATUS "${_variableName}=${${_variableName}}")
endforeach()
get_cmake_property(_variableNames VARIABLES)
foreach (_variableName ${_variableNames})
message(STATUS "${_variableName}=${${_variableName}}")
endforeach()
endif()
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
< BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer.
Copyright (C) 2016 Pilar Cossio, David Rohr and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution.
......@@ -104,7 +104,7 @@ int bioem::configure(int ac, char* av[])
HighResTimer timer;
std::string infile, modelfile, mapfile,Inputanglefile,Inputbestmap;
std::string infile, modelfile, mapfile,Inputanglefile,Inputbestmap;
if (mpi_rank == 0)
{
// *** Inizialzing default variables ***
......@@ -290,32 +290,34 @@ int bioem::configure(int ac, char* av[])
exit(1);
}
if(!param.printModel)param.CalculateGridsParam(Inputanglefile.c_str());
// Generating Grids of orientations
if(!param.printModel)param.CalculateGridsParam(Inputanglefile.c_str());
}
#ifdef WITH_MPI
// 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******************
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 !!!!!!!!!!!!!!
//We have to reinitialize all pointers !!!!!!!!!!!!
if (mpi_rank != 0) param.angprior = NULL;
// cout << "HERE MPI " << param.nTotGridAngles << " " << param.nTotGridAngles * sizeof(myfloat3_t*) << "\n";
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);
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);
//Allocating memory for prior
// if (mpi_rank != 0){//param.angprior = (myfloat_t*) mallocchk(3* param.nTotGridAngles * sizeof (myfloat_t*));
// delete[] param.angprior;
// param.angprior = new myfloat_t[param.NotUn_angles] ;
//
#ifdef PILAR_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";}
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";}
#endif
//refCtf, CtfParam, angles automatically filled by precalculare function below
//****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);
......@@ -328,6 +330,7 @@ int bioem::configure(int ac, char* av[])
}
#endif
// ****************** Precalculating Necessary Stuff *********************
if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
param.PrepareFFTs();
......@@ -339,6 +342,7 @@ int bioem::configure(int ac, char* av[])
}
precalculate();
// ****************** For debugging *********************
if (getenv("BIOEM_DEBUG_BREAK"))
{
const int cut = atoi(getenv("BIOEM_DEBUG_BREAK"));
......@@ -358,7 +362,11 @@ int bioem::configure(int ac, char* av[])
printf("Time Init Probabilities %f\n", timer.GetCurrentElapsedTime());
timer.ResetStart();
}
// ****************** Initializng pointers *********************
deviceInit();
if (DebugOutput >= 2 && mpi_rank == 0) printf("Time Device Init %f\n", timer.GetCurrentElapsedTime());
return(0);
......@@ -366,7 +374,7 @@ int bioem::configure(int ac, char* av[])
void bioem::cleanup()
{
//Deleting allocated pointers
//Deleting allocated pointers
free_device_host(pProb.ptr);
RefMap.freePointers();
}
......@@ -377,10 +385,6 @@ int bioem::precalculate()
// **Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels**
// **************************************************************************************
HighResTimer timer;
// Generating Grids of orientations
// param.CalculateGridsParam();
if (DebugOutput >= 3)
{
printf("\tTime Precalculate Grids Param: %f\n", timer.GetCurrentElapsedTime());
......@@ -477,13 +481,14 @@ int bioem::run()
cc++;
}
}
if(!FFTAlgo){cout << "Cross correlation calculation must be with enviormental variable FFTALGO=1\n"; exit(1);}
if(!FFTAlgo){cout << "Cross correlation calculation must be with enviormental variable FFTALGO=1\n"; exit(1);}
}
}
if(!FFTAlgo){cout << "Remark: Not using FFT algorithm. Not using Prior in B-Env.";}
// **************************************************************************************
deviceStartRun();
// ******************************** MAIN CYCLE ******************************************
......@@ -495,6 +500,7 @@ int bioem::run()
//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
//******** Alocating Vectors *************
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);
......@@ -504,13 +510,15 @@ int bioem::run()
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);
if(mpi_size > param.nTotGridAngles){
cout << "EXIT: Wrong MPI setup More MPI processes than orientations\n"; exit(1);
}
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;
// **************************Loop Over orientations***************************************
for (int iOrientAtOnce = iOrientStart; iOrientAtOnce < iOrientEnd; iOrientAtOnce += nProjectionsAtOnce)
{
// ***************************************************************************************
......@@ -518,6 +526,9 @@ int bioem::run()
if (DebugOutput >= 1) timer2.ResetStart();
if (DebugOutput >= 2) timer.ResetStart();
int iTmpEnd = std::min(iOrientEnd, iOrientAtOnce + nProjectionsAtOnce);
// **************************Parallel orientations for projections at once***************
#pragma omp parallel for
for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
{
......@@ -528,8 +539,10 @@ int bioem::run()
for (int iOrient = iOrientAtOnce; iOrient < iTmpEnd;iOrient++)
{
mycomplex_t* proj_mapFFT = &proj_mapsFFT[(iOrient - iOrientAtOnce) * ProjMapSize];
// ***************************************************************************************
// ***** **** Internal Loop over convolutions **** *****
// ***** **** Internal Loop over PSF/CTF convolutions **** *****
for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
{
// *** Calculating convolutions of projection map and crosscorrelations ***
......@@ -538,8 +551,7 @@ int bioem::run()
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);
// ***************************************************************************************
// *** Comparing each calculated convoluted map with all experimental maps ***
if (DebugOutput >= 2) timer.ResetStart();
myfloat_t amp,pha,env;
......@@ -547,6 +559,9 @@ int bioem::run()
pha=param.CtfParam[iConv].pos[1];
env=param.CtfParam[iConv].pos[2];
// ******************Internal loop over Reference images CUDA or OpenMP******************
// *** Comparing each calculated convoluted map with all experimental maps ***
compareRefMaps(iOrient, iConv, amp, pha, env, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
if (DebugOutput >= 2)
......@@ -576,9 +591,9 @@ int bioem::run()
deviceFinishRun();
// ************* Writing Out Probabilities ***************
// *** Angular Probability ***
// ************* Collecing all the probabilities from MPI replicas ***************
#ifdef WITH_MPI
if (mpi_size > 1)
......@@ -598,7 +613,11 @@ int bioem::run()
for (int i = 0;i < RefMap.ntotRefMap;i++)
{
bioem_Probability_map& pProbMap = pProb.getProbMap(i);
#ifdef PILAR_DEBUG
cout << "Reduction " << mpi_rank << " Map " << i << " Prob " << pProbMap.Total << " Const " << pProbMap.Constoadd << "\n";
#endif
tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
}
MPI_Reduce(tmp1, tmp3, RefMap.ntotRefMap, MY_MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD);
......@@ -661,13 +680,13 @@ int bioem::run()
myfloat_t* tmp3 = new myfloat_t[count];
for (int i = 0;i < RefMap.ntotRefMap;i++)
{
for (int j = 0;j < param.nTotGridAngles;j++)
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;
// 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++)
......@@ -698,6 +717,8 @@ int bioem::run()
}
#endif
// ************* Writing Out Probabilities ***************
if (mpi_rank == 0)
{
......@@ -710,8 +731,8 @@ int bioem::run()
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";
// angProbfile <<"Model Used: " << modelfile.c_str() << "\n";
// angProbfile <<"Input Used: " << infile.c_str() << "\n";
}
// Output for Cross Correlation File
ofstream ccProbfile;
......@@ -746,6 +767,7 @@ int bioem::run()
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";
// Loop over reference maps
// ************* Over all maps ***************
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
......@@ -757,27 +779,34 @@ int bioem::run()
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)) << " ";
// *** 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";
}else{
outputProbFile << "PROBLEM!! with Map " << iRefMap << "Numerical Integrated Probability = 0.0; Options for error: \n - Try re-normalizing experimental map \n -looking at your model density. \n - Maybe radius is less than pixel size (If so try KEYWORD NO_PROJECT_RADIUS in parameter inputfile)\n";
outputProbFile << "RefMap: " << iRefMap << "Most probabily Numerical Constant Out of Float-point range: " << pProbMap.Constoadd << "\n"; }
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";
// Writing out CTF parameters if requiered
if(param.writeCTF && param.usepsf){
......@@ -792,28 +821,40 @@ int bioem::run()
outputProbFile << 4*M_PI*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi << " [A²] \n";
}
// Writing Individual Angle probabilities
//*************** Writing Individual Angle probabilities
if(param.param_device.writeAngles)
{
for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
{
bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
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.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";
}
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";
}
}
}
}
// Writing Cross-Correlations if requiered
//************* Writing Cross-Correlations if requiered
if(param.param_device.writeCC){
int cc=0;
......@@ -872,11 +913,11 @@ int bioem::run()
{
for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry = ry + param.param_device.CCdisplace)
{
if(localcc[ rx * param.param_device.NumberPixels + ry ] > -FLT_MAX){
if(localcc[ rx * param.param_device.NumberPixels + ry ] > -FLT_MAX){
ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
}else{
}else{
ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << -FLT_MAX << "\n" ;
}
}
// cout << " cc " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] <<"\n" ;
}
ccProbfile << "\n";
......@@ -886,11 +927,11 @@ int bioem::run()
{
for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry++)
{
if(localcc[ rx * param.param_device.NumberPixels + ry ] > -FLT_MAX){
if(localcc[ rx * param.param_device.NumberPixels + ry ] > -FLT_MAX){
ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
}else{
}else{
ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << -FLT_MAX << "\n" ;
}
}
}
ccProbfile << "\n";
}
......@@ -947,6 +988,12 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_
//***************************************************************************************
//***** Calculating cross correlation in FFTALGOrithm *****
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++)
{
......@@ -955,9 +1002,12 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_
}
myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
// 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);
doRefMapFFT(iRefMap, iOrient, iConv, amp, pha, env, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
#ifdef PILAR_DEBUG
if (param.param_device.writeCC)