Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
On Thursday, 7th July from 1 to 3 pm there will be a maintenance with a short downtime of GitLab.
Open sidebar
MPIBP-Hummer
BioEM
Commits
0392963f
Commit
0392963f
authored
Apr 16, 2014
by
Pilar Cossio
Browse files
Included Input Checks
parent
92e7d74a
Changes
5
Hide whitespace changes
Inline
Side-by-side
bioem.cpp
View file @
0392963f
...
@@ -47,11 +47,11 @@ bioem::~bioem()
...
@@ -47,11 +47,11 @@ bioem::~bioem()
int
bioem
::
configure
(
int
ac
,
char
*
av
[])
int
bioem
::
configure
(
int
ac
,
char
*
av
[])
{
{
/**************************************************************************************
/
/
/**************************************************************************************
/**** Configuration Routine using boost for extracting parameters, models and maps ****
/
/
/**** Configuration Routine using boost for extracting parameters, models and maps ****
/**************************************************************************************
/
/
/**************************************************************************************
/****** And Precalculating necessary grids, map crosscorrelations and kernels ********
/
/
/****** And Precalculating necessary grids, map crosscorrelations and kernels ********
/*************************************************************************************
/
/
/*************************************************************************************
/*** Inizialzing default variables ***/
/*** Inizialzing default variables ***/
std
::
string
infile
,
modelfile
,
mapfile
;
std
::
string
infile
,
modelfile
,
mapfile
;
...
@@ -61,19 +61,20 @@ int bioem::configure(int ac, char* av[])
...
@@ -61,19 +61,20 @@ int bioem::configure(int ac, char* av[])
RefMap
.
loadMap
=
false
;
RefMap
.
loadMap
=
false
;
RefMap
.
readMRC
=
false
;
RefMap
.
readMRC
=
false
;
RefMap
.
readMultMRC
=
false
;
RefMap
.
readMultMRC
=
false
;
/*************************************************************************************
/
/
/*************************************************************************************
cout
<<
"
++++++++++++ FROM COMMAND LINE +++++++++++
\n\n
"
;
cout
<<
"
+
++++++++++++ FROM COMMAND LINE +++++++++++
\n\n
"
;
/*************************************************************************************
/
/
/*************************************************************************************
/********************* Command line reading input with BOOST ************************
/
/
/********************* Command line reading input with BOOST ************************
try
{
try
{
po
::
options_description
desc
(
"Command line inputs"
);
po
::
options_description
desc
(
"Command line inputs"
);
desc
.
add_options
()
desc
.
add_options
()
(
"Inputfile"
,
po
::
value
<
std
::
string
>
(),
"Name of input parameter file"
)
(
"Inputfile"
,
po
::
value
<
std
::
string
>
(),
"
(Mandatory)
Name of input parameter file"
)
(
"Modelfile"
,
po
::
value
<
std
::
string
>
()
,
"Name of model file"
)
(
"Modelfile"
,
po
::
value
<
std
::
string
>
()
,
"
(Mandatory)
Name of model file"
)
(
"Particlesfile"
,
po
::
value
<
std
::
string
>
(),
"Name of paricles file"
)
(
"Particlesfile"
,
po
::
value
<
std
::
string
>
(),
"
(Mandatory)
Name of paricles file"
)
(
"ReadPDB"
,
"(Optional) If reading model file in PDB format"
)
(
"ReadPDB"
,
"(Optional) If reading model file in PDB format"
)
(
"ReadMRC"
,
"(Optional) If reading particle file in MRC format"
)
(
"ReadMRC"
,
"(Optional) If reading particle file in MRC format"
)
(
"ReadMultipleMRC"
,
"(Optional) If reading Multiple MRCs"
)
(
"ReadMultipleMRC"
,
"(Optional) If reading Multiple MRCs"
)
...
@@ -98,8 +99,9 @@ int bioem::configure(int ac, char* av[])
...
@@ -98,8 +99,9 @@ int bioem::configure(int ac, char* av[])
po
::
notify
(
vm
);
po
::
notify
(
vm
);
if
((
ac
<
6
))
{
if
((
ac
<
6
))
{
std
::
cout
<<
desc
<<
std
::
endl
;
std
::
cout
<<
desc
<<
std
::
endl
;
return
0
;
exit
(
1
);
return
0
;
}
}
if
(
vm
.
count
(
"help"
))
{
if
(
vm
.
count
(
"help"
))
{
cout
<<
"Usage: options_description [options]
\n
"
;
cout
<<
"Usage: options_description [options]
\n
"
;
...
@@ -111,7 +113,7 @@ int bioem::configure(int ac, char* av[])
...
@@ -111,7 +113,7 @@ int bioem::configure(int ac, char* av[])
{
{
cout
<<
"Input file is: "
;
cout
<<
"Input file is: "
;
cout
<<
vm
[
"Inputfile"
].
as
<
std
::
string
>
()
<<
"
\n
"
;
cout
<<
vm
[
"Inputfile"
].
as
<
std
::
string
>
()
<<
"
\n
"
;
infile
=
vm
[
"Inputfile"
].
as
<
std
::
string
>
();
infile
=
vm
[
"Inputfile"
].
as
<
std
::
string
>
();
}
}
if
(
vm
.
count
(
"Modelfile"
))
if
(
vm
.
count
(
"Modelfile"
))
{
{
...
@@ -162,23 +164,23 @@ int bioem::configure(int ac, char* av[])
...
@@ -162,23 +164,23 @@ int bioem::configure(int ac, char* av[])
return
1
;
return
1
;
}
}
/********************* Reading Parameter Input ***************************
/
/
/********************* Reading Parameter Input ***************************
// copying inputfile to param class
// copying inputfile to param class
param
.
fileinput
=
infile
.
c_str
();
param
.
fileinput
=
infile
.
c_str
();
param
.
readParameters
();
param
.
readParameters
();
/********************* Reading Model Input ******************************
/
/
/********************* Reading Model Input ******************************
// copying modelfile to model class
// copying modelfile to model class
Model
.
filemodel
=
modelfile
.
c_str
();
Model
.
filemodel
=
modelfile
.
c_str
();
Model
.
readModel
();
Model
.
readModel
();
/********************* Reading Particle Maps Input **********************
/
/
/********************* Reading Particle Maps Input **********************
/********* HERE: PROBLEM if maps dont fit on the memory!! ***************
/
/
/********* HERE: PROBLEM if maps dont fit on the memory!! ***************
// copying mapfile to ref map class
// copying mapfile to ref map class
RefMap
.
filemap
=
mapfile
.
c_str
();
RefMap
.
filemap
=
mapfile
.
c_str
();
RefMap
.
readRefMaps
(
param
);
RefMap
.
readRefMaps
(
param
);
/****************** Precalculating Necessary Stuff *********************
/
/
/****************** Precalculating Necessary Stuff *********************
precalculate
();
precalculate
();
if
(
getenv
(
"BIOEM_DEBUG_BREAK"
))
if
(
getenv
(
"BIOEM_DEBUG_BREAK"
))
...
@@ -194,9 +196,9 @@ int bioem::configure(int ac, char* av[])
...
@@ -194,9 +196,9 @@ int bioem::configure(int ac, char* av[])
int
bioem
::
precalculate
()
int
bioem
::
precalculate
()
{
{
/**************************************************************************************
/
/
/**************************************************************************************
/* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels *
/
/
/* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels *
/**************************************************************************************
/
/
/**************************************************************************************
// Generating Grids of orientations
// Generating Grids of orientations
param
.
CalculateGridsParam
();
param
.
CalculateGridsParam
();
...
@@ -224,12 +226,12 @@ int bioem::precalculate()
...
@@ -224,12 +226,12 @@ int bioem::precalculate()
int
bioem
::
run
()
int
bioem
::
run
()
{
{
/**************************************************************************************
/
/
/**************************************************************************************
/**** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
/
/
/**** 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); ******
/
/
/**** If we want to control the number of threads -> omp_set_num_threads(XX); ******
/****************** Declarying class of Probability Pointer *************************
/
/
/****************** Declarying class of Probability Pointer *************************
pProb
=
new
bioem_Probability
[
RefMap
.
ntotRefMap
];
pProb
=
new
bioem_Probability
[
RefMap
.
ntotRefMap
];
printf
(
"
\t
Initializing
\n
"
);
printf
(
"
\t
Initializing
\n
"
);
...
@@ -245,12 +247,12 @@ int bioem::run()
...
@@ -245,12 +247,12 @@ int bioem::run()
pProb
[
iRefMap
].
ConstAngle
[
iOrient
]
=-
99999999
;
pProb
[
iRefMap
].
ConstAngle
[
iOrient
]
=-
99999999
;
}
}
}
}
/**************************************************************************************
/
/
/**************************************************************************************
deviceStartRun
();
deviceStartRun
();
/******************************** MAIN CYCLE ******************************************
/
/
/******************************** MAIN CYCLE ******************************************
/*** Declaring Private variables for each thread *****
/
/
/*** Declaring Private variables for each thread *****
mycomplex_t
*
proj_mapFFT
;
mycomplex_t
*
proj_mapFFT
;
bioem_map
conv_map
;
bioem_map
conv_map
;
mycomplex_t
*
conv_mapFFT
;
mycomplex_t
*
conv_mapFFT
;
...
@@ -266,25 +268,25 @@ int bioem::run()
...
@@ -266,25 +268,25 @@ int bioem::run()
printf
(
"
\t
Inner 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
)));
printf
(
"
\t
Inner Loop Count (%d %d %d) %lld
\n
"
,
param
.
param_device
.
maxDisplaceCenter
,
param
.
param_device
.
GridSpaceCenter
,
param
.
param_device
.
NumberPixels
,
(
long
long
int
)
(
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberPixels
*
(
2
*
param
.
param_device
.
maxDisplaceCenter
/
param
.
param_device
.
GridSpaceCenter
+
1
)
*
(
2
*
param
.
param_device
.
maxDisplaceCenter
/
param
.
param_device
.
GridSpaceCenter
+
1
)));
for
(
int
iProjectionOut
=
0
;
iProjectionOut
<
param
.
nTotGridAngles
;
iProjectionOut
++
)
for
(
int
iProjectionOut
=
0
;
iProjectionOut
<
param
.
nTotGridAngles
;
iProjectionOut
++
)
{
{
/***************************************************************************************
/
/
/***************************************************************************************
/***** Creating Projection for given orientation and transforming to Fourier space *****
/
/
/***** Creating Projection for given orientation and transforming to Fourier space *****
timer
.
ResetStart
();
timer
.
ResetStart
();
createProjection
(
iProjectionOut
,
proj_mapFFT
);
createProjection
(
iProjectionOut
,
proj_mapFFT
);
printf
(
"Time Projection %d: %f
\n
"
,
iProjectionOut
,
timer
.
GetCurrentElapsedTime
());
printf
(
"Time Projection %d: %f
\n
"
,
iProjectionOut
,
timer
.
GetCurrentElapsedTime
());
/***************************************************************************************
/
/
/***************************************************************************************
/***** **** Internal Loop over convolutions **** *****
/
/
/***** **** Internal Loop over convolutions **** *****
for
(
int
iConv
=
0
;
iConv
<
param
.
nTotCTFs
;
iConv
++
)
for
(
int
iConv
=
0
;
iConv
<
param
.
nTotCTFs
;
iConv
++
)
{
{
printf
(
"
\t\t
Convolution %d %d
\n
"
,
iProjectionOut
,
iConv
);
printf
(
"
\t\t
Convolution %d %d
\n
"
,
iProjectionOut
,
iConv
);
/*** Calculating convolutions of projection map and crosscorrelations ***
/
/
/*** Calculating convolutions of projection map and crosscorrelations ***
timer
.
ResetStart
();
timer
.
ResetStart
();
createConvolutedProjectionMap
(
iProjectionOut
,
iConv
,
proj_mapFFT
,
conv_map
,
conv_mapFFT
,
sumCONV
,
sumsquareCONV
);
createConvolutedProjectionMap
(
iProjectionOut
,
iConv
,
proj_mapFFT
,
conv_map
,
conv_mapFFT
,
sumCONV
,
sumsquareCONV
);
printf
(
"Time Convolution %d %d: %f
\n
"
,
iProjectionOut
,
iConv
,
timer
.
GetCurrentElapsedTime
());
printf
(
"Time Convolution %d %d: %f
\n
"
,
iProjectionOut
,
iConv
,
timer
.
GetCurrentElapsedTime
());
/***************************************************************************************
/
/
/***************************************************************************************
/*** Comparing each calculated convoluted map with all experimental maps ***
/
/
/*** Comparing each calculated convoluted map with all experimental maps ***
timer
.
ResetStart
();
timer
.
ResetStart
();
compareRefMaps
(
iProjectionOut
,
iConv
,
conv_map
,
conv_mapFFT
,
sumCONV
,
sumsquareCONV
);
compareRefMaps
(
iProjectionOut
,
iConv
,
conv_map
,
conv_mapFFT
,
sumCONV
,
sumsquareCONV
);
...
@@ -305,9 +307,9 @@ int bioem::run()
...
@@ -305,9 +307,9 @@ int bioem::run()
deviceFinishRun
();
deviceFinishRun
();
/************* Writing Out Probabilities ***************
/
/
/************* Writing Out Probabilities ***************
/*** Angular Probability ***
/
/
/*** Angular Probability ***
// if(param.writeAngles){
// if(param.writeAngles){
ofstream
angProbfile
;
ofstream
angProbfile
;
...
@@ -319,12 +321,12 @@ int bioem::run()
...
@@ -319,12 +321,12 @@ int bioem::run()
for
(
int
iRefMap
=
0
;
iRefMap
<
RefMap
.
ntotRefMap
;
iRefMap
++
)
for
(
int
iRefMap
=
0
;
iRefMap
<
RefMap
.
ntotRefMap
;
iRefMap
++
)
{
{
/**** Total Probability ***
/
/
/**** 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
<<
" 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: "
;
outputProbFile
<<
"RefMap "
<<
iRefMap
<<
" Maximizing Param: "
;
/*** Param that maximize probability****
/
/
/*** 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
<<
(
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
[
0
]
<<
" "
;
outputProbFile
<<
param
.
angles
[
pProb
[
iRefMap
].
max_prob_orient
].
pos
[
1
]
<<
" "
;
outputProbFile
<<
param
.
angles
[
pProb
[
iRefMap
].
max_prob_orient
].
pos
[
1
]
<<
" "
;
...
@@ -336,7 +338,7 @@ int bioem::run()
...
@@ -336,7 +338,7 @@ int bioem::run()
outputProbFile
<<
pProb
[
iRefMap
].
max_prob_cent_y
;
outputProbFile
<<
pProb
[
iRefMap
].
max_prob_cent_y
;
outputProbFile
<<
"
\n
"
;
outputProbFile
<<
"
\n
"
;
/*** For individual files***
/
//angProbfile.open ("ANG_PROB_"iRefMap);
/
/*** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap);
if
(
param
.
writeAngles
)
if
(
param
.
writeAngles
)
{
{
...
@@ -375,7 +377,11 @@ int bioem::run()
...
@@ -375,7 +377,11 @@ int bioem::run()
int
bioem
::
compareRefMaps
(
int
iProjectionOut
,
int
iConv
,
const
bioem_map
&
conv_map
,
mycomplex_t
*
localmultFFT
,
myfloat_t
sumC
,
myfloat_t
sumsquareC
,
const
int
startMap
)
int
bioem
::
compareRefMaps
(
int
iProjectionOut
,
int
iConv
,
const
bioem_map
&
conv_map
,
mycomplex_t
*
localmultFFT
,
myfloat_t
sumC
,
myfloat_t
sumsquareC
,
const
int
startMap
)
{
{
if
(
FFTAlgo
)
//***************************************************************************************
//***** BioEM routine for comparing reference maps to convoluted maps *****
if
(
FFTAlgo
)
//IF using the fft algorithm
{
{
#pragma omp parallel
#pragma omp parallel
{
{
...
@@ -411,6 +417,10 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
...
@@ -411,6 +417,10 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
inline
void
bioem
::
calculateCCFFT
(
int
iRefMap
,
int
iOrient
,
int
iConv
,
myfloat_t
sumC
,
myfloat_t
sumsquareC
,
mycomplex_t
*
localConvFFT
,
mycomplex_t
*
localCCT
,
myfloat_t
*
lCC
)
inline
void
bioem
::
calculateCCFFT
(
int
iRefMap
,
int
iOrient
,
int
iConv
,
myfloat_t
sumC
,
myfloat_t
sumsquareC
,
mycomplex_t
*
localConvFFT
,
mycomplex_t
*
localCCT
,
myfloat_t
*
lCC
)
{
{
//***************************************************************************************
//***** Calculating cross correlation in FFTALGOrithm *****
const
mycomplex_t
*
RefMapFFT
=
&
RefMap
.
RefMapsFFT
[
iRefMap
*
param
.
RefMapSize
];
const
mycomplex_t
*
RefMapFFT
=
&
RefMap
.
RefMapsFFT
[
iRefMap
*
param
.
RefMapSize
];
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
;
i
++
)
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
;
i
++
)
{
{
...
@@ -425,10 +435,10 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t
...
@@ -425,10 +435,10 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t
int
bioem
::
createProjection
(
int
iMap
,
mycomplex_t
*
mapFFT
)
int
bioem
::
createProjection
(
int
iMap
,
mycomplex_t
*
mapFFT
)
{
{
/**************************************************************************************
/
/
/**************************************************************************************
/**** BioEM Create Projection routine in Euler angle predefined grid****************
/
/**** BioEM Create Projection routine in Euler angle predefined grid****************
********************* and turns projection into Fourier space **********************
/
//
********************* and turns projection into Fourier space **********************
/**************************************************************************************
/
/
/**************************************************************************************
myfloat3_t
RotatedPointsModel
[
Model
.
nPointsModel
];
myfloat3_t
RotatedPointsModel
[
Model
.
nPointsModel
];
myfloat_t
rotmat
[
3
][
3
];
myfloat_t
rotmat
[
3
][
3
];
...
@@ -442,9 +452,9 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
...
@@ -442,9 +452,9 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
beta
=
param
.
angles
[
iMap
].
pos
[
1
];
beta
=
param
.
angles
[
iMap
].
pos
[
1
];
gam
=
param
.
angles
[
iMap
].
pos
[
2
];
gam
=
param
.
angles
[
iMap
].
pos
[
2
];
/**** To see how things are going: cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***
/
/
/**** 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**********
/
/
/********** Creat Rotation with pre-defiend grid of orientations**********
rotmat
[
0
][
0
]
=
cos
(
gam
)
*
cos
(
alpha
)
-
cos
(
beta
)
*
sin
(
alpha
)
*
sin
(
gam
);
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
][
1
]
=
cos
(
gam
)
*
sin
(
alpha
)
+
cos
(
beta
)
*
cos
(
alpha
)
*
sin
(
gam
);
...
@@ -475,7 +485,7 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
...
@@ -475,7 +485,7 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
int
i
,
j
;
int
i
,
j
;
/************ Projection over the Z axis********************
/
/
/************ Projection over the Z axis********************
for
(
int
n
=
0
;
n
<
Model
.
nPointsModel
;
n
++
)
for
(
int
n
=
0
;
n
<
Model
.
nPointsModel
;
n
++
)
{
{
//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
...
@@ -485,7 +495,7 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
...
@@ -485,7 +495,7 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
localproj
[
i
*
param
.
param_device
.
NumberPixels
+
j
]
+=
Model
.
densityPointsModel
[
n
]
/
Model
.
NormDen
;
localproj
[
i
*
param
.
param_device
.
NumberPixels
+
j
]
+=
Model
.
densityPointsModel
[
n
]
/
Model
.
NormDen
;
}
}
/**** Output Just to check****
/
/
/**** Output Just to check****
if
(
iMap
==
10
)
if
(
iMap
==
10
)
{
{
ofstream
myexamplemap
;
ofstream
myexamplemap
;
...
@@ -503,8 +513,8 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
...
@@ -503,8 +513,8 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
myexampleRot
.
close
();
myexampleRot
.
close
();
}
}
/***** Converting projection to Fourier Space for Convolution later with kernel****
/
/
/***** Converting projection to Fourier Space for Convolution later with kernel****
/********** Omp Critical is necessary with FFTW*******
/
/
/********** Omp Critical is necessary with FFTW*******
myfftw_execute_dft_r2c
(
param
.
fft_plan_r2c_forward
,
localproj
,
mapFFT
);
myfftw_execute_dft_r2c
(
param
.
fft_plan_r2c_forward
,
localproj
,
mapFFT
);
return
(
0
);
return
(
0
);
...
@@ -512,18 +522,18 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
...
@@ -512,18 +522,18 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
int
bioem
::
createConvolutedProjectionMap
(
int
iMap
,
int
iConv
,
mycomplex_t
*
lproj
,
bioem_map
&
Mapconv
,
mycomplex_t
*
localmultFFT
,
myfloat_t
&
sumC
,
myfloat_t
&
sumsquareC
)
int
bioem
::
createConvolutedProjectionMap
(
int
iMap
,
int
iConv
,
mycomplex_t
*
lproj
,
bioem_map
&
Mapconv
,
mycomplex_t
*
localmultFFT
,
myfloat_t
&
sumC
,
myfloat_t
&
sumsquareC
)
{
{
/**************************************************************************************
/
/
/**************************************************************************************
/**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
/
/**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
**************** calculated Projection with convoluted precalculated Kernel**********
//
**************** calculated Projection with convoluted precalculated Kernel**********
*************** and Backtransforming it to real Space ******************************
/
//
*************** and Backtransforming it to real Space ******************************
/**************************************************************************************
/
/
/**************************************************************************************
myfloat_t
*
localconvFFT
;
myfloat_t
*
localconvFFT
;
localconvFFT
=
(
myfloat_t
*
)
myfftw_malloc
(
sizeof
(
myfloat_t
)
*
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberPixels
);
localconvFFT
=
(
myfloat_t
*
)
myfftw_malloc
(
sizeof
(
myfloat_t
)
*
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberPixels
);
mycomplex_t
*
tmp
;
mycomplex_t
*
tmp
;
tmp
=
(
mycomplex_t
*
)
myfftw_malloc
(
sizeof
(
mycomplex_t
)
*
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
);
tmp
=
(
mycomplex_t
*
)
myfftw_malloc
(
sizeof
(
mycomplex_t
)
*
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
);
/**** Multiplying FFTmap with corresponding kernel ****
/
/
/**** Multiplying FFTmap with corresponding kernel ****
const
mycomplex_t
*
refCTF
=
&
param
.
refCTF
[
iConv
*
param
.
RefMapSize
];
const
mycomplex_t
*
refCTF
=
&
param
.
refCTF
[
iConv
*
param
.
RefMapSize
];
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
;
i
++
)
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
;
i
++
)
...
@@ -536,10 +546,10 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
...
@@ -536,10 +546,10 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
//FFTW_C2R will destroy the input array, so we have to work on a copy here
//FFTW_C2R will destroy the input array, so we have to work on a copy here
memcpy
(
tmp
,
localmultFFT
,
sizeof
(
mycomplex_t
)
*
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
);
memcpy
(
tmp
,
localmultFFT
,
sizeof
(
mycomplex_t
)
*
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
);
/**** Bringing convoluted Map to real Space ****
/
/
/**** Bringing convoluted Map to real Space ****
myfftw_execute_dft_c2r
(
param
.
fft_plan_c2r_backward
,
tmp
,
localconvFFT
);
myfftw_execute_dft_c2r
(
param
.
fft_plan_c2r_backward
,
tmp
,
localconvFFT
);
/****Asigning convolution fftw_complex to bioem_map ****
/
/
/****Asigning convolution fftw_complex to bioem_map ****
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
;
i
++
)
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
;
i
++
)
{
{
for
(
int
j
=
0
;
j
<
param
.
param_device
.
NumberPixels
;
j
++
)
for
(
int
j
=
0
;
j
<
param
.
param_device
.
NumberPixels
;
j
++
)
...
@@ -548,7 +558,7 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
...
@@ -548,7 +558,7 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
}
}
}
}
/*** Calculating Cross-correlations of cal-convoluted map with its self *****
/
/
/*** Calculating Cross-correlations of cal-convoluted map with its self *****
sumC
=
0
;
sumC
=
0
;
sumsquareC
=
0
;
sumsquareC
=
0
;
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberPixels
;
i
++
)
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberPixels
;
i
++
)
...
@@ -556,14 +566,14 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
...
@@ -556,14 +566,14 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
sumC
+=
localconvFFT
[
i
];
sumC
+=
localconvFFT
[
i
];
sumsquareC
+=
localconvFFT
[
i
]
*
localconvFFT
[
i
];
sumsquareC
+=
localconvFFT
[
i
]
*
localconvFFT
[
i
];
}
}
/*** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***
/
/
/*** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***
// Normalizing
// Normalizing
myfloat_t
norm2
=
(
myfloat_t
)
(
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberPixels
);
myfloat_t
norm2
=
(
myfloat_t
)
(
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberPixels
);
myfloat_t
norm4
=
norm2
*
norm2
;
myfloat_t
norm4
=
norm2
*
norm2
;
sumC
=
sumC
/
norm2
;
sumC
=
sumC
/
norm2
;
sumsquareC
=
sumsquareC
/
norm4
;
sumsquareC
=
sumsquareC
/
norm4
;
/**** Freeing fftw_complex created (dont know if omp critical is necessary) ****
/
/
/**** Freeing fftw_complex created (dont know if omp critical is necessary) ****
myfftw_free
(
localconvFFT
);
myfftw_free
(
localconvFFT
);
myfftw_free
(
tmp
);
myfftw_free
(
tmp
);
...
@@ -572,7 +582,7 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
...
@@ -572,7 +582,7 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
int
bioem
::
calcross_cor
(
bioem_map
&
localmap
,
myfloat_t
&
sum
,
myfloat_t
&
sumsquare
)
int
bioem
::
calcross_cor
(
bioem_map
&
localmap
,
myfloat_t
&
sum
,
myfloat_t
&
sumsquare
)
{
{
/*********************** Routine to calculate Cross correlations***********************
/
/
/*********************** Routine to calculate Cross correlations***********************
sum
=
0.0
;
sum
=
0.0
;
sumsquare
=
0.0
;
sumsquare
=
0.0
;
...
...
main.cpp
View file @
0392963f
...
@@ -21,9 +21,9 @@
...
@@ -21,9 +21,9 @@
int
main
(
int
argc
,
char
*
argv
[])
int
main
(
int
argc
,
char
*
argv
[])
{
{
/
*
*************************************************************************************
/
/
/
*************************************************************************************
/********************************* Main BioEM code **********************************
/
/
/********************************* Main BioEM code **********************************
/************************************************************************************
/
/
/*
************************************************************************************
#pragma omp parallel
#pragma omp parallel
{
{
...
@@ -43,17 +43,17 @@ int main(int argc, char* argv[])
...
@@ -43,17 +43,17 @@ int main(int argc, char* argv[])
bio
=
new
bioem
;
bio
=
new
bioem
;
}
}
/************ Configuration and Pre-calculating necessary objects *****************
/
/
/************ Configuration and Pre-calculating necessary objects *****************
****
printf
(
"Configuring
\n
"
);
printf
(
"Configuring
\n
"
);
bio
->
configure
(
argc
,
argv
);
bio
->
configure
(
argc
,
argv
);
/******************************* Run BioEM routine ******************************
/
/
/******************************* Run BioEM routine ******************************
******
printf
(
"Running
\n
"
);
printf
(
"Running
\n
"
);
timer
.
Start
();
timer
.
Start
();
bio
->
run
();
bio
->
run
();
timer
.
Stop
();
timer
.
Stop
();
/************************************ End **********************************
/
/
/************************************ End **********************************
*************
printf
(
"The code ran for %f seconds.
\n
"
,
timer
.
GetElapsedTime
());
printf
(
"The code ran for %f seconds.
\n
"
,
timer
.
GetElapsedTime
());
delete
bio
;
delete
bio
;
...
...
map.cpp
View file @
0392963f
...
@@ -13,9 +13,9 @@ using namespace std;
...
@@ -13,9 +13,9 @@ using namespace std;
int
bioem_RefMap
::
readRefMaps
(
bioem_param
&
param
)
int
bioem_RefMap
::
readRefMaps
(
bioem_param
&
param
)
{
{
/**************************************************************************************
/
/
/**************************************************************************************
/***********************Reading reference Particle Maps************************
/
/
/***********************Reading reference Particle Maps************************
/**************************************************************************************
/
/
/**************************************************************************************
if
(
loadMap
)
if
(
loadMap
)
{
{
FILE
*
fp
=
fopen
(
"maps.dump"
,
"rb"
);
FILE
*
fp
=
fopen
(
"maps.dump"
,
"rb"
);
...
@@ -34,7 +34,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
...
@@ -34,7 +34,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
fclose
(
fp
);
fclose
(
fp
);
cout
<<
"Particle Maps read from Map Dump
\n
Total Number of particles: "
<<
ntotRefMap
;
cout
<<
"Particle Maps read from Map Dump
\n
Total Number of particles: "
<<
ntotRefMap
;
cout
<<
"
\n
+++++++++++++++++++++++++++++++++++++++
++
\n
"
;
cout
<<
"
\n
+++++++++++++++++++++++++++++++++++++++
\n
"
;
}
}
else
else
{
{
...
@@ -79,11 +79,13 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
...
@@ -79,11 +79,13 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
for
(
int
i
=
3
;
i
<
100
;
i
++
)
mapname
[
i
]
=
{
0
};
for
(
int
i
=
3
;
i
<
100
;
i
++
)
mapname
[
i
]
=
{
0
};
}
}
cout
<<
"
\n
+++++++++++++++++++++++++++++++++++++++++++
\n
"
;
cout
<<
"Particle Maps read from MULTIPLE MRC Files in: "
<<
filemap
<<
"
\n
"
;
cout
<<
"Particle Maps read from MULTIPLE MRC Files in: "
<<
filemap
<<
"
\n
"
;
}
else
{
}
else
{
read_MRC
(
filemap
,
param
);
read_MRC
(
filemap
,
param
);
cout
<<
"
\n
++++++++++++++++++++++++++++++++++++++++++
\n
"
;
cout
<<
"Particle Maps read from ONE MRC File: "
<<
filemap
<<
"
\n
"
;
cout
<<
"Particle Maps read from ONE MRC File: "
<<
filemap
<<
"
\n
"
;
}
}
...
@@ -160,6 +162,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
...
@@ -160,6 +162,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
}
}
cout
<<
"."
;
cout
<<
"."
;
ntotRefMap
=
nummap
+
1
;
ntotRefMap
=
nummap
+
1
;
cout
<<
"
\n
++++++++++++++++++++++++++++++++++++++++++
\n
"
;
cout
<<
"Particle Maps read from Standard File
\n
"
;
cout
<<
"Particle Maps read from Standard File
\n
"
;
}
}
cout
<<
"Total Number of particles: "
<<
ntotRefMap
;
cout
<<
"Total Number of particles: "
<<
ntotRefMap
;
...
@@ -184,9 +187,9 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
...
@@ -184,9 +187,9 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
int
bioem_RefMap
::
PreCalculateMapsFFT
(
bioem_param
&
param
)
int
bioem_RefMap
::
PreCalculateMapsFFT
(
bioem_param
&
param
)
{