From c9754f9e8d492633c6320dc4cc080fc64e41f05c Mon Sep 17 00:00:00 2001
From: qon <qon@jwdt.org>
Date: Sun, 6 Apr 2014 22:25:31 +0200
Subject: [PATCH] Change all line endings in files to unix style
---
bioem.cpp | 1094 +++++++++++++++++++++++------------------------
bioem_cuda.cu | 420 +++++++++---------
include/bioem.h | 96 ++---
include/defs.h | 84 ++--
include/map.h | 176 ++++----
include/model.h | 54 +--
main.cpp | 118 ++---
map.cpp | 232 +++++-----
model.cpp | 472 ++++++++++----------
9 files changed, 1373 insertions(+), 1373 deletions(-)
diff --git a/bioem.cpp b/bioem.cpp
index a552886..48642d5 100644
--- a/bioem.cpp
+++ b/bioem.cpp
@@ -1,547 +1,547 @@
-#include <fstream>
-#include <boost/program_options.hpp>
-#include <iostream>
-#include <algorithm>
-#include <iterator>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string>
-#include <cmath>
-#include <omp.h>
-
-#include <fftw3.h>
-#include <math.h>
-#include "cmodules/timer.h"
-
-#include "param.h"
-#include "bioem.h"
-#include "model.h"
-#include "map.h"
-
-
-#include "bioem_algorithm.h"
-
-
-using namespace boost;
-namespace po = boost::program_options;
-
-using namespace std;
-
-// A helper function of Boost
-template<class T>
-ostream& operator<<(ostream& os, const vector<T>& v)
-{
- copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
- return os;
-}
-
-bioem::bioem()
-{
-
-}
-
-bioem::~bioem()
-{
-
-}
-
-int bioem::configure(int ac, char* av[])
-{
- /**************************************************************************************/
- /**** Configuration Routine using boost for extracting parameters, models and maps ****/
- /**************************************************************************************/
- /****** And Precalculating necessary grids, map crosscorrelations and kernels ********/
- /*************************************************************************************/
-
- /*** Inizialzing default variables ***/
- std::string infile,modelfile,mapfile;
- Model.readPDB=false;
- param.writeAngles=false;
- RefMap.dumpMap = false;
- RefMap.loadMap = false;
-
- /*************************************************************************************/
- cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
- /*************************************************************************************/
-
- /********************* Command line reading input with BOOST ************************/
-
- try {
- po::options_description desc("Command line inputs");
- desc.add_options()
- ("Inputfile", po::value<std::string>(), "Name of input parameter file")
- ("Modelfile", po::value< std::string>() , "Name of model file")
- ("Particlesfile", po::value< std::string>(), "Name of paricles file")
- ("ReadPDB", "(Optional) If reading model file in PDB format")
- ("DumpMaps", "(Optional) Dump maps after they were red from maps file")
- ("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
- ("help", "(Optional) Produce help message")
- ;
-
- po::positional_options_description p;
- p.add("Inputfile", -1);
- p.add("Modelfile", -1);
- p.add("Particlesfile", -1);
- p.add("ReadPDB", -1);
- p.add("DumpMaps", -1);
- p.add("LoadMapDump", -1);
-
- po::variables_map vm;
- po::store(po::command_line_parser(ac, av).
- options(desc).positional(p).run(), vm);
- po::notify(vm);
-
- if((ac < 6)) {
- std::cout << desc << std::endl;
- return 0;
- }
- if (vm.count("help")) {
- cout << "Usage: options_description [options]\n";
- cout << desc;
- return 0;
- }
-
- if (vm.count("Inputfile"))
- {
- cout << "Input file is: ";
- cout << vm["Inputfile"].as< std::string >()<< "\n";
- infile=vm["Inputfile"].as< std::string >();
- }
- if (vm.count("Modelfile"))
- {
- cout << "Model file is: "
- << vm["Modelfile"].as< std::string >() << "\n";
- modelfile=vm["Modelfile"].as< std::string >();
- }
-
- if (vm.count("ReadPDB"))
- {
- cout << "Reading model file in PDB format.\n";
- Model.readPDB=true;
- }
-
- if (vm.count("DumpMaps"))
- {
- cout << "Dumping Maps after reading from file.\n";
- RefMap.dumpMap = true;
- }
-
- if (vm.count("LoadMapDump"))
- {
- cout << "Loading Map dump.\n";
- RefMap.loadMap = true;
- }
-
- if (vm.count("Particlesfile"))
- {
- cout << "Paricle file is: "
- << vm["Particlesfile"].as< std::string >() << "\n";
- mapfile=vm["Particlesfile"].as< std::string >();
- }
- }
- catch(std::exception& e)
- {
- cout << e.what() << "\n";
- return 1;
- }
-
- /********************* Reading Parameter Input ***************************/
- // copying inputfile to param class
- param.fileinput = infile.c_str();
- param.readParameters();
-
- /********************* Reading Model Input ******************************/
- // copying modelfile to model class
- Model.filemodel = modelfile.c_str();
- Model.readModel();
-
- /********************* Reading Particle Maps Input **********************/
- /********* HERE: PROBLEM if maps dont fit on the memory!! ***************/
- // copying mapfile to ref map class
- RefMap.filemap = mapfile.c_str();
- RefMap.readRefMaps();
-
- /****************** Precalculating Necessary Stuff *********************/
- precalculate();
-
- param.nTotGridAngles = 10;
- param.nTotCTFs = 10;
- //param.param_device.maxDisplaceCenter = 0;
-
- deviceInit();
-
- return(0);
-}
-
-int bioem::precalculate()
-{
- /**************************************************************************************/
- /* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels */
- /**************************************************************************************/
-
- // Generating Grids of orientations
- param.CalculateGridsParam();
-
- //Inizialzing crosscorrelations of Maps
- memset(RefMap.sum_RefMap, 0, BIOEM_MAX_MAPS * sizeof(*RefMap.sum_RefMap));
- memset(RefMap.sumsquare_RefMap, 0, BIOEM_MAX_MAPS * sizeof(*RefMap.sum_RefMap));
-
- myfloat_t sum,sumsquare;
-
- //Precalculating cross-correlations of maps
- for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap ; iRefMap++)
- {
- calcross_cor(RefMap.Ref[iRefMap],sum,sumsquare);
- //Storing Crosscorrelations in Map class
- RefMap.sum_RefMap[iRefMap]=sum;
- RefMap.sumsquare_RefMap[iRefMap]=sumsquare;
- }
-
- // Precalculating CTF Kernels stored in class Param
- param.CalculateRefCTF();
-
- return(0);
-}
-
-
-int bioem::run()
-{
- /**************************************************************************************/
- /**** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****/
- /**************************************************************************************/
-
- /**** If we want to control the number of threads -> omp_set_num_threads(XX); ******/
- /****************** Declarying class of Probability Pointer *************************/
- pProb = new bioem_Probability[RefMap.ntotRefMap];
-
- printf("\tInitializing\n");
- // Inizialzing Probabilites to zero and constant to -Infinity
- for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
- {
- pProb[iRefMap].Total=0.0;
- pProb[iRefMap].Constoadd=-9999999;
- pProb[iRefMap].max_prob=-9999999;
- for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
- {
- pProb[iRefMap].forAngles[iOrient]=0.0;
- pProb[iRefMap].ConstAngle[iOrient]=-99999999;
- }
- }
- /**************************************************************************************/
- deviceStartRun();
-
- /******************************** MAIN CYCLE ******************************************/
-
- /*** Declaring Private variables for each thread *****/
- mycomplex_t* proj_mapFFT;
- bioem_map conv_map;
-
- //allocating fftw_complex vector
- proj_mapFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t) *4*param.param_device.NumberPixels*param.param_device.NumberPixels);
-
- 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);
- 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)));
- //#pragma omp parallel for
- for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
- {
- /***************************************************************************************/
- /***** Creating Projection for given orientation and transforming to Fourier space *****/
- timer.ResetStart();
- createProjection(iProjectionOut, proj_mapFFT);
- printf("Time Projection %d: %f\n", iProjectionOut, timer.GetCurrentElapsedTime());
-
- /***************************************************************************************/
- /***** **** Internal Loop over convolutions **** *****/
- for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
- {
- printf("\t\tConvolution %d %d\n", iProjectionOut, iConv);
- /*** Calculating convolutions of projection map and crosscorrelations ***/
- timer.ResetStart();
- createConvolutedProjectionMap(iProjectionOut,iConv,proj_mapFFT,conv_map);
- printf("Time Convolution %d %d: %f\n", iProjectionOut, iConv, timer.GetCurrentElapsedTime());
-
- /***************************************************************************************/
- /*** Comparing each calculated convoluted map with all experimental maps ***/
- timer.ResetStart();
- compareRefMaps(iProjectionOut, iConv, conv_map);
- 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", iProjectionOut, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000.);
- }
-
- }
- //deallocating fftw_complex vector
- fftw_free(proj_mapFFT);
-
- deviceFinishRun();
-
- /************* Writing Out Probabilities ***************/
-
- /*** Angular Probability ***/
-
- // if(param.writeAngles){
- ofstream angProbfile;
- angProbfile.open ("ANG_PROB_iRefMap");
- // }
-
- ofstream outputProbFile;
- outputProbFile.open ("Output_Probabilities");
-
- for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
- {
- /**** Total Probability ***/
- outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProb[iRefMap].Total)+pProb[iRefMap].Constoadd+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd << "\n";
-
- outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
-
- /*** Param that maximize probability****/
- outputProbFile << (pProb[iRefMap].max_prob + 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[pProb[iRefMap].max_prob_orient].pos[0] << " ";
- outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[1] << " ";
- outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[2] << " ";
- outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[0] << " ";
- outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[1] << " ";
- outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[2] << " ";
- outputProbFile << pProb[iRefMap].max_prob_cent_x << " ";
- outputProbFile << pProb[iRefMap].max_prob_cent_y;
- outputProbFile << "\n";
-
- /*** For individual files***/ //angProbfile.open ("ANG_PROB_"iRefMap);
-
- if(param.writeAngles)
- {
- for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
- {
- angProbfile << " " << iRefMap << " " << param.angles[iProjectionOut].pos[0] << " " << param.angles[iProjectionOut].pos[1] << " " << param.angles[iProjectionOut].pos[2] << " " << log(pProb[iRefMap].forAngles[iProjectionOut])+pProb[iRefMap].ConstAngle[iProjectionOut]+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
-
- }
- }
-
- }
-
- angProbfile.close();
- outputProbFile.close();
-
- //Deleting allocated pointers
-
- if (pProb)
- {
- delete[] pProb;
- pProb = NULL;
- }
-
- if (param.refCTF)
- {
- delete[] param.refCTF;
- param.refCTF =NULL;
- }
-
- return(0);
-}
-
-int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map)
-{
-#pragma omp parallel for
- for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
- {
- compareRefMapShifted<-1>(iRefMap,iProjectionOut,iConv,conv_map, pProb, param.param_device, RefMap);
- }
- return(0);
-}
-
-int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
-{
-
- /**************************************************************************************/
- /**** BioEM Create Projection routine in Euler angle predefined grid****************
- ********************* and turns projection into Fourier space **********************/
- /**************************************************************************************/
-
- myfloat3_t RotatedPointsModel[Model.nPointsModel];
- myfloat_t rotmat[3][3];
- myfloat_t alpha, gam,beta;
- fftw_plan plan;
- mycomplex_t* localproj;
- int totnumPixFFT=2*param.param_device.NumberPixels;
-
- localproj= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t) *4*param.param_device.NumberPixels*param.param_device.NumberPixels);
- memset(localproj,0,4*param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localproj));
-
- alpha=param.angles[iMap].pos[0];
- beta=param.angles[iMap].pos[1];
- gam=param.angles[iMap].pos[2];
-
- /**** To see how things are going: cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***/
-
- /********** Creat Rotation with pre-defiend grid of orientations**********/
-
- rotmat[0][0]=cos(gam)*cos(alpha)-cos(beta)*sin(alpha)*sin(gam);
- rotmat[0][1]=cos(gam)*sin(alpha)+cos(beta)*cos(alpha)*sin(gam);
- rotmat[0][2]=sin(gam)*sin(beta);
- rotmat[1][0]=-sin(gam)*cos(alpha)-cos(beta)*sin(alpha)*cos(gam);
- rotmat[1][1]=-sin(gam)*sin(alpha)+cos(beta)*cos(alpha)*cos(gam);
- rotmat[1][2]=cos(gam)*sin(beta);
- rotmat[2][0]=sin(beta)*sin(alpha);
- rotmat[2][1]=-sin(beta)*cos(alpha);
- rotmat[2][2]=cos(beta);
-
-
- for(int n=0; n< Model.nPointsModel; n++)
- {
- RotatedPointsModel[n].pos[0]=0.0;
- RotatedPointsModel[n].pos[1]=0.0;
- RotatedPointsModel[n].pos[2]=0.0;
- }
- for(int n=0; n< Model.nPointsModel; n++)
- {
- for(int k=0; k< 3; k++)
- {
- for(int j=0; j< 3; j++)
- {
- RotatedPointsModel[n].pos[k]+=rotmat[k][j]*Model.PointsModel[n].pos[j];
- }
- }
- }
-
- int i, j;
-
- /************ Projection over the Z axis********************/
- for(int n=0; n< Model.nPointsModel; n++)
- {
- //Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
- i=floor(RotatedPointsModel[n].pos[0]/(*param.pixelSize)+float(param.param_device.NumberPixels)/2.0+0.5);
- j=floor(RotatedPointsModel[n].pos[1]/(*param.pixelSize)+float(param.param_device.NumberPixels)/2.0+0.5);
-
- localproj[i*2*param.param_device.NumberPixels+j+param.param_device.NumberPixels*param.param_device.NumberPixels+int(param.param_device.NumberPixels/2.0)][0]+=Model.densityPointsModel[n];
-
- }
-
- /**** Output Just to check****/
- if(iMap==10)
- {
- ofstream myexamplemap;
- ofstream myexampleRot;
- myexamplemap.open ("MAP_i10");
- myexampleRot.open ("Rot_i10");
- myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
- for(int k=0; k<2*param.param_device.NumberPixels; k++)
- {
- for(int j=0; j<2*param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*2*param.param_device.NumberPixels+j][0];
- }
- myexamplemap << " \n";
- for(int n=0; n< Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
- myexamplemap.close();
- myexampleRot.close();
- }
-
- /***** Converting projection to Fourier Space for Convolution later with kernel****/
- /********** Omp Critical is necessary with FFTW*******/
- //#pragma omp critical
- {
- plan = fftw_plan_dft_2d(totnumPixFFT,totnumPixFFT,localproj,mapFFT,FFTW_FORWARD,FFTW_ESTIMATE);
- fftw_execute(plan);
- fftw_destroy_plan(plan);
- fftw_free(localproj);
- }
-
- return(0);
-}
-
-
-int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv)
-{
- /**************************************************************************************/
- /**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
- **************** calculated Projection with convoluted precalculated Kernel**********
- *************** and Backtransforming it to real Space ******************************/
- /**************************************************************************************/
-
- fftw_plan plan;
- mycomplex_t* localmultFFT;
- mycomplex_t* localconvFFT;
- int totnumPixFFT=2*param.param_device.NumberPixels;
- localmultFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t)*totnumPixFFT*totnumPixFFT);
- localconvFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t)*totnumPixFFT*totnumPixFFT);
-
-
- /**** Multiplying FFTmap with corresponding kernel ****/
-
- for(int i=0; i < 2*param.param_device.NumberPixels ; i++ )
- {
- for(int j=0; j < 2*param.param_device.NumberPixels ; j++ )
- { //Projection*CONJ(KERNEL)
- localmultFFT[i*2*param.param_device.NumberPixels+j][0]=lproj[i*2*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0]+lproj[i*2*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1];
- localmultFFT[i*2*param.param_device.NumberPixels+j][1]=lproj[i*2*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0]-lproj[i*2*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1];
- // cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0] << " " <<param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1] <<" " <<lproj[i*2*param.param_device.NumberPixels+j][0] <<" " <<lproj[i*2*param.param_device.NumberPixels+j][1] << "\n";
- }
- }
-
- /**** Bringing convoluted Map to real Space ****/
- //#pragma omp critical
- {
- plan = fftw_plan_dft_2d(totnumPixFFT,totnumPixFFT,localmultFFT,localconvFFT,FFTW_BACKWARD,FFTW_ESTIMATE);
- fftw_execute(plan);
- }
-
-
- /****Asigning convolution fftw_complex to bioem_map ****/
- for(int i=0; i < param.param_device.NumberPixels ; i++ )
- {
- for(int j=0; j < param.param_device.NumberPixels ; j++ )
- {
- Mapconv.points[i][j]=localconvFFT[i*2*param.param_device.NumberPixels+j+param.param_device.NumberPixels*param.param_device.NumberPixels+int(param.param_device.NumberPixels/2.0)][0];
- }
- }
-
- /**** Freeing fftw_complex created (dont know if omp critical is necessary) ****/
- //#pragma omp critical
- {
- fftw_destroy_plan(plan);
- fftw_free(localconvFFT);
- fftw_free(localmultFFT);
- }
-
- return(0);
-}
-
-int bioem::calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare)
-{
- /*********************** Routine to calculate Cross correlations***********************/
-
- sum=0.0;
- sumsquare=0.0;
- for (int i = 0; i < param.param_device.NumberPixels; i++)
- {
- for (int j = 0; j < param.param_device.NumberPixels; j++)
- {
- // Calculate Sum of pixels
- sum += localmap.points[i][j];
- // Calculate Sum of pixels squared
- sumsquare += localmap.points[i][j]*localmap.points[i][j];
- }
- }
- return(0);
-}
-
-int bioem::deviceInit()
-{
- return(0);
-}
-
-int bioem::deviceStartRun()
-{
- return(0);
-}
-
-int bioem::deviceFinishRun()
-{
- return(0);
-}
+#include <fstream>
+#include <boost/program_options.hpp>
+#include <iostream>
+#include <algorithm>
+#include <iterator>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <cmath>
+#include <omp.h>
+
+#include <fftw3.h>
+#include <math.h>
+#include "cmodules/timer.h"
+
+#include "param.h"
+#include "bioem.h"
+#include "model.h"
+#include "map.h"
+
+
+#include "bioem_algorithm.h"
+
+
+using namespace boost;
+namespace po = boost::program_options;
+
+using namespace std;
+
+// A helper function of Boost
+template<class T>
+ostream& operator<<(ostream& os, const vector<T>& v)
+{
+ copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
+ return os;
+}
+
+bioem::bioem()
+{
+
+}
+
+bioem::~bioem()
+{
+
+}
+
+int bioem::configure(int ac, char* av[])
+{
+ /**************************************************************************************/
+ /**** Configuration Routine using boost for extracting parameters, models and maps ****/
+ /**************************************************************************************/
+ /****** And Precalculating necessary grids, map crosscorrelations and kernels ********/
+ /*************************************************************************************/
+
+ /*** Inizialzing default variables ***/
+ std::string infile,modelfile,mapfile;
+ Model.readPDB=false;
+ param.writeAngles=false;
+ RefMap.dumpMap = false;
+ RefMap.loadMap = false;
+
+ /*************************************************************************************/
+ cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
+ /*************************************************************************************/
+
+ /********************* Command line reading input with BOOST ************************/
+
+ try {
+ po::options_description desc("Command line inputs");
+ desc.add_options()
+ ("Inputfile", po::value<std::string>(), "Name of input parameter file")
+ ("Modelfile", po::value< std::string>() , "Name of model file")
+ ("Particlesfile", po::value< std::string>(), "Name of paricles file")
+ ("ReadPDB", "(Optional) If reading model file in PDB format")
+ ("DumpMaps", "(Optional) Dump maps after they were red from maps file")
+ ("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
+ ("help", "(Optional) Produce help message")
+ ;
+
+ po::positional_options_description p;
+ p.add("Inputfile", -1);
+ p.add("Modelfile", -1);
+ p.add("Particlesfile", -1);
+ p.add("ReadPDB", -1);
+ p.add("DumpMaps", -1);
+ p.add("LoadMapDump", -1);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(ac, av).
+ options(desc).positional(p).run(), vm);
+ po::notify(vm);
+
+ if((ac < 6)) {
+ std::cout << desc << std::endl;
+ return 0;
+ }
+ if (vm.count("help")) {
+ cout << "Usage: options_description [options]\n";
+ cout << desc;
+ return 0;
+ }
+
+ if (vm.count("Inputfile"))
+ {
+ cout << "Input file is: ";
+ cout << vm["Inputfile"].as< std::string >()<< "\n";
+ infile=vm["Inputfile"].as< std::string >();
+ }
+ if (vm.count("Modelfile"))
+ {
+ cout << "Model file is: "
+ << vm["Modelfile"].as< std::string >() << "\n";
+ modelfile=vm["Modelfile"].as< std::string >();
+ }
+
+ if (vm.count("ReadPDB"))
+ {
+ cout << "Reading model file in PDB format.\n";
+ Model.readPDB=true;
+ }
+
+ if (vm.count("DumpMaps"))
+ {
+ cout << "Dumping Maps after reading from file.\n";
+ RefMap.dumpMap = true;
+ }
+
+ if (vm.count("LoadMapDump"))
+ {
+ cout << "Loading Map dump.\n";
+ RefMap.loadMap = true;
+ }
+
+ if (vm.count("Particlesfile"))
+ {
+ cout << "Paricle file is: "
+ << vm["Particlesfile"].as< std::string >() << "\n";
+ mapfile=vm["Particlesfile"].as< std::string >();
+ }
+ }
+ catch(std::exception& e)
+ {
+ cout << e.what() << "\n";
+ return 1;
+ }
+
+ /********************* Reading Parameter Input ***************************/
+ // copying inputfile to param class
+ param.fileinput = infile.c_str();
+ param.readParameters();
+
+ /********************* Reading Model Input ******************************/
+ // copying modelfile to model class
+ Model.filemodel = modelfile.c_str();
+ Model.readModel();
+
+ /********************* Reading Particle Maps Input **********************/
+ /********* HERE: PROBLEM if maps dont fit on the memory!! ***************/
+ // copying mapfile to ref map class
+ RefMap.filemap = mapfile.c_str();
+ RefMap.readRefMaps();
+
+ /****************** Precalculating Necessary Stuff *********************/
+ precalculate();
+
+ param.nTotGridAngles = 10;
+ param.nTotCTFs = 10;
+ //param.param_device.maxDisplaceCenter = 0;
+
+ deviceInit();
+
+ return(0);
+}
+
+int bioem::precalculate()
+{
+ /**************************************************************************************/
+ /* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels */
+ /**************************************************************************************/
+
+ // Generating Grids of orientations
+ param.CalculateGridsParam();
+
+ //Inizialzing crosscorrelations of Maps
+ memset(RefMap.sum_RefMap, 0, BIOEM_MAX_MAPS * sizeof(*RefMap.sum_RefMap));
+ memset(RefMap.sumsquare_RefMap, 0, BIOEM_MAX_MAPS * sizeof(*RefMap.sum_RefMap));
+
+ myfloat_t sum,sumsquare;
+
+ //Precalculating cross-correlations of maps
+ for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap ; iRefMap++)
+ {
+ calcross_cor(RefMap.Ref[iRefMap],sum,sumsquare);
+ //Storing Crosscorrelations in Map class
+ RefMap.sum_RefMap[iRefMap]=sum;
+ RefMap.sumsquare_RefMap[iRefMap]=sumsquare;
+ }
+
+ // Precalculating CTF Kernels stored in class Param
+ param.CalculateRefCTF();
+
+ return(0);
+}
+
+
+int bioem::run()
+{
+ /**************************************************************************************/
+ /**** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****/
+ /**************************************************************************************/
+
+ /**** If we want to control the number of threads -> omp_set_num_threads(XX); ******/
+ /****************** Declarying class of Probability Pointer *************************/
+ pProb = new bioem_Probability[RefMap.ntotRefMap];
+
+ printf("\tInitializing\n");
+ // Inizialzing Probabilites to zero and constant to -Infinity
+ for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
+ {
+ pProb[iRefMap].Total=0.0;
+ pProb[iRefMap].Constoadd=-9999999;
+ pProb[iRefMap].max_prob=-9999999;
+ for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
+ {
+ pProb[iRefMap].forAngles[iOrient]=0.0;
+ pProb[iRefMap].ConstAngle[iOrient]=-99999999;
+ }
+ }
+ /**************************************************************************************/
+ deviceStartRun();
+
+ /******************************** MAIN CYCLE ******************************************/
+
+ /*** Declaring Private variables for each thread *****/
+ mycomplex_t* proj_mapFFT;
+ bioem_map conv_map;
+
+ //allocating fftw_complex vector
+ proj_mapFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t) *4*param.param_device.NumberPixels*param.param_device.NumberPixels);
+
+ 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);
+ 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)));
+ //#pragma omp parallel for
+ for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
+ {
+ /***************************************************************************************/
+ /***** Creating Projection for given orientation and transforming to Fourier space *****/
+ timer.ResetStart();
+ createProjection(iProjectionOut, proj_mapFFT);
+ printf("Time Projection %d: %f\n", iProjectionOut, timer.GetCurrentElapsedTime());
+
+ /***************************************************************************************/
+ /***** **** Internal Loop over convolutions **** *****/
+ for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
+ {
+ printf("\t\tConvolution %d %d\n", iProjectionOut, iConv);
+ /*** Calculating convolutions of projection map and crosscorrelations ***/
+ timer.ResetStart();
+ createConvolutedProjectionMap(iProjectionOut,iConv,proj_mapFFT,conv_map);
+ printf("Time Convolution %d %d: %f\n", iProjectionOut, iConv, timer.GetCurrentElapsedTime());
+
+ /***************************************************************************************/
+ /*** Comparing each calculated convoluted map with all experimental maps ***/
+ timer.ResetStart();
+ compareRefMaps(iProjectionOut, iConv, conv_map);
+ 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", iProjectionOut, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000.);
+ }
+
+ }
+ //deallocating fftw_complex vector
+ fftw_free(proj_mapFFT);
+
+ deviceFinishRun();
+
+ /************* Writing Out Probabilities ***************/
+
+ /*** Angular Probability ***/
+
+ // if(param.writeAngles){
+ ofstream angProbfile;
+ angProbfile.open ("ANG_PROB_iRefMap");
+ // }
+
+ ofstream outputProbFile;
+ outputProbFile.open ("Output_Probabilities");
+
+ for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
+ {
+ /**** Total Probability ***/
+ outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProb[iRefMap].Total)+pProb[iRefMap].Constoadd+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd << "\n";
+
+ outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
+
+ /*** Param that maximize probability****/
+ outputProbFile << (pProb[iRefMap].max_prob + 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[pProb[iRefMap].max_prob_orient].pos[0] << " ";
+ outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[1] << " ";
+ outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[2] << " ";
+ outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[0] << " ";
+ outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[1] << " ";
+ outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[2] << " ";
+ outputProbFile << pProb[iRefMap].max_prob_cent_x << " ";
+ outputProbFile << pProb[iRefMap].max_prob_cent_y;
+ outputProbFile << "\n";
+
+ /*** For individual files***/ //angProbfile.open ("ANG_PROB_"iRefMap);
+
+ if(param.writeAngles)
+ {
+ for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
+ {
+ angProbfile << " " << iRefMap << " " << param.angles[iProjectionOut].pos[0] << " " << param.angles[iProjectionOut].pos[1] << " " << param.angles[iProjectionOut].pos[2] << " " << log(pProb[iRefMap].forAngles[iProjectionOut])+pProb[iRefMap].ConstAngle[iProjectionOut]+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
+
+ }
+ }
+
+ }
+
+ angProbfile.close();
+ outputProbFile.close();
+
+ //Deleting allocated pointers
+
+ if (pProb)
+ {
+ delete[] pProb;
+ pProb = NULL;
+ }
+
+ if (param.refCTF)
+ {
+ delete[] param.refCTF;
+ param.refCTF =NULL;
+ }
+
+ return(0);
+}
+
+int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map)
+{
+#pragma omp parallel for
+ for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
+ {
+ compareRefMapShifted<-1>(iRefMap,iProjectionOut,iConv,conv_map, pProb, param.param_device, RefMap);
+ }
+ return(0);
+}
+
+int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
+{
+
+ /**************************************************************************************/
+ /**** BioEM Create Projection routine in Euler angle predefined grid****************
+ ********************* and turns projection into Fourier space **********************/
+ /**************************************************************************************/
+
+ myfloat3_t RotatedPointsModel[Model.nPointsModel];
+ myfloat_t rotmat[3][3];
+ myfloat_t alpha, gam,beta;
+ fftw_plan plan;
+ mycomplex_t* localproj;
+ int totnumPixFFT=2*param.param_device.NumberPixels;
+
+ localproj= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t) *4*param.param_device.NumberPixels*param.param_device.NumberPixels);
+ memset(localproj,0,4*param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localproj));
+
+ alpha=param.angles[iMap].pos[0];
+ beta=param.angles[iMap].pos[1];
+ gam=param.angles[iMap].pos[2];
+
+ /**** To see how things are going: cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***/
+
+ /********** Creat Rotation with pre-defiend grid of orientations**********/
+
+ rotmat[0][0]=cos(gam)*cos(alpha)-cos(beta)*sin(alpha)*sin(gam);
+ rotmat[0][1]=cos(gam)*sin(alpha)+cos(beta)*cos(alpha)*sin(gam);
+ rotmat[0][2]=sin(gam)*sin(beta);
+ rotmat[1][0]=-sin(gam)*cos(alpha)-cos(beta)*sin(alpha)*cos(gam);
+ rotmat[1][1]=-sin(gam)*sin(alpha)+cos(beta)*cos(alpha)*cos(gam);
+ rotmat[1][2]=cos(gam)*sin(beta);
+ rotmat[2][0]=sin(beta)*sin(alpha);
+ rotmat[2][1]=-sin(beta)*cos(alpha);
+ rotmat[2][2]=cos(beta);
+
+
+ for(int n=0; n< Model.nPointsModel; n++)
+ {
+ RotatedPointsModel[n].pos[0]=0.0;
+ RotatedPointsModel[n].pos[1]=0.0;
+ RotatedPointsModel[n].pos[2]=0.0;
+ }
+ for(int n=0; n< Model.nPointsModel; n++)
+ {
+ for(int k=0; k< 3; k++)
+ {
+ for(int j=0; j< 3; j++)
+ {
+ RotatedPointsModel[n].pos[k]+=rotmat[k][j]*Model.PointsModel[n].pos[j];
+ }
+ }
+ }
+
+ int i, j;
+
+ /************ Projection over the Z axis********************/
+ for(int n=0; n< Model.nPointsModel; n++)
+ {
+ //Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
+ i=floor(RotatedPointsModel[n].pos[0]/(*param.pixelSize)+float(param.param_device.NumberPixels)/2.0+0.5);
+ j=floor(RotatedPointsModel[n].pos[1]/(*param.pixelSize)+float(param.param_device.NumberPixels)/2.0+0.5);
+
+ localproj[i*2*param.param_device.NumberPixels+j+param.param_device.NumberPixels*param.param_device.NumberPixels+int(param.param_device.NumberPixels/2.0)][0]+=Model.densityPointsModel[n];
+
+ }
+
+ /**** Output Just to check****/
+ if(iMap==10)
+ {
+ ofstream myexamplemap;
+ ofstream myexampleRot;
+ myexamplemap.open ("MAP_i10");
+ myexampleRot.open ("Rot_i10");
+ myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
+ for(int k=0; k<2*param.param_device.NumberPixels; k++)
+ {
+ for(int j=0; j<2*param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*2*param.param_device.NumberPixels+j][0];
+ }
+ myexamplemap << " \n";
+ for(int n=0; n< Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
+ myexamplemap.close();
+ myexampleRot.close();
+ }
+
+ /***** Converting projection to Fourier Space for Convolution later with kernel****/
+ /********** Omp Critical is necessary with FFTW*******/
+ //#pragma omp critical
+ {
+ plan = fftw_plan_dft_2d(totnumPixFFT,totnumPixFFT,localproj,mapFFT,FFTW_FORWARD,FFTW_ESTIMATE);
+ fftw_execute(plan);
+ fftw_destroy_plan(plan);
+ fftw_free(localproj);
+ }
+
+ return(0);
+}
+
+
+int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv)
+{
+ /**************************************************************************************/
+ /**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
+ **************** calculated Projection with convoluted precalculated Kernel**********
+ *************** and Backtransforming it to real Space ******************************/
+ /**************************************************************************************/
+
+ fftw_plan plan;
+ mycomplex_t* localmultFFT;
+ mycomplex_t* localconvFFT;
+ int totnumPixFFT=2*param.param_device.NumberPixels;
+ localmultFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t)*totnumPixFFT*totnumPixFFT);
+ localconvFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t)*totnumPixFFT*totnumPixFFT);
+
+
+ /**** Multiplying FFTmap with corresponding kernel ****/
+
+ for(int i=0; i < 2*param.param_device.NumberPixels ; i++ )
+ {
+ for(int j=0; j < 2*param.param_device.NumberPixels ; j++ )
+ { //Projection*CONJ(KERNEL)
+ localmultFFT[i*2*param.param_device.NumberPixels+j][0]=lproj[i*2*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0]+lproj[i*2*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1];
+ localmultFFT[i*2*param.param_device.NumberPixels+j][1]=lproj[i*2*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0]-lproj[i*2*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1];
+ // cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0] << " " <<param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1] <<" " <<lproj[i*2*param.param_device.NumberPixels+j][0] <<" " <<lproj[i*2*param.param_device.NumberPixels+j][1] << "\n";
+ }
+ }
+
+ /**** Bringing convoluted Map to real Space ****/
+ //#pragma omp critical
+ {
+ plan = fftw_plan_dft_2d(totnumPixFFT,totnumPixFFT,localmultFFT,localconvFFT,FFTW_BACKWARD,FFTW_ESTIMATE);
+ fftw_execute(plan);
+ }
+
+
+ /****Asigning convolution fftw_complex to bioem_map ****/
+ for(int i=0; i < param.param_device.NumberPixels ; i++ )
+ {
+ for(int j=0; j < param.param_device.NumberPixels ; j++ )
+ {
+ Mapconv.points[i][j]=localconvFFT[i*2*param.param_device.NumberPixels+j+param.param_device.NumberPixels*param.param_device.NumberPixels+int(param.param_device.NumberPixels/2.0)][0];
+ }
+ }
+
+ /**** Freeing fftw_complex created (dont know if omp critical is necessary) ****/
+ //#pragma omp critical
+ {
+ fftw_destroy_plan(plan);
+ fftw_free(localconvFFT);
+ fftw_free(localmultFFT);
+ }
+
+ return(0);
+}
+
+int bioem::calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare)
+{
+ /*********************** Routine to calculate Cross correlations***********************/
+
+ sum=0.0;
+ sumsquare=0.0;
+ for (int i = 0; i < param.param_device.NumberPixels; i++)
+ {
+ for (int j = 0; j < param.param_device.NumberPixels; j++)
+ {
+ // Calculate Sum of pixels
+ sum += localmap.points[i][j];
+ // Calculate Sum of pixels squared
+ sumsquare += localmap.points[i][j]*localmap.points[i][j];
+ }
+ }
+ return(0);
+}
+
+int bioem::deviceInit()
+{
+ return(0);
+}
+
+int bioem::deviceStartRun()
+{
+ return(0);
+}
+
+int bioem::deviceFinishRun()
+{
+ return(0);
+}
diff --git a/bioem_cuda.cu b/bioem_cuda.cu
index c1c4230..2af9817 100644
--- a/bioem_cuda.cu
+++ b/bioem_cuda.cu
@@ -1,210 +1,210 @@
-#define BIOEM_GPUCODE
-
-#if defined(_WIN32)
-#include <windows.h>
-#endif
-
-#include <iostream>
-using namespace std;
-
-#include "bioem_cuda_internal.h"
-#include "bioem_algorithm.h"
-
-//__constant__ bioem_map gConvMap;
-
-static inline void checkCudaErrors(cudaError_t error)
-{
- if (error != cudaSuccess)
- {
- printf("CUDA Error %d / %s\n", error, cudaGetErrorString(error));
- exit(1);
- }
-}
-
-bioem_cuda::bioem_cuda()
-{
- deviceInitialized = 0;
- GPUAlgo = getenv("GPUALGO") == NULL ? 2 : atoi(getenv("GPUALGO"));
- GPUAsync = getenv("GPUASYNC") == NULL ? 0 : atoi(getenv("GPUASYNC"));
-}
-
-bioem_cuda::~bioem_cuda()
-{
- deviceExit();
-}
-
-__global__ void compareRefMap_kernel(const int iOrient, const int iConv, const bioem_map* pMap, bioem_Probability* pProb, const bioem_param_device param, const bioem_RefMap_Mod* RefMap, const int cent_x, const int cent_y)
-{
- const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
- if (iRefMap < RefMap->ntotRefMap)
- {
- compareRefMap<0>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap, cent_x, cent_y);
- }
-}
-
-__global__ void compareRefMapShifted_kernel(const int iOrient, const int iConv, const bioem_map* pMap, bioem_Probability* pProb, const bioem_param_device param, const bioem_RefMap_Mod* RefMap)
-{
- const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
- if (iRefMap < RefMap->ntotRefMap)
- {
- compareRefMapShifted<1>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap);
- }
-}
-
-__global__ void compareRefMapLoopShifts_kernel(const int iOrient, const int iConv, const bioem_map* pMap, bioem_Probability* pProb, const bioem_param_device param, const bioem_RefMap* RefMap, const int blockoffset, const int nShifts, const int nShiftBits)
-{
- const size_t myid = (myBlockIdxX + blockoffset) * myBlockDimX + myThreadIdxX;
- const int iRefMap = myid >> (nShiftBits << 1);
- const int myRef = myThreadIdxX >> (nShiftBits << 1);
- const int myShiftIdx = (myid >> nShiftBits) & (nShifts - 1);
- const int myShiftIdy = myid & (nShifts - 1);
- const int myShift = myid & (nShifts * nShifts - 1);
- const int cent_x = myShiftIdx * param.GridSpaceCenter - param.maxDisplaceCenter;
- const int cent_y = myShiftIdy * param.GridSpaceCenter - param.maxDisplaceCenter;
-
- compareRefMap<2>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap, cent_x, cent_y, myShift, nShifts * nShifts, myRef);
-}
-
-template <class T> static inline T divup(T num, T divider) {return((num + divider - 1) / divider);}
-static inline bool IsPowerOf2(int x) {return ((x > 0) && ((x & (x - 1)) == 0));}
-#if defined(_WIN32)
-static inline int ilog2 (int value)
-{
- DWORD index;
- _BitScanReverse (&index, value);
- return(value);
-}
-#else
-static inline int ilog2(int value) {return 31 - __builtin_clz(value);}
-#endif
-
-int bioem_cuda::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map)
-{
- //cudaMemcpyToSymbolAsync(&gConvMap, &conv_map, sizeof(bioem_map), 0, cudaMemcpyHostToDevice, cudaStream);
- if (GPUAsync)
- {
- checkCudaErrors(cudaEventSynchronize(cudaEvent[iConv & 1]));
- }
- checkCudaErrors(cudaMemcpyAsync(pConvMap_device[iConv & 1], &conv_map, sizeof(bioem_map), cudaMemcpyHostToDevice, cudaStream));
-
- if (GPUAlgo == 2) //Loop over shifts
- {
- const int nShifts = 2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1;
- if (!IsPowerOf2(nShifts))
- {
- cout << "Invalid number of displacements, no power of two\n";
- exit(1);
- }
- if (CUDA_THREAD_COUNT % (nShifts * nShifts))
- {
- cout << "CUDA Thread count (" << CUDA_THREAD_COUNT << ") is no multiple of number of shifts (" << (nShifts * nShifts) << ")\n";
- exit(1);
- }
- if (nShifts > CUDA_MAX_SHIFT_REDUCE)
- {
- cout << "Too many displacements for CUDA reduction\n";
- exit(1);
- }
- const int nShiftBits = ilog2(nShifts);
- size_t totalBlocks = divup((size_t) RefMap.ntotRefMap * (size_t) nShifts * (size_t) nShifts, (size_t) CUDA_THREAD_COUNT);
- size_t nBlocks = CUDA_BLOCK_COUNT;
- for (size_t i = 0;i < totalBlocks;i += nBlocks)
- {
- compareRefMapLoopShifts_kernel <<<min(nBlocks, totalBlocks - i), CUDA_THREAD_COUNT, (CUDA_THREAD_COUNT * 2 + CUDA_THREAD_COUNT / (nShifts * nShifts) * 4) * sizeof(myfloat_t), cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device, i, nShifts, nShiftBits);
- }
- }
- else if (GPUAlgo == 1) //Split shifts in multiple kernels
- {
- for (int cent_x = -param.param_device.maxDisplaceCenter; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter)
- {
- for (int cent_y = -param.param_device.maxDisplaceCenter; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
- {
- compareRefMap_kernel <<<divup(RefMap.ntotRefMap, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, cent_x, cent_y);
- }
- }
- }
- else if (GPUAlgo == 0) //All shifts in one kernel
- {
- compareRefMapShifted_kernel <<<divup(RefMap.ntotRefMap, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod);
- }
- else
- {
- cout << "Invalid GPU Algorithm selected\n";
- exit(1);
- }
- if (GPUAsync)
- {
- checkCudaErrors(cudaEventRecord(cudaEvent[iConv & 1], cudaStream));
- }
- else
- {
- checkCudaErrors(cudaStreamSynchronize(cudaStream));
- }
- return(0);
-}
-
-int bioem_cuda::deviceInit()
-{
- deviceExit();
-
- checkCudaErrors(cudaStreamCreate(&cudaStream));
- checkCudaErrors(cudaMalloc(&pRefMap_device, sizeof(bioem_RefMap)));
- checkCudaErrors(cudaMalloc(&pProb_device, sizeof(bioem_Probability) * RefMap.ntotRefMap));
- for (int i = 0;i < 2;i++)
- {
- checkCudaErrors(cudaEventCreate(&cudaEvent[i]));
- checkCudaErrors(cudaMalloc(&pConvMap_device[i], sizeof(bioem_map)));
- }
- pRefMap_device_Mod = (bioem_RefMap_Mod*) pRefMap_device;
-
- if (GPUAlgo == 0 || GPUAlgo == 1)
- {
- bioem_RefMap_Mod* RefMapGPU = new bioem_RefMap_Mod(RefMap);
- cudaMemcpy(pRefMap_device_Mod, RefMapGPU, sizeof(bioem_RefMap), cudaMemcpyHostToDevice);
- delete RefMapGPU;
- }
- else
- {
- cudaMemcpy(pRefMap_device, &RefMap, sizeof(bioem_RefMap), cudaMemcpyHostToDevice);
- }
- deviceInitialized = 1;
- return(0);
-}
-
-int bioem_cuda::deviceExit()
-{
- if (deviceInitialized == 0) return(0);
-
- cudaStreamDestroy(cudaStream);
- cudaFree(pRefMap_device);
- cudaFree(pProb_device);
- for (int i = 0;i < 2;i++)
- {
- cudaEventDestroy(cudaEvent[i]);
- cudaFree(pConvMap_device);
- }
- cudaThreadExit();
-
- deviceInitialized = 0;
- return(0);
-}
-
-int bioem_cuda::deviceStartRun()
-{
-
- cudaMemcpy(pProb_device, pProb, sizeof(bioem_Probability) * RefMap.ntotRefMap, cudaMemcpyHostToDevice);
- return(0);
-}
-
-int bioem_cuda::deviceFinishRun()
-{
- if (GPUAsync) cudaStreamSynchronize(cudaStream);
- cudaMemcpy(pProb, pProb_device, sizeof(bioem_Probability) * RefMap.ntotRefMap, cudaMemcpyDeviceToHost);
-
- return(0);
-}
-
-bioem* bioem_cuda_create()
-{
- return new bioem_cuda;
-}
+#define BIOEM_GPUCODE
+
+#if defined(_WIN32)
+#include <windows.h>
+#endif
+
+#include <iostream>
+using namespace std;
+
+#include "bioem_cuda_internal.h"
+#include "bioem_algorithm.h"
+
+//__constant__ bioem_map gConvMap;
+
+static inline void checkCudaErrors(cudaError_t error)
+{
+ if (error != cudaSuccess)
+ {
+ printf("CUDA Error %d / %s\n", error, cudaGetErrorString(error));
+ exit(1);
+ }
+}
+
+bioem_cuda::bioem_cuda()
+{
+ deviceInitialized = 0;
+ GPUAlgo = getenv("GPUALGO") == NULL ? 2 : atoi(getenv("GPUALGO"));
+ GPUAsync = getenv("GPUASYNC") == NULL ? 0 : atoi(getenv("GPUASYNC"));
+}
+
+bioem_cuda::~bioem_cuda()
+{
+ deviceExit();
+}
+
+__global__ void compareRefMap_kernel(const int iOrient, const int iConv, const bioem_map* pMap, bioem_Probability* pProb, const bioem_param_device param, const bioem_RefMap_Mod* RefMap, const int cent_x, const int cent_y)
+{
+ const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
+ if (iRefMap < RefMap->ntotRefMap)
+ {
+ compareRefMap<0>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap, cent_x, cent_y);
+ }
+}
+
+__global__ void compareRefMapShifted_kernel(const int iOrient, const int iConv, const bioem_map* pMap, bioem_Probability* pProb, const bioem_param_device param, const bioem_RefMap_Mod* RefMap)
+{
+ const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
+ if (iRefMap < RefMap->ntotRefMap)
+ {
+ compareRefMapShifted<1>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap);
+ }
+}
+
+__global__ void compareRefMapLoopShifts_kernel(const int iOrient, const int iConv, const bioem_map* pMap, bioem_Probability* pProb, const bioem_param_device param, const bioem_RefMap* RefMap, const int blockoffset, const int nShifts, const int nShiftBits)
+{
+ const size_t myid = (myBlockIdxX + blockoffset) * myBlockDimX + myThreadIdxX;
+ const int iRefMap = myid >> (nShiftBits << 1);
+ const int myRef = myThreadIdxX >> (nShiftBits << 1);
+ const int myShiftIdx = (myid >> nShiftBits) & (nShifts - 1);
+ const int myShiftIdy = myid & (nShifts - 1);
+ const int myShift = myid & (nShifts * nShifts - 1);
+ const int cent_x = myShiftIdx * param.GridSpaceCenter - param.maxDisplaceCenter;
+ const int cent_y = myShiftIdy * param.GridSpaceCenter - param.maxDisplaceCenter;
+
+ compareRefMap<2>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap, cent_x, cent_y, myShift, nShifts * nShifts, myRef);
+}
+
+template <class T> static inline T divup(T num, T divider) {return((num + divider - 1) / divider);}
+static inline bool IsPowerOf2(int x) {return ((x > 0) && ((x & (x - 1)) == 0));}
+#if defined(_WIN32)
+static inline int ilog2 (int value)
+{
+ DWORD index;
+ _BitScanReverse (&index, value);
+ return(value);
+}
+#else
+static inline int ilog2(int value) {return 31 - __builtin_clz(value);}
+#endif
+
+int bioem_cuda::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map)
+{
+ //cudaMemcpyToSymbolAsync(&gConvMap, &conv_map, sizeof(bioem_map), 0, cudaMemcpyHostToDevice, cudaStream);
+ if (GPUAsync)
+ {
+ checkCudaErrors(cudaEventSynchronize(cudaEvent[iConv & 1]));
+ }
+ checkCudaErrors(cudaMemcpyAsync(pConvMap_device[iConv & 1], &conv_map, sizeof(bioem_map), cudaMemcpyHostToDevice, cudaStream));
+
+ if (GPUAlgo == 2) //Loop over shifts
+ {
+ const int nShifts = 2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1;
+ if (!IsPowerOf2(nShifts))
+ {
+ cout << "Invalid number of displacements, no power of two\n";
+ exit(1);
+ }
+ if (CUDA_THREAD_COUNT % (nShifts * nShifts))
+ {
+ cout << "CUDA Thread count (" << CUDA_THREAD_COUNT << ") is no multiple of number of shifts (" << (nShifts * nShifts) << ")\n";
+ exit(1);
+ }
+ if (nShifts > CUDA_MAX_SHIFT_REDUCE)
+ {
+ cout << "Too many displacements for CUDA reduction\n";
+ exit(1);
+ }
+ const int nShiftBits = ilog2(nShifts);
+ size_t totalBlocks = divup((size_t) RefMap.ntotRefMap * (size_t) nShifts * (size_t) nShifts, (size_t) CUDA_THREAD_COUNT);
+ size_t nBlocks = CUDA_BLOCK_COUNT;
+ for (size_t i = 0;i < totalBlocks;i += nBlocks)
+ {
+ compareRefMapLoopShifts_kernel <<<min(nBlocks, totalBlocks - i), CUDA_THREAD_COUNT, (CUDA_THREAD_COUNT * 2 + CUDA_THREAD_COUNT / (nShifts * nShifts) * 4) * sizeof(myfloat_t), cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device, i, nShifts, nShiftBits);
+ }
+ }
+ else if (GPUAlgo == 1) //Split shifts in multiple kernels
+ {
+ for (int cent_x = -param.param_device.maxDisplaceCenter; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter)
+ {
+ for (int cent_y = -param.param_device.maxDisplaceCenter; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
+ {
+ compareRefMap_kernel <<<divup(RefMap.ntotRefMap, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, cent_x, cent_y);
+ }
+ }
+ }
+ else if (GPUAlgo == 0) //All shifts in one kernel
+ {
+ compareRefMapShifted_kernel <<<divup(RefMap.ntotRefMap, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod);
+ }
+ else
+ {
+ cout << "Invalid GPU Algorithm selected\n";
+ exit(1);
+ }
+ if (GPUAsync)
+ {
+ checkCudaErrors(cudaEventRecord(cudaEvent[iConv & 1], cudaStream));
+ }
+ else
+ {
+ checkCudaErrors(cudaStreamSynchronize(cudaStream));
+ }
+ return(0);
+}
+
+int bioem_cuda::deviceInit()
+{
+ deviceExit();
+
+ checkCudaErrors(cudaStreamCreate(&cudaStream));
+ checkCudaErrors(cudaMalloc(&pRefMap_device, sizeof(bioem_RefMap)));
+ checkCudaErrors(cudaMalloc(&pProb_device, sizeof(bioem_Probability) * RefMap.ntotRefMap));
+ for (int i = 0;i < 2;i++)
+ {
+ checkCudaErrors(cudaEventCreate(&cudaEvent[i]));
+ checkCudaErrors(cudaMalloc(&pConvMap_device[i], sizeof(bioem_map)));
+ }
+ pRefMap_device_Mod = (bioem_RefMap_Mod*) pRefMap_device;
+
+ if (GPUAlgo == 0 || GPUAlgo == 1)
+ {
+ bioem_RefMap_Mod* RefMapGPU = new bioem_RefMap_Mod(RefMap);
+ cudaMemcpy(pRefMap_device_Mod, RefMapGPU, sizeof(bioem_RefMap), cudaMemcpyHostToDevice);
+ delete RefMapGPU;
+ }
+ else
+ {
+ cudaMemcpy(pRefMap_device, &RefMap, sizeof(bioem_RefMap), cudaMemcpyHostToDevice);
+ }
+ deviceInitialized = 1;
+ return(0);
+}
+
+int bioem_cuda::deviceExit()
+{
+ if (deviceInitialized == 0) return(0);
+
+ cudaStreamDestroy(cudaStream);
+ cudaFree(pRefMap_device);
+ cudaFree(pProb_device);
+ for (int i = 0;i < 2;i++)
+ {
+ cudaEventDestroy(cudaEvent[i]);
+ cudaFree(pConvMap_device);
+ }
+ cudaThreadExit();
+
+ deviceInitialized = 0;
+ return(0);
+}
+
+int bioem_cuda::deviceStartRun()
+{
+
+ cudaMemcpy(pProb_device, pProb, sizeof(bioem_Probability) * RefMap.ntotRefMap, cudaMemcpyHostToDevice);
+ return(0);
+}
+
+int bioem_cuda::deviceFinishRun()
+{
+ if (GPUAsync) cudaStreamSynchronize(cudaStream);
+ cudaMemcpy(pProb, pProb_device, sizeof(bioem_Probability) * RefMap.ntotRefMap, cudaMemcpyDeviceToHost);
+
+ return(0);
+}
+
+bioem* bioem_cuda_create()
+{
+ return new bioem_cuda;
+}
diff --git a/include/bioem.h b/include/bioem.h
index b61cffc..afed15d 100644
--- a/include/bioem.h
+++ b/include/bioem.h
@@ -1,48 +1,48 @@
-#ifndef BIOEM_H
-#define BIOEM_H
-#include <fstream>
-#include <iostream>
-#include <complex>
-
-#include "defs.h"
-#include "bioem.h"
-#include "model.h"
-#include "map.h"
-#include "param.h"
-
-class bioem
-{
-public:
- bioem();
- virtual ~bioem();
-
- int configure(int ac, char* av[]);
- int precalculate(); // Is it better to pass directly the input File names?
- int dopreCalCrossCorrelation(int iRefMap, int iRefMapLocal);
- int run();
- int doProjections(int iMap);
- int createConvolutedProjectionMap(int iOreint,int iMap, mycomplex_t* lproj,bioem_map& Mapconv);
-
- virtual int compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map);
- int createProjection(int iMap, mycomplex_t* map);
- int calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare);
-
- bioem_Probability* pProb;
-
-protected:
- virtual int deviceInit();
- virtual int deviceStartRun();
- virtual int deviceFinishRun();
-
- bioem_param param;
- bioem_model Model;
- bioem_RefMap RefMap;
-
- int nReferenceMaps; //Maps in memory at a time
- int nReferenceMapsTotal; //Maps in total
-
- int nProjectionMaps; //Maps in memory at a time
- int nProjectionMapsTotal; //Maps in total
-};
-
-#endif
+#ifndef BIOEM_H
+#define BIOEM_H
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+#include "defs.h"
+#include "bioem.h"
+#include "model.h"
+#include "map.h"
+#include "param.h"
+
+class bioem
+{
+public:
+ bioem();
+ virtual ~bioem();
+
+ int configure(int ac, char* av[]);
+ int precalculate(); // Is it better to pass directly the input File names?
+ int dopreCalCrossCorrelation(int iRefMap, int iRefMapLocal);
+ int run();
+ int doProjections(int iMap);
+ int createConvolutedProjectionMap(int iOreint,int iMap, mycomplex_t* lproj,bioem_map& Mapconv);
+
+ virtual int compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map);
+ int createProjection(int iMap, mycomplex_t* map);
+ int calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare);
+
+ bioem_Probability* pProb;
+
+protected:
+ virtual int deviceInit();
+ virtual int deviceStartRun();
+ virtual int deviceFinishRun();
+
+ bioem_param param;
+ bioem_model Model;
+ bioem_RefMap RefMap;
+
+ int nReferenceMaps; //Maps in memory at a time
+ int nReferenceMapsTotal; //Maps in total
+
+ int nProjectionMaps; //Maps in memory at a time
+ int nProjectionMapsTotal; //Maps in total
+};
+
+#endif
diff --git a/include/defs.h b/include/defs.h
index 6381287..7d8c82a 100644
--- a/include/defs.h
+++ b/include/defs.h
@@ -1,42 +1,42 @@
-#ifndef BIOEM_DEFS_H
-#define BIOEM_DEFS_H
-
-typedef float myfloat_t;
-typedef double mycomplex_t[2];
-
-#define BIOEM_FLOAT_3_PHYSICAL_SIZE 3 //Possible set to 4 for GPU
-#define BIOEM_MAP_SIZE_X 224
-#define BIOEM_MAP_SIZE_Y 224
-#define BIOEM_MODEL_SIZE 120000
-#define BIOEM_MAX_MAPS 12000
-#define MAX_REF_CTF 200
-#define MAX_ORIENT 20000
-
-struct myfloat3_t
-{
- myfloat_t pos[BIOEM_FLOAT_3_PHYSICAL_SIZE];
-};
-
-#ifdef BIOEM_GPUCODE
-#define myThreadIdxX threadIdx.x
-#define myThreadIdxY threadIdx.y
-#define myBlockDimX blockDim.x
-#define myBlockDimY blockDim.y
-#define myBlockIdxX blockIdx.x
-#define myBlockIdxY blockIdx.y
-#else
-#define __device__
-#define __host__
-#define myThreadIdxX 0
-#define myThreadIdxY 0
-#define myBlockDimX 1
-#define myBlockDimY 1
-#define myBlockIdxX 0
-#define myBlockIdxY 0
-#endif
-
-#define CUDA_THREAD_COUNT 256
-#define CUDA_BLOCK_COUNT 1024 * 16
-#define CUDA_MAX_SHIFT_REDUCE 1024
-
-#endif
+#ifndef BIOEM_DEFS_H
+#define BIOEM_DEFS_H
+
+typedef float myfloat_t;
+typedef double mycomplex_t[2];
+
+#define BIOEM_FLOAT_3_PHYSICAL_SIZE 3 //Possible set to 4 for GPU
+#define BIOEM_MAP_SIZE_X 224
+#define BIOEM_MAP_SIZE_Y 224
+#define BIOEM_MODEL_SIZE 120000
+#define BIOEM_MAX_MAPS 12000
+#define MAX_REF_CTF 200
+#define MAX_ORIENT 20000
+
+struct myfloat3_t
+{
+ myfloat_t pos[BIOEM_FLOAT_3_PHYSICAL_SIZE];
+};
+
+#ifdef BIOEM_GPUCODE
+#define myThreadIdxX threadIdx.x
+#define myThreadIdxY threadIdx.y
+#define myBlockDimX blockDim.x
+#define myBlockDimY blockDim.y
+#define myBlockIdxX blockIdx.x
+#define myBlockIdxY blockIdx.y
+#else
+#define __device__
+#define __host__
+#define myThreadIdxX 0
+#define myThreadIdxY 0
+#define myBlockDimX 1
+#define myBlockDimY 1
+#define myBlockIdxX 0
+#define myBlockIdxY 0
+#endif
+
+#define CUDA_THREAD_COUNT 256
+#define CUDA_BLOCK_COUNT 1024 * 16
+#define CUDA_MAX_SHIFT_REDUCE 1024
+
+#endif
diff --git a/include/map.h b/include/map.h
index 4c8e5ad..c053851 100644
--- a/include/map.h
+++ b/include/map.h
@@ -1,88 +1,88 @@
-#ifndef BIOEM_MAP_H
-#define BIOEM_MAP_H
-
-#include "defs.h"
-#include <complex>
-#include <math.h>
-
-class bioem_param;
-
-class bioem_map
-{
-public:
- myfloat_t points[BIOEM_MAP_SIZE_X][BIOEM_MAP_SIZE_Y];
-};
-
-class bioem_map_forFFT
-{
-public:
- mycomplex_t cpoints[2*BIOEM_MAP_SIZE_X*2*BIOEM_MAP_SIZE_Y];
-
-};
-
-class bioem_RefMap
-{
-public:
- int readRefMaps();
-
- const char* filemap;
- int ntotRefMap;
- bioem_map Ref[BIOEM_MAX_MAPS];
- myfloat_t sum_RefMap[BIOEM_MAX_MAPS];
- myfloat_t sumsquare_RefMap[BIOEM_MAX_MAPS];
- myfloat_t ForLogProbfromRef[BIOEM_MAX_MAPS];
-
- bool dumpMap, loadMap;
-
- __host__ __device__ inline myfloat_t get(int map, int x, int y) const {return(Ref[map].points[x][y]);}
- __host__ __device__ inline const myfloat_t* getp(int map, int x, int y) const {return(&Ref[map].points[x][y]);}
-};
-
-
-class bioem_RefMap_Mod
-{
-public:
-
- const char* filemap;
- int ntotRefMap;
- myfloat_t Ref[BIOEM_MAP_SIZE_X][BIOEM_MAP_SIZE_Y][BIOEM_MAX_MAPS];
- myfloat_t sum_RefMap[BIOEM_MAX_MAPS];
- myfloat_t sumsquare_RefMap[BIOEM_MAX_MAPS];
- myfloat_t ForLogProbfromRef[BIOEM_MAX_MAPS];
-
- __host__ __device__ inline myfloat_t get(int map, int x, int y) const {return(Ref[x][y][map]);}
-
- bioem_RefMap_Mod() {ntotRefMap = 0;}
-
- bioem_RefMap_Mod(const bioem_RefMap& map)
- {
- ntotRefMap = map.ntotRefMap;
- memcpy(sum_RefMap, map.sum_RefMap, sizeof(sum_RefMap));
- memcpy(sumsquare_RefMap, map.sumsquare_RefMap, sizeof(sumsquare_RefMap));
- memcpy(ForLogProbfromRef, map.ForLogProbfromRef, sizeof(ForLogProbfromRef));
-#pragma omp parallel for
- for (int i = 0;i < ntotRefMap;i++)
- {
- for (int j = 0;j < BIOEM_MAP_SIZE_X;j++)
- {
- for (int k = 0;k < BIOEM_MAP_SIZE_Y;k++)
- {
- Ref[j][k][i] = map.get(i, j, k);
- }
- }
- }
- }
-};
-
-
-class bioem_Probability
-{
-public:
- myfloat_t forAngles[MAX_ORIENT];
- myfloat_t Total;
- myfloat_t Constoadd;
- myfloat_t ConstAngle[MAX_ORIENT];
- myfloat_t max_prob;
- int max_prob_cent_x, max_prob_cent_y, max_prob_orient, max_prob_conv;
-};
-#endif
+#ifndef BIOEM_MAP_H
+#define BIOEM_MAP_H
+
+#include "defs.h"
+#include <complex>
+#include <math.h>
+
+class bioem_param;
+
+class bioem_map
+{
+public:
+ myfloat_t points[BIOEM_MAP_SIZE_X][BIOEM_MAP_SIZE_Y];
+};
+
+class bioem_map_forFFT
+{
+public:
+ mycomplex_t cpoints[2*BIOEM_MAP_SIZE_X*2*BIOEM_MAP_SIZE_Y];
+
+};
+
+class bioem_RefMap
+{
+public:
+ int readRefMaps();
+
+ const char* filemap;
+ int ntotRefMap;
+ bioem_map Ref[BIOEM_MAX_MAPS];
+ myfloat_t sum_RefMap[BIOEM_MAX_MAPS];
+ myfloat_t sumsquare_RefMap[BIOEM_MAX_MAPS];
+ myfloat_t ForLogProbfromRef[BIOEM_MAX_MAPS];
+
+ bool dumpMap, loadMap;
+
+ __host__ __device__ inline myfloat_t get(int map, int x, int y) const {return(Ref[map].points[x][y]);}
+ __host__ __device__ inline const myfloat_t* getp(int map, int x, int y) const {return(&Ref[map].points[x][y]);}
+};
+
+
+class bioem_RefMap_Mod
+{
+public:
+
+ const char* filemap;
+ int ntotRefMap;
+ myfloat_t Ref[BIOEM_MAP_SIZE_X][BIOEM_MAP_SIZE_Y][BIOEM_MAX_MAPS];
+ myfloat_t sum_RefMap[BIOEM_MAX_MAPS];
+ myfloat_t sumsquare_RefMap[BIOEM_MAX_MAPS];
+ myfloat_t ForLogProbfromRef[BIOEM_MAX_MAPS];
+
+ __host__ __device__ inline myfloat_t get(int map, int x, int y) const {return(Ref[x][y][map]);}
+
+ bioem_RefMap_Mod() {ntotRefMap = 0;}
+
+ bioem_RefMap_Mod(const bioem_RefMap& map)
+ {
+ ntotRefMap = map.ntotRefMap;
+ memcpy(sum_RefMap, map.sum_RefMap, sizeof(sum_RefMap));
+ memcpy(sumsquare_RefMap, map.sumsquare_RefMap, sizeof(sumsquare_RefMap));
+ memcpy(ForLogProbfromRef, map.ForLogProbfromRef, sizeof(ForLogProbfromRef));
+#pragma omp parallel for
+ for (int i = 0;i < ntotRefMap;i++)
+ {
+ for (int j = 0;j < BIOEM_MAP_SIZE_X;j++)
+ {
+ for (int k = 0;k < BIOEM_MAP_SIZE_Y;k++)
+ {
+ Ref[j][k][i] = map.get(i, j, k);
+ }
+ }
+ }
+ }
+};
+
+
+class bioem_Probability
+{
+public:
+ myfloat_t forAngles[MAX_ORIENT];
+ myfloat_t Total;
+ myfloat_t Constoadd;
+ myfloat_t ConstAngle[MAX_ORIENT];
+ myfloat_t max_prob;
+ int max_prob_cent_x, max_prob_cent_y, max_prob_orient, max_prob_conv;
+};
+#endif
diff --git a/include/model.h b/include/model.h
index 55b6a2e..9e94708 100644
--- a/include/model.h
+++ b/include/model.h
@@ -1,27 +1,27 @@
-#ifndef BIOEM_MODEL_H
-#define BIOEM_MODEL_H
-
-#include "defs.h"
-
-class bioem_model
-{
-public:
- //bioem_model();
- //~bioem_model();
-
- int readModel();
-
- bool readPDB;
-
- myfloat_t getAminoAcidRad(char *name);
- myfloat_t getAminoAcidDensity(char *name);
-
- const char* filemodel;
-
- int nPointsModel;
- myfloat3_t PointsModel[BIOEM_MODEL_SIZE];
- myfloat_t radiusPointsModel[BIOEM_MODEL_SIZE];
- myfloat_t densityPointsModel[BIOEM_MODEL_SIZE];
-};
-
-#endif
+#ifndef BIOEM_MODEL_H
+#define BIOEM_MODEL_H
+
+#include "defs.h"
+
+class bioem_model
+{
+public:
+ //bioem_model();
+ //~bioem_model();
+
+ int readModel();
+
+ bool readPDB;
+
+ myfloat_t getAminoAcidRad(char *name);
+ myfloat_t getAminoAcidDensity(char *name);
+
+ const char* filemodel;
+
+ int nPointsModel;
+ myfloat3_t PointsModel[BIOEM_MODEL_SIZE];
+ myfloat_t radiusPointsModel[BIOEM_MODEL_SIZE];
+ myfloat_t densityPointsModel[BIOEM_MODEL_SIZE];
+};
+
+#endif
diff --git a/main.cpp b/main.cpp
index 744c50e..4454603 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,59 +1,59 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <omp.h>
-#include <time.h>
-#include <fenv.h>
-
-#ifdef _WIN32
-#include <Windows.h>
-#include <WinBase.h>
-#endif
-
-#include <iostream>
-#include <algorithm>
-#include <iterator>
-#include "bioem.h"
-#include "bioem_cuda.h"
-#include <omp.h>
-
-#include "cmodules/timer.h"
-
-int main(int argc, char* argv[])
-{
- /**************************************************************************************/
- /********************************* Main BioEM code **********************************/
- /************************************************************************************/
-
-#pragma omp parallel
- {
- _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); //Flush denormals to zero in all OpenMP threads
- }
- HighResTimer timer;
-
- bioem* bio;
- if (getenv("GPU") && atoi(getenv("GPU")))
- {
- bio = bioem_cuda_create();
- }
- else
- {
- bio = new bioem;
- }
-
- /************ Configuration and Pre-calculating necessary objects *****************/
- printf("Configuring\n");
- bio->configure(argc,argv);
-
- /******************************* Run BioEM routine ******************************/
- printf("Running\n");
- timer.Start();
- bio->run();
- timer.Stop();
-
- /************************************ End **********************************/
- printf ("The code ran for %f seconds.\n", timer.GetElapsedTime());
- delete bio;
-
- return(0);
-}
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <omp.h>
+#include <time.h>
+#include <fenv.h>
+
+#ifdef _WIN32
+#include <Windows.h>
+#include <WinBase.h>
+#endif
+
+#include <iostream>
+#include <algorithm>
+#include <iterator>
+#include "bioem.h"
+#include "bioem_cuda.h"
+#include <omp.h>
+
+#include "cmodules/timer.h"
+
+int main(int argc, char* argv[])
+{
+ /**************************************************************************************/
+ /********************************* Main BioEM code **********************************/
+ /************************************************************************************/
+
+#pragma omp parallel
+ {
+ _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); //Flush denormals to zero in all OpenMP threads
+ }
+ HighResTimer timer;
+
+ bioem* bio;
+ if (getenv("GPU") && atoi(getenv("GPU")))
+ {
+ bio = bioem_cuda_create();
+ }
+ else
+ {
+ bio = new bioem;
+ }
+
+ /************ Configuration and Pre-calculating necessary objects *****************/
+ printf("Configuring\n");
+ bio->configure(argc,argv);
+
+ /******************************* Run BioEM routine ******************************/
+ printf("Running\n");
+ timer.Start();
+ bio->run();
+ timer.Stop();
+
+ /************************************ End **********************************/
+ printf ("The code ran for %f seconds.\n", timer.GetElapsedTime());
+ delete bio;
+
+ return(0);
+}
diff --git a/map.cpp b/map.cpp
index 4dd479d..c2136f5 100644
--- a/map.cpp
+++ b/map.cpp
@@ -1,116 +1,116 @@
-#include <fstream>
-#include <iostream>
-#include <stdio.h>
-#include <stdlib.h>
-#include <cstring>
-
-#include "map.h"
-
-using namespace std;
-
-int bioem_RefMap::readRefMaps()
-{
- /**************************************************************************************/
- /***********************Reading reference Particle Maps************************/
- /**************************************************************************************/
- if (loadMap)
- {
- FILE* fp = fopen("maps.dump", "rb");
- if (fp == NULL)
- {
- cout << "Error opening dump file\n";
- exit(1);
- }
- fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
- if (ntotRefMap > BIOEM_MAX_MAPS)
- {
- cout << "BIOEM_MAX_MAPS too small\n";
- exit(1);
- }
- fread(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
- fclose(fp);
-
- cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
- cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
- }
- else
- {
- int nummap=-1;
- ifstream input(filemap);
- if (!input.good())
- {
- cout << "Particle Maps Failed to open file" << endl ;
- exit(1);
- }
-
- char line[512] = {' '};
- char tmpLine[512] = {' '};
-
- while (!input.eof())
- {
- input.getline(line,512);
-
- strncpy(tmpLine,line,strlen(line));
- char *token = strtok(tmpLine," ");
-
- if (strcmp(token,"PARTICLE")==0) // to count the number of maps
- {
- nummap++;
- if (nummap % 128 == 0)
- {
- cout << "..." << nummap << "\n";
- }
- if (nummap == BIOEM_MAX_MAPS)
- {
- cout << "BIOEM_MAX_MAPS too small\n";
- exit(1);
- }
- }
- else
- {
- int i,j;
- myfloat_t z;
-
- char tmpVals[36] = {' '};
-
- strncpy (tmpVals,line,8);
- sscanf (tmpVals,"%d",&i);
-
- strncpy (tmpVals,line+8,8);
- sscanf (tmpVals,"%d",&j);
-
- strncpy (tmpVals,line+16,12);
- sscanf (tmpVals,"%f",&z);
- //checking for Map limits
- if(i>0 && i-1 <BIOEM_MAP_SIZE_X&& j>0 && j-1<BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
- {
- Ref[nummap].points[i-1][j-1]=z;
- }
- else
- {
- cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
- exit(1);
- }
- }
- }
- cout << ".";
- ntotRefMap=nummap+1;
- cout << "Particle Maps read from Standard File \nTotal Number of particles: " << ntotRefMap ;
- cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
-
- if (dumpMap)
- {
- FILE* fp = fopen("maps.dump", "w+b");
- if (fp == NULL)
- {
- cout << "Error opening dump file\n";
- exit(1);
- }
- fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
- fwrite(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
- fclose(fp);
- }
- }
-
- return(0);
-}
+#include <fstream>
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <cstring>
+
+#include "map.h"
+
+using namespace std;
+
+int bioem_RefMap::readRefMaps()
+{
+ /**************************************************************************************/
+ /***********************Reading reference Particle Maps************************/
+ /**************************************************************************************/
+ if (loadMap)
+ {
+ FILE* fp = fopen("maps.dump", "rb");
+ if (fp == NULL)
+ {
+ cout << "Error opening dump file\n";
+ exit(1);
+ }
+ fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
+ if (ntotRefMap > BIOEM_MAX_MAPS)
+ {
+ cout << "BIOEM_MAX_MAPS too small\n";
+ exit(1);
+ }
+ fread(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
+ fclose(fp);
+
+ cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
+ cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
+ }
+ else
+ {
+ int nummap=-1;
+ ifstream input(filemap);
+ if (!input.good())
+ {
+ cout << "Particle Maps Failed to open file" << endl ;
+ exit(1);
+ }
+
+ char line[512] = {' '};
+ char tmpLine[512] = {' '};
+
+ while (!input.eof())
+ {
+ input.getline(line,512);
+
+ strncpy(tmpLine,line,strlen(line));
+ char *token = strtok(tmpLine," ");
+
+ if (strcmp(token,"PARTICLE")==0) // to count the number of maps
+ {
+ nummap++;
+ if (nummap % 128 == 0)
+ {
+ cout << "..." << nummap << "\n";
+ }
+ if (nummap == BIOEM_MAX_MAPS)
+ {
+ cout << "BIOEM_MAX_MAPS too small\n";
+ exit(1);
+ }
+ }
+ else
+ {
+ int i,j;
+ myfloat_t z;
+
+ char tmpVals[36] = {' '};
+
+ strncpy (tmpVals,line,8);
+ sscanf (tmpVals,"%d",&i);
+
+ strncpy (tmpVals,line+8,8);
+ sscanf (tmpVals,"%d",&j);
+
+ strncpy (tmpVals,line+16,12);
+ sscanf (tmpVals,"%f",&z);
+ //checking for Map limits
+ if(i>0 && i-1 <BIOEM_MAP_SIZE_X&& j>0 && j-1<BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
+ {
+ Ref[nummap].points[i-1][j-1]=z;
+ }
+ else
+ {
+ cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
+ exit(1);
+ }
+ }
+ }
+ cout << ".";
+ ntotRefMap=nummap+1;
+ cout << "Particle Maps read from Standard File \nTotal Number of particles: " << ntotRefMap ;
+ cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
+
+ if (dumpMap)
+ {
+ FILE* fp = fopen("maps.dump", "w+b");
+ if (fp == NULL)
+ {
+ cout << "Error opening dump file\n";
+ exit(1);
+ }
+ fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
+ fwrite(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
+ fclose(fp);
+ }
+ }
+
+ return(0);
+}
diff --git a/model.cpp b/model.cpp
index a86d69e..95a7bcb 100644
--- a/model.cpp
+++ b/model.cpp
@@ -1,236 +1,236 @@
-#include <fstream>
-#include <iostream>
-#include <stdio.h>
-#include <stdlib.h>
-#include <cstring>
-
-#include "model.h"
-
-using namespace std;
-
-
-int bioem_model::readModel()
-{
- /**************************************************************************************/
- /***************Reading reference Models either PDB or x,y,z,r,d format****************/
- /**************************************************************************************/
-
- ofstream exampleReadCoor;
- exampleReadCoor.open ("COORDREAD");
-
- std::ifstream input(filemodel);
- if(readPDB)
- {
- ifstream input(filemodel);
- if (!input.good())
- {
- cout << "PDB Failed to open file" << endl ; // pdbfilename << " ("<<filename<<")\n";
- exit(0);
- }
-
- char line[512] = {' '};
- char tmpLine[512] = {' '};
- int numres=0;
-
- // cout << " HERE " << filemodel ;
- // for eachline in the file
- while (!input.eof())
- {
- input.getline(line,512);
-
- strncpy(tmpLine,line,strlen(line));
- char *token = strtok(tmpLine," ");
-
- if (strcmp(token,"ATOM")==0) // Optional,Mandatory if standard residues exist
- {
- /*
- 1-6 "ATOM "
- 7 - 11 Integer serial Atom serial number.
- 13 - 16 Atom name Atom name.
- 17 Character altLoc Alternate location indicator.
- 18 - 20 Residue name resName Residue name.
- 22 Character chainID Chain identifier.
- 23 - 26 Integer resSeq Residue sequence number.
- 27 AChar iCode Code for insertion of residues.
- 31 - 38 Real(8.3) x Orthogonal coordinates for X in
- 39 - 46 Real(8.3) y Orthogonal coordinates for Y in
- 47 - 54 Real(8.3) z Orthogonal coordinates for Z in
- */
-
- char name[5] = {"PPP"};
- char resName[3] = {' '};
- float x=0.0;
- float y=0.0;
- float z=0.0;
- char tmp[6] = {' '};
-
-
- // parse name
- strncpy(tmp,line+12,4);
- sscanf (tmp,"%s",name);
-
- // parse resName
- strncpy(tmp,line+17,3);
- sscanf (tmp,"%s",resName);
-
- // parse x, y, z
- char tmpVals[36] = {' '};
-
- strncpy (tmpVals,line+30,8);
- sscanf (tmpVals,"%f",&x);
-
- strncpy (tmpVals,line+38,8);
- sscanf (tmpVals,"%f",&y);
-
- strncpy (tmpVals,line+46,8);
- sscanf (tmpVals,"%f",&z);
-
- if (strcmp(name,"CA") == 0)
- {
- if (numres >= BIOEM_MODEL_SIZE)
- {
- cout << "BIOEM_MODEL_SIZE too small\n";
- exit(1);
- }
- //Getting residue Radius and electron density
- radiusPointsModel[numres]=getAminoAcidRad(resName);
- densityPointsModel[numres]=getAminoAcidDensity(resName);
-
- //Getting the coordinates
- PointsModel[numres].pos[0]=x;
- PointsModel[numres].pos[1]=y;
- PointsModel[numres].pos[2]=z;
- exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " << PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres] << "\n";
- numres++;
- }
- }
-
-
- }
- nPointsModel=numres;
- cout << "Protein structure read from PDB\nTotal Number of Residues " << nPointsModel;
- cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
- }
- else //Reading model from FILE FORMAT x,y,z,rad,density
- {
- char line[128];
- int numres=0;
- FILE *file = fopen ( filemodel , "r" );
- if ( file != NULL )
- {
- while ( fgets ( line, sizeof line, file ) != NULL )
- {
- if (numres >= BIOEM_MODEL_SIZE)
- {
- cout << "BIOEM_MODEL_SIZE too small\n";
- exit(1);
- }
- sscanf(line, "%f %f %f %f %f",&PointsModel[numres].pos[0],&PointsModel[numres].pos[1],&PointsModel[numres].pos[2],&radiusPointsModel[numres],&densityPointsModel[numres]);
- exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " << PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres]<< "\n";
- numres++;
- }
- nPointsModel=numres;
- cout << "Protein structure read from Standard File \nTotal Number of Residues " << nPointsModel ;
- cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
- }
- }
- exampleReadCoor.close();
-
- //Moving to Model to its center of mass:
- myfloat3_t r_cm;
-
- for(int n=0; n < 3; n++)r_cm.pos[n] = 0.0;
-
- for(int n=0; n < nPointsModel; n++)
- {
- r_cm.pos[0] += PointsModel[n].pos[0];
- r_cm.pos[1] += PointsModel[n].pos[1];
- r_cm.pos[2] += PointsModel[n].pos[2];
- }
- r_cm.pos[0]=r_cm.pos[0]/float(nPointsModel);
- r_cm.pos[1]=r_cm.pos[1]/float(nPointsModel);
- r_cm.pos[2]=r_cm.pos[2]/float(nPointsModel);
- for(int n=0; n < nPointsModel; n++)
- {
- PointsModel[n].pos[0]-= r_cm.pos[0];
- PointsModel[n].pos[1]-= r_cm.pos[1];
- PointsModel[n].pos[2]-= r_cm.pos[2];
-
- }
-
- return(0);
-}
-
-
-
-myfloat_t bioem_model::getAminoAcidRad(char *name)
-{
- /*************** Function that gets the radius for each amino acid ****************/
- myfloat_t iaa=0;
-
- if(std::strcmp(name,"CYS")==0)iaa=2.75;
- if(std::strcmp(name,"PHE")==0)iaa=3.2;
- if(std::strcmp(name,"LEU")==0)iaa=3.1;
- if(std::strcmp(name,"TRP")==0)iaa=3.4;
- if(std::strcmp(name,"VAL")==0)iaa=2.95;
- if(std::strcmp(name,"ILE")==0)iaa=3.1;
- if(std::strcmp(name,"MET")==0)iaa=3.1;
- if(std::strcmp(name,"HIS")==0)iaa=3.05;
- if(std::strcmp(name,"TYR")==0)iaa=3.25;
- if(std::strcmp(name,"ALA")==0)iaa=2.5;
- if(std::strcmp(name,"GLY")==0)iaa=2.25;
- if(std::strcmp(name,"PRO")==0)iaa=2.8;
- if(std::strcmp(name,"ASN")==0)iaa=2.85;
- if(std::strcmp(name,"THR")==0)iaa=2.8;
- if(std::strcmp(name,"SER")==0)iaa=2.6;
- if(std::strcmp(name,"ARG")==0)iaa=3.3;
- if(std::strcmp(name,"GLN")==0)iaa=3.0;
- if(std::strcmp(name,"ASP")==0)iaa=2.8;
- if(std::strcmp(name,"LYS")==0)iaa=3.2;
- if(std::strcmp(name,"GLU")==0)iaa=2.95;
-
- if(iaa == 0)
- {
- cout << "PROBLEM WITH AMINO ACID " << name << endl;
- exit(0);
- }
- return iaa;
-
-}
-
-
-myfloat_t bioem_model::getAminoAcidDensity(char *name)
-{
- /*************** Function that gets the number of electrons for each amino acid ****************/
- myfloat_t iaa=0.0;
-
- if(std::strcmp(name,"CYS")==0)iaa=64.0;
- if(std::strcmp(name,"PHE")==0)iaa=88.0;
- if(std::strcmp(name,"LEU")==0)iaa=72.0;
- if(std::strcmp(name,"TRP")==0)iaa=108.0;
- if(std::strcmp(name,"VAL")==0)iaa=64.0;
- if(std::strcmp(name,"ILE")==0)iaa=72.0;
- if(std::strcmp(name,"MET")==0)iaa=80.0;
- if(std::strcmp(name,"HIS")==0)iaa=82.0;
- if(std::strcmp(name,"TYR")==0)iaa=96.0;
- if(std::strcmp(name,"ALA")==0)iaa=48.0;
- if(std::strcmp(name,"GLY")==0)iaa=40.0;
- if(std::strcmp(name,"PRO")==0)iaa=62.0;
- if(std::strcmp(name,"ASN")==0)iaa=66.0;
- if(std::strcmp(name,"THR")==0)iaa=64.0;
- if(std::strcmp(name,"SER")==0)iaa=56.0;
- if(std::strcmp(name,"ARG")==0)iaa=93.0;
- if(std::strcmp(name,"GLN")==0)iaa=78.0;
- if(std::strcmp(name,"ASP")==0)iaa=59.0;
- if(std::strcmp(name,"LYS")==0)iaa=79.0;
- if(std::strcmp(name,"GLU")==0)iaa=53.0;
-
- if(iaa == 0.0)
- {
- cout << "PROBLEM WITH AMINO ACID " << name << endl;
- exit(0);
- }
- return iaa;
-
-}
-
+#include <fstream>
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <cstring>
+
+#include "model.h"
+
+using namespace std;
+
+
+int bioem_model::readModel()
+{
+ /**************************************************************************************/
+ /***************Reading reference Models either PDB or x,y,z,r,d format****************/
+ /**************************************************************************************/
+
+ ofstream exampleReadCoor;
+ exampleReadCoor.open ("COORDREAD");
+
+ std::ifstream input(filemodel);
+ if(readPDB)
+ {
+ ifstream input(filemodel);
+ if (!input.good())
+ {
+ cout << "PDB Failed to open file" << endl ; // pdbfilename << " ("<<filename<<")\n";
+ exit(0);
+ }
+
+ char line[512] = {' '};
+ char tmpLine[512] = {' '};
+ int numres=0;
+
+ // cout << " HERE " << filemodel ;
+ // for eachline in the file
+ while (!input.eof())
+ {
+ input.getline(line,512);
+
+ strncpy(tmpLine,line,strlen(line));
+ char *token = strtok(tmpLine," ");
+
+ if (strcmp(token,"ATOM")==0) // Optional,Mandatory if standard residues exist
+ {
+ /*
+ 1-6 "ATOM "
+ 7 - 11 Integer serial Atom serial number.
+ 13 - 16 Atom name Atom name.
+ 17 Character altLoc Alternate location indicator.
+ 18 - 20 Residue name resName Residue name.
+ 22 Character chainID Chain identifier.
+ 23 - 26 Integer resSeq Residue sequence number.
+ 27 AChar iCode Code for insertion of residues.
+ 31 - 38 Real(8.3) x Orthogonal coordinates for X in
+ 39 - 46 Real(8.3) y Orthogonal coordinates for Y in
+ 47 - 54 Real(8.3) z Orthogonal coordinates for Z in
+ */
+
+ char name[5] = {"PPP"};
+ char resName[3] = {' '};
+ float x=0.0;
+ float y=0.0;
+ float z=0.0;
+ char tmp[6] = {' '};
+
+
+ // parse name
+ strncpy(tmp,line+12,4);
+ sscanf (tmp,"%s",name);
+
+ // parse resName
+ strncpy(tmp,line+17,3);
+ sscanf (tmp,"%s",resName);
+
+ // parse x, y, z
+ char tmpVals[36] = {' '};
+
+ strncpy (tmpVals,line+30,8);
+ sscanf (tmpVals,"%f",&x);
+
+ strncpy (tmpVals,line+38,8);
+ sscanf (tmpVals,"%f",&y);
+
+ strncpy (tmpVals,line+46,8);
+ sscanf (tmpVals,"%f",&z);
+
+ if (strcmp(name,"CA") == 0)
+ {
+ if (numres >= BIOEM_MODEL_SIZE)
+ {
+ cout << "BIOEM_MODEL_SIZE too small\n";
+ exit(1);
+ }
+ //Getting residue Radius and electron density
+ radiusPointsModel[numres]=getAminoAcidRad(resName);
+ densityPointsModel[numres]=getAminoAcidDensity(resName);
+
+ //Getting the coordinates
+ PointsModel[numres].pos[0]=x;
+ PointsModel[numres].pos[1]=y;
+ PointsModel[numres].pos[2]=z;
+ exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " << PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres] << "\n";
+ numres++;
+ }
+ }
+
+
+ }
+ nPointsModel=numres;
+ cout << "Protein structure read from PDB\nTotal Number of Residues " << nPointsModel;
+ cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
+ }
+ else //Reading model from FILE FORMAT x,y,z,rad,density
+ {
+ char line[128];
+ int numres=0;
+ FILE *file = fopen ( filemodel , "r" );
+ if ( file != NULL )
+ {
+ while ( fgets ( line, sizeof line, file ) != NULL )
+ {
+ if (numres >= BIOEM_MODEL_SIZE)
+ {
+ cout << "BIOEM_MODEL_SIZE too small\n";
+ exit(1);
+ }
+ sscanf(line, "%f %f %f %f %f",&PointsModel[numres].pos[0],&PointsModel[numres].pos[1],&PointsModel[numres].pos[2],&radiusPointsModel[numres],&densityPointsModel[numres]);
+ exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " << PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres]<< "\n";
+ numres++;
+ }
+ nPointsModel=numres;
+ cout << "Protein structure read from Standard File \nTotal Number of Residues " << nPointsModel ;
+ cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
+ }
+ }
+ exampleReadCoor.close();
+
+ //Moving to Model to its center of mass:
+ myfloat3_t r_cm;
+
+ for(int n=0; n < 3; n++)r_cm.pos[n] = 0.0;
+
+ for(int n=0; n < nPointsModel; n++)
+ {
+ r_cm.pos[0] += PointsModel[n].pos[0];
+ r_cm.pos[1] += PointsModel[n].pos[1];
+ r_cm.pos[2] += PointsModel[n].pos[2];
+ }
+ r_cm.pos[0]=r_cm.pos[0]/float(nPointsModel);
+ r_cm.pos[1]=r_cm.pos[1]/float(nPointsModel);
+ r_cm.pos[2]=r_cm.pos[2]/float(nPointsModel);
+ for(int n=0; n < nPointsModel; n++)
+ {
+ PointsModel[n].pos[0]-= r_cm.pos[0];
+ PointsModel[n].pos[1]-= r_cm.pos[1];
+ PointsModel[n].pos[2]-= r_cm.pos[2];
+
+ }
+
+ return(0);
+}
+
+
+
+myfloat_t bioem_model::getAminoAcidRad(char *name)
+{
+ /*************** Function that gets the radius for each amino acid ****************/
+ myfloat_t iaa=0;
+
+ if(std::strcmp(name,"CYS")==0)iaa=2.75;
+ if(std::strcmp(name,"PHE")==0)iaa=3.2;
+ if(std::strcmp(name,"LEU")==0)iaa=3.1;
+ if(std::strcmp(name,"TRP")==0)iaa=3.4;
+ if(std::strcmp(name,"VAL")==0)iaa=2.95;
+ if(std::strcmp(name,"ILE")==0)iaa=3.1;
+ if(std::strcmp(name,"MET")==0)iaa=3.1;
+ if(std::strcmp(name,"HIS")==0)iaa=3.05;
+ if(std::strcmp(name,"TYR")==0)iaa=3.25;
+ if(std::strcmp(name,"ALA")==0)iaa=2.5;
+ if(std::strcmp(name,"GLY")==0)iaa=2.25;
+ if(std::strcmp(name,"PRO")==0)iaa=2.8;
+ if(std::strcmp(name,"ASN")==0)iaa=2.85;
+ if(std::strcmp(name,"THR")==0)iaa=2.8;
+ if(std::strcmp(name,"SER")==0)iaa=2.6;
+ if(std::strcmp(name,"ARG")==0)iaa=3.3;
+ if(std::strcmp(name,"GLN")==0)iaa=3.0;
+ if(std::strcmp(name,"ASP")==0)iaa=2.8;
+ if(std::strcmp(name,"LYS")==0)iaa=3.2;
+ if(std::strcmp(name,"GLU")==0)iaa=2.95;
+
+ if(iaa == 0)
+ {
+ cout << "PROBLEM WITH AMINO ACID " << name << endl;
+ exit(0);
+ }
+ return iaa;
+
+}
+
+
+myfloat_t bioem_model::getAminoAcidDensity(char *name)
+{
+ /*************** Function that gets the number of electrons for each amino acid ****************/
+ myfloat_t iaa=0.0;
+
+ if(std::strcmp(name,"CYS")==0)iaa=64.0;
+ if(std::strcmp(name,"PHE")==0)iaa=88.0;
+ if(std::strcmp(name,"LEU")==0)iaa=72.0;
+ if(std::strcmp(name,"TRP")==0)iaa=108.0;
+ if(std::strcmp(name,"VAL")==0)iaa=64.0;
+ if(std::strcmp(name,"ILE")==0)iaa=72.0;
+ if(std::strcmp(name,"MET")==0)iaa=80.0;
+ if(std::strcmp(name,"HIS")==0)iaa=82.0;
+ if(std::strcmp(name,"TYR")==0)iaa=96.0;
+ if(std::strcmp(name,"ALA")==0)iaa=48.0;
+ if(std::strcmp(name,"GLY")==0)iaa=40.0;
+ if(std::strcmp(name,"PRO")==0)iaa=62.0;
+ if(std::strcmp(name,"ASN")==0)iaa=66.0;
+ if(std::strcmp(name,"THR")==0)iaa=64.0;
+ if(std::strcmp(name,"SER")==0)iaa=56.0;
+ if(std::strcmp(name,"ARG")==0)iaa=93.0;
+ if(std::strcmp(name,"GLN")==0)iaa=78.0;
+ if(std::strcmp(name,"ASP")==0)iaa=59.0;
+ if(std::strcmp(name,"LYS")==0)iaa=79.0;
+ if(std::strcmp(name,"GLU")==0)iaa=53.0;
+
+ if(iaa == 0.0)
+ {
+ cout << "PROBLEM WITH AMINO ACID " << name << endl;
+ exit(0);
+ }
+ return iaa;
+
+}
+
--
GitLab