8 #ifndef IMPBAYESIANEM_UTILITIES_H
9 #define IMPBAYESIANEM_UTILITIES_H
11 #include <IMP/bayesianem/bayesianem_config.h>
15 #include <Eigen/Geometry>
17 #include <Eigen/Eigenvalues>
24 IMPBAYESIANEM_BEGIN_NAMESPACE
33 for (
int i=0; i<m1.size(); ++i) {
40 return std::sqrt(rmsd / N);
48 for (
int i=0; i<m1.size(); ++i) {
54 return ((
double) N)/m1.size();
60 const double percentage){
61 std::vector<double> sq_distances(m1.size());
62 for (
int i=0; i<m1.size(); ++i) {
65 std::sort(sq_distances.begin(), sq_distances.end());
68 for (
int i=0; ((double) i) < (percentage*m1.size()); ++i) {
72 return std::sqrt(sd/N);
77 const double percentage){
85 for (
int i=0; i<ps1.size(); ++i){ds1.push_back(
IMP::core::XYZ(ps1[i]));}
86 for (
int i=0; i<ps2.size(); ++i){ds2.push_back(
IMP::core::XYZ(ps2[i]));}
88 return get_rmsd_of_best_population(ds1,ds2,percentage);
100 for (
int i=0; i<ps1.size(); ++i){ds1.push_back(
IMP::core::XYZ(ps1[i]));}
101 for (
int i=0; i<ps2.size(); ++i){ds2.push_back(
IMP::core::XYZ(ps2[i]));}
110 const double percentage){
111 std::vector<double> sq_distances(m1.size());
112 for (
int i=0; i<m1.size(); ++i) {
115 std::sort(sq_distances.begin(), sq_distances.end());
118 for (
int i=0; ((double) i) < (percentage*m1.size()); ++i) {
122 return std::sqrt(sd/N);
125 double gem_score_cc(Particles model_ps, Particles density_ps){
126 Eigen::Vector3d deriv;
127 double mm_score(0.0);
128 double md_score(0.0);
129 double dd_score(0.0);
131 int nm = model_ps.size();
132 int nd = density_ps.size();
135 for(
int nd1=0 ; nd1<nd ; ++nd1){
136 for(
int nd2=0 ; nd2<nd ; ++nd2){
137 dd_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(density_ps[nd1]->
get_index(), density_ps[nd2]->
get_index()), &deriv);
140 for(
int nm1=0 ; nm1<nm ; ++nm1){
141 for(
int nm2=0 ; nm2<nm ; ++nm2){
142 mm_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(model_ps[nm1]->
get_index(), model_ps[nm2]->
get_index()), &deriv);
144 for(
int nd1=0 ; nd1<nd ; ++nd1){
145 md_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(model_ps[nm1]->
get_index(), density_ps[nd1]->
get_index()), &deriv);
148 double cc = md_score/std::sqrt(mm_score*dd_score);
154 double get_gaussian_eval_prefactor(
double determinant) {
155 return 1.0 / pow(2 *
algebra::PI, 2.0 / 3.0) / std::sqrt(determinant);
157 Eigen::Vector3d get_vec_from_center(IMP::algebra::GridIndex3D i,
158 DensityGrid
const &g,
159 Eigen::Vector3d
const ¢er) {
161 Eigen::Vector3d loc(aloc[0], aloc[1], aloc[2]);
162 Eigen::Vector3d r(loc - center);
165 void update_value(DensityGrid *g,
166 DensityGrid::Index i, Eigen::Vector3d r,
167 Eigen::Matrix3d inverse,
double pre,
Float weight) {
168 double d(r.transpose() * (inverse * r));
169 double score(pre * weight * std::exp(-0.5 * (d)));
180 "X voxels don't match");
183 "Y voxels don't match");
186 "Z voxels don't match");
187 for (
unsigned int i = 0; i < ret.get_number_of_voxels(0); ++i) {
188 for (
unsigned int j = 0; j < ret.get_number_of_voxels(1); ++j) {
189 for (
unsigned int k = 0; k < ret.get_number_of_voxels(2); ++k) {
190 DensityGrid::ExtendedIndex ei(i, j, k);
191 DensityGrid::Index gi = ret.get_index(ei);
206 for (
const DensityGrid::Index & i : mask.get_all_indexes()) {
210 for (
unsigned int ng = 0; ng < gmm.size(); ng++) {
212 Eigen::Matrix3d inverse = Eigen::Matrix3d::Zero(3, 3);
216 covar.computeInverseAndDetWithCheck(inverse, determinant, invertible);
218 "\n\n\n->>>>not proper matrix!!\n\n\n");
219 double pre(get_gaussian_eval_prefactor(determinant));
220 Eigen::Vector3d evals = covar.eigenvalues().real();
221 double maxeval = sqrt(evals.maxCoeff());
222 double cutoff = 2.5 * maxeval;
223 double cutoff2 = cutoff * cutoff;
229 IMP::algebra::GridIndex3D lowerindex = mask.get_nearest_index(lower);
230 IMP::algebra::GridIndex3D upperindex = mask.get_nearest_index(upper);
231 Eigen::Vector3d center(c.get_data());
234 lowerindex, upperindex, {
235 IMP::algebra::GridIndex3D i(voxel_index[0], voxel_index[1],
237 Eigen::Vector3d r(get_vec_from_center(i, mask, center));
238 if (r.squaredNorm() < cutoff2) {
239 update_value(&mask, i, r, inverse, pre, weights[ng]);
245 return ret.release();
253 IMP::em::emreal *data1 = sub_gmm->get_data();
254 IMP::em::emreal *data2 = gmm->get_data();
255 IMP::em::emreal *new_data = m_map->get_data();
256 for (
long i = 0; i < header->get_number_of_voxels(); i++) {
257 if (data2[i] != 0.0) {
258 double const w(data1[i] / data2[i]);
264 return m_map.release();
269 double factor = 2.5) {
289 for (
const DensityGrid::Index & i : orig.get_all_indexes()) {
291 DensityGrid::Index j(rasterized.get_nearest_index(position));
293 double x = rasterized[j] - orig[i];
302 double factor = 2.5,
int Nbin=100) {
305 Floats densities(Nbin,0.);
309 double min(std::numeric_limits<double>::max());
318 if ((rho>0.0) && (rho<min)){
322 double const dx( (max-min)/(Nbin) );
325 for(
int i(0) ; i<Nbin ; ++i){
326 densities[i]=min+dx*i;
336 size_t const index = floor((rho-min)/dx);
338 double const gmm_val = rasterized_map->get_value(i);
339 double const x = gmm_val - rho;
340 score[index] += x * x;
345 for(
int i(0) ; i<Nbin ; ++i){
349 ret[i][0]=densities[i];
363 float sgn(
double v) {
364 return ( v < 0.0 ) ? -1.0 : ( ( v > 0.0 ) ? 1.0 : 0 );
367 IMP::algebra::PrincipalComponentAnalysis NormalizePCA(
const IMP::algebra::PrincipalComponentAnalysis& pca,
373 for (
unsigned int i = 0; i < 3; ++i) {
379 for (
unsigned int j = 0 ; j < ps.size() ; ++j ) {
382 double ip = x.get_scalar_product(y);
384 ips += sgn(ip) * ip * ip;
387 if ( ips < 0. ) x *= -1.0;
392 IMP::algebra::PrincipalComponentAnalysis newpcas =
393 IMP::algebra::PrincipalComponentAnalysis(newpcs, pca.get_principal_values(), pca.get_centroid());
402 const IMP::algebra::PrincipalComponentAnalysis& pca2)
418 rotations[0] = IMP::bayesianem::get_rotation_matrix(xi, yi);
419 for (
unsigned int k = 0; k < rotations.size(); k++) {
422 rotations[k], -(rotations[k] * pca1.get_centroid()));
424 inverse_map_trans * points_trans;
425 transforms.push_back(ps2dens);
454 IMPBAYESIANEM_END_NAMESPACE
const DensityHeader * get_header() const
Returns a read-only pointer to the header of the map.
static const double PI
the constant pi
#define IMP_CHECK_OBJECT(obj)
Perform some basic validity checks on the object for memory debugging.
#define IMP_FUNCTION_LOG
Beginning logging for a non-member function.
IMP::Vector< Float > Floats
Standard way to pass a bunch of Float values.
emreal get_value(float x, float y, float z) const
Gets the value of the voxel located at (x,y,z)
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the squared distance between two vectors.
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Class for storing model, its restraints, constraints, and particles.
Class for handling density maps.
ParticlesTemp get_selected_particles(bool with_representation=true) const
Get the selected particles.
IMP::algebra::Transformation3D get_transformation_aligning_first_to_second(const Vector3DsOrXYZs0 &source, const Vector3DsOrXYZs1 &target)
Compute the rigid transform bringing the first point set to the second.
Vector3D get_vector_product(const Vector3D &p1, const Vector3D &p2)
Return the vector product (cross product) of two vectors.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
DensityMap * multiply(const DensityMap *m1, const DensityMap *m2)
Return a density map for which voxel i contains the result of m1[i]*m2[i].
Include all non-deprecated headers in IMP.isd.
BoundingBoxD< 3 > get_bounding_box(const Cone3D &g)
A decorator for a particle with x,y,z coordinates.
DensityMap * create_density_map(const DensityMap *other)
Create a copy of another map.
const algebra::Vector3D & get_coordinates() const
Convert it to a vector.
long get_number_of_voxels() const
Get the number of map voxels.
long get_voxel_by_location(float x, float y, float z) const
Calculate the voxel of a given location.
DensityGrid get_grid(IMP::em::DensityMap *in)
Return a dense grid containing the voxels of the passed density map.
Rotation3D get_rotation_from_matrix(double m00, double m01, double m02, double m10, double m11, double m12, double m20, double m21, double m22)
Generate a Rotation3D object from a rotation matrix.
IMP::Vector< Transformation3D > Transformation3Ds
Eigen::Matrix3d get_covariance(const Gaussian3D &g)
Include all non-deprecated headers in IMP.em.
IMP::algebra::Transformation3D get_transformation_aligning_first_to_second(const IMP::atom::Selection &s1, const IMP::atom::Selection &s2)
Get the transformation to align two selections.
double Float
Basic floating-point value (could be float, double...)
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
IMP::Vector< Int > Ints
Standard way to pass a bunch of Int values.
DenseGrid3D< double > get_rasterized_fast(const Gaussian3Ds &gmm, const Floats &weights, double cell_width, const BoundingBox3D &bb, double factor=2.5)
Rasterize the Gaussians to a grid.
#define IMP_GRID3D_FOREACH_SMALLER_EXTENDED_INDEX_RANGE(grid, center, lower_corner, upper_corner, action)
Iterate over each voxel in a subset of the grid that are less than center.
Select hierarchy particles identified by the biological name.