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 (
size_t i=0; i<m1.size(); ++i) {
40 return std::sqrt(rmsd / N);
48 for (
size_t 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 (
size_t 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 (
size_t i=0; i<ps1.size(); ++i) {
88 for (
size_t i=0; i<ps2.size(); ++i) {
92 return get_rmsd_of_best_population(ds1,ds2,percentage);
104 for (
size_t i=0; i<ps1.size(); ++i) {
107 for (
size_t i=0; i<ps2.size(); ++i) {
118 const double percentage){
119 std::vector<double> sq_distances(m1.size());
120 for (
size_t i=0; i<m1.size(); ++i) {
123 std::sort(sq_distances.begin(), sq_distances.end());
126 for (
int i=0; ((double) i) < (percentage*m1.size()); ++i) {
130 return std::sqrt(sd/N);
133 double gem_score_cc(Particles model_ps, Particles density_ps){
134 Eigen::Vector3d deriv;
135 double mm_score(0.0);
136 double md_score(0.0);
137 double dd_score(0.0);
139 int nm = model_ps.size();
140 int nd = density_ps.size();
143 for(
int nd1=0 ; nd1<nd ; ++nd1){
144 for(
int nd2=0 ; nd2<nd ; ++nd2){
145 dd_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(density_ps[nd1]->
get_index(), density_ps[nd2]->
get_index()), &deriv);
148 for(
int nm1=0 ; nm1<nm ; ++nm1){
149 for(
int nm2=0 ; nm2<nm ; ++nm2){
150 mm_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(model_ps[nm1]->
get_index(), model_ps[nm2]->
get_index()), &deriv);
152 for(
int nd1=0 ; nd1<nd ; ++nd1){
153 md_score += IMP::isd::score_gaussian_overlap(mdl, ParticleIndexPair(model_ps[nm1]->
get_index(), density_ps[nd1]->
get_index()), &deriv);
156 double cc = md_score/std::sqrt(mm_score*dd_score);
162 double get_gaussian_eval_prefactor(
double determinant) {
163 return 1.0 / pow(2 *
algebra::PI, 2.0 / 3.0) / std::sqrt(determinant);
165 Eigen::Vector3d get_vec_from_center(IMP::algebra::GridIndex3D i,
166 DensityGrid
const &g,
167 Eigen::Vector3d
const ¢er) {
169 Eigen::Vector3d loc(aloc[0], aloc[1], aloc[2]);
170 Eigen::Vector3d r(loc - center);
173 void update_value(DensityGrid *g,
174 DensityGrid::Index i, Eigen::Vector3d r,
175 Eigen::Matrix3d inverse,
double pre,
Float weight) {
176 double d(r.transpose() * (inverse * r));
177 double score(pre * weight * std::exp(-0.5 * (d)));
188 "X voxels don't match");
191 "Y voxels don't match");
194 "Z voxels don't match");
195 for (
unsigned int i = 0; i < ret.get_number_of_voxels(0); ++i) {
196 for (
unsigned int j = 0; j < ret.get_number_of_voxels(1); ++j) {
197 for (
unsigned int k = 0; k < ret.get_number_of_voxels(2); ++k) {
198 DensityGrid::ExtendedIndex ei(i, j, k);
199 DensityGrid::Index gi = ret.get_index(ei);
214 for (
const DensityGrid::Index & i : mask.get_all_indexes()) {
218 for (
unsigned int ng = 0; ng < gmm.size(); ng++) {
220 Eigen::Matrix3d inverse = Eigen::Matrix3d::Zero(3, 3);
224 covar.computeInverseAndDetWithCheck(inverse, determinant, invertible);
226 "\n\n\n->>>>not proper matrix!!\n\n\n");
227 double pre(get_gaussian_eval_prefactor(determinant));
228 Eigen::Vector3d evals = covar.eigenvalues().real();
229 double maxeval = sqrt(evals.maxCoeff());
230 double cutoff = 2.5 * maxeval;
231 double cutoff2 = cutoff * cutoff;
237 IMP::algebra::GridIndex3D lowerindex = mask.get_nearest_index(lower);
238 IMP::algebra::GridIndex3D upperindex = mask.get_nearest_index(upper);
239 Eigen::Vector3d center(c.get_data());
242 lowerindex, upperindex, {
243 IMP::algebra::GridIndex3D i(voxel_index[0], voxel_index[1],
245 Eigen::Vector3d r(get_vec_from_center(i, mask, center));
246 if (r.squaredNorm() < cutoff2) {
247 update_value(&mask, i, r, inverse, pre, weights[ng]);
253 return ret.release();
261 double *data1 = sub_gmm->get_data();
262 double *data2 = gmm->get_data();
263 double *new_data = m_map->get_data();
264 for (
long i = 0; i < header->get_number_of_voxels(); i++) {
265 if (data2[i] != 0.0) {
266 double const w(data1[i] / data2[i]);
272 return m_map.release();
277 double factor = 2.5) {
297 for (
const DensityGrid::Index & i : orig.get_all_indexes()) {
299 DensityGrid::Index j(rasterized.get_nearest_index(position));
301 double x = rasterized[j] - orig[i];
310 double factor = 2.5,
int Nbin=100) {
313 Floats densities(Nbin,0.);
317 double min(std::numeric_limits<double>::max());
326 if ((rho>0.0) && (rho<min)){
330 double const dx( (max-min)/(Nbin) );
333 for(
int i(0) ; i<Nbin ; ++i){
334 densities[i]=min+dx*i;
344 size_t const index = floor((rho-min)/dx);
346 double const gmm_val = rasterized_map->get_value(i);
347 double const x = gmm_val - rho;
348 score[index] += x * x;
353 for(
int i(0) ; i<Nbin ; ++i){
357 ret[i][0]=densities[i];
371 float sgn(
double v) {
372 return ( v < 0.0 ) ? -1.0 : ( ( v > 0.0 ) ? 1.0 : 0 );
375 IMP::algebra::PrincipalComponentAnalysis NormalizePCA(
const IMP::algebra::PrincipalComponentAnalysis& pca,
381 for (
unsigned int i = 0; i < 3; ++i) {
387 for (
unsigned int j = 0 ; j < ps.size() ; ++j ) {
390 double ip = x.get_scalar_product(y);
392 ips += sgn(ip) * ip * ip;
395 if ( ips < 0. ) x *= -1.0;
400 IMP::algebra::PrincipalComponentAnalysis newpcas =
401 IMP::algebra::PrincipalComponentAnalysis(newpcs, pca.get_principal_values(), pca.get_centroid());
410 const IMP::algebra::PrincipalComponentAnalysis& pca2)
426 rotations[0] = IMP::bayesianem::get_rotation_matrix(xi, yi);
427 for (
unsigned int k = 0; k < rotations.size(); k++) {
430 rotations[k], -(rotations[k] * pca1.get_centroid()));
432 inverse_map_trans * points_trans;
433 transforms.push_back(ps2dens);
462 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.
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...)
double get_value(float x, float y, float z) const
Gets the value of the voxel located at (x,y,z)
#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.