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.