1 """@namespace IMP.bayesianem
2 Restraints for handling electron microscopy maps.
5 from __future__
import print_function
20 """Fit Gaussian-decorated particles to an EM map
21 (also represented with a set of Gaussians)
22 \note This class wraps an isd::GaussianEMRestraint
27 cutoff_dist_model_model=10.0,
28 cutoff_dist_model_data=10.0,
29 target_mass_scale=1.0,
31 target_radii_scale=3.0,
32 model_radii_scale=1.0,
34 spherical_gaussians=
False,
35 close_pair_container=
None,
37 scale_target_to_mass=
False,
39 target_is_rigid_body=
False):
41 @param densities The Gaussian-decorated particles to be restrained
42 @param target_fn GMM file of the target density map
43 (alternatively, pass the ps)
44 @param target_ps List of Gaussians of the target map
45 (alternatively, pass the filename)
46 @param cutoff_dist_model_model Distance in model-model close
48 @param cutoff_dist_model_data Distance in model-data close pair
49 container. Usually can set to zero because we multiply the
51 @param target_mass_scale Scale up the target densities so that
53 Needed if the GMM you generated was not already scaled.
54 To make it the same as model mass, set scale_to_target_mass=True
55 @param target_mass Sets the mass of the target density to the given value. Default is None. This
56 will override target_mass_scale argument
57 @param target_radii_scale Scale the target density radii -
58 only used for the close pair container.
59 If you keep this at 3.0 or so you don't have to use cutoff dist.
60 @param model_radii_scale Scale the model density radii - only used
61 for the close pair container
62 @param slope Linear term added to help bring model into the density
63 @param spherical_gaussians Set to True for a speed bonus when
64 the model densities are spheres. (This means you don't have
65 to do a matrix multiplication if they rotate.)
66 @param close_pair_container Pass a close pair container for
67 the model if you already have one (e.g. for an excluded
68 volume restraint.) May give a speed bonus.
69 @param backbone_slope Only apply slope to backbone particles -
70 only matters for atomic
71 @param scale_target_to_mass Set True if you would need to scale
72 target to EXACTLY the model mass
73 @param weight The restraint weight
74 @param target_is_rigid_body Set True if you want to put the target density particles
75 into a rigid body that need to be sampled (e.g.,when you need to fit one density
76 against another one). Default is False.
80 self.sigmaissampled =
True
81 self.sigmamaxtrans = 0.3
82 self.sigmamin = 0.0001
86 self.densities = densities
87 self.em_root_hier =
None
90 self.m = self.densities[0].get_model()
92 if scale_target_to_mass:
93 def hierarchy_mass(h):
96 target_mass = sum(hierarchy_mass(h)
for h
in densities)
97 print(
'will scale target mass by', target_mass_scale)
100 self._set_dataset(target_fn)
106 radius_scale=target_radii_scale,
107 mass_scale=target_mass_scale)
108 elif target_ps != []:
109 self.target_ps = target_ps
111 print(
'Gaussian EM restraint: must provide target density file or properly set up target densities')
116 scale=target_mass/tmass
117 print(
'will set target mass to', target_mass,tmass,scale)
118 for p
in self.target_ps:
122 for p, state
in IMP.pmi.tools._all_protocol_outputs([], densities[0]):
123 p.add_em3d_restraint(state, self.target_ps, self.densities, self)
127 for h
in self.densities:
129 if model_radii_scale != 1.0:
130 for p
in self.model_ps:
138 if target_is_rigid_body:
146 self.sigmaglobal = IMP.pmi.tools.SetupNuisance(self.m, self.sigmainit,
147 self.sigmamin, self.sigmamax,
148 self.sigmaissampled).get_particle()
151 print(
'target num particles', len(self.target_ps), \
153 for p
in self.target_ps]))
154 print(
'model num particles', len(self.model_ps), \
156 for p
in self.model_ps]))
173 update_model=
not spherical_gaussians
179 self.sigmaglobal.get_particle().
get_index(),
180 cutoff_dist_model_model,
181 cutoff_dist_model_data,
183 update_model, backbone_slope)
185 print(
'done EM setup')
187 self.rs.add_restraint(self.gaussianEM_restraint)
190 def _set_dataset(self, target_fn):
191 """Set the dataset to point to the input file"""
193 self.dataset = p.parse_file(target_fn)[
'dataset']
197 aligns the center of mass of the target GMM on the center of mass of the model
201 for p
in self.target_ps:
204 target_com += pos * mass
206 target_com /= target_mass
207 print(
'target com', target_com)
210 for p
in self.model_ps:
213 model_com += pos * mass
215 model_com /= model_mass
216 print(
'model com', model_com)
218 v = target_com - model_com
219 print(
'translating with', -v)
221 for p
in self.target_ps:
226 '''Returns the geometric center of the GMM particles
227 @param target = True - returns target map gmm COM
228 @param target = False - returns model gmm COM'''
245 aligns the center of mass of the target GMM on the origin
248 print(
'target com', target_com)
250 print(
'model com', model_com)
251 v = target_com - model_com
252 print(
'translating with', -v)
254 for p
in self.target_ps:
260 aligns the model on the target density
261 @param input_objects IMP.pmi.representation.Representation or IMP.pmi.topology.State
263 if type(input_object)
is IMP.pmi.representation.Representation:
264 hier = input_object.prot
266 hier = input_object.get_hierarchy()
268 raise Exception(
"Input must be a Representation or topology.State object")
270 print(
'target com', target_com)
272 print(
'model com', model_com)
273 v = target_com - model_com
274 print(
'translating with', v)
287 for rb
in list(rigid_bodies):
290 for p
in list(XYZRs):
295 align the model on target GMM
298 print(
'target com', target_com)
300 print(
'model com', model_com)
301 v = target_com - model_com
302 print(
'translating with', v)
308 for p
in self.model_ps:
315 for rb
in list(rigid_bodies):
318 for p
in list(XYZRs):
324 set the weight of the restraint
328 self.rs.set_weight(weight)
330 def set_label(self, label):
333 def add_to_model(self):
336 def get_particles_to_sample(self):
338 if self.sigmaissampled:
339 ps[
"Nuisances_GaussianEMRestraint_sigma_" +
340 self.label] = ([self.sigmaglobal], self.sigmamaxtrans)
342 ps[
"Rigid_Bodies_GaussianEMRestraint_"+self.label] = (
348 def get_rigid_body(self):
350 raise Exception(
"No rigid body created for GMM particles. Ensure target_is_rigid_body is set to True")
355 returns a hierarchy whose leaves are the gaussian particles of the target GMM
357 if self.em_root_hier
is None:
359 self.em_root_hier.set_name(
"GaussianEMRestraint_density_"+self.label)
360 for p
in self.target_ps:
361 self.em_root_hier.add_child(p)
362 return self.em_root_hier
365 ''' Can add a target GMM to a Hierarchy.
366 For PMI2 a state object may also be passed'''
372 raise Exception(
"Can only add a density to a PMI State object or IMP.atom.Hierarchy. You passed a", type(inp))
374 def get_restraint_set(self):
377 def get_output(self):
380 score = self.weight * self.rs.unprotected_evaluate(
None)
381 output[
"_TotalScore"] = str(score)
382 output[
"GaussianEMRestraint_" +
383 self.label] = str(score)
384 output[
"GaussianEMRestraint_sigma_" +
385 self.label] = str(self.sigmaglobal.get_scale())
389 return self.weight * self.rs.unprotected_evaluate(
None)
392 '''Writes target GMM file to MRC'''
394 fileout=
"Gaussian_map_" + self.label +
".mrc"
Support for the mmCIF file format.
def get_center_of_mass
Returns the geometric center of the GMM particles.
Restraints for handling electron microscopy maps.
A member of a rigid body, it has internal (local) coordinates.
def center_model_on_target_density
aligns the model on the target density
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Various classes to hold sets of particles.
static XYZR setup_particle(Model *m, ParticleIndex pi)
def center_target_density_on_origin
aligns the center of mass of the target GMM on the origin
double get_mass(ResidueType c)
Get the mass from the residue type.
def add_target_density_to_hierarchy
Can add a target GMM to a Hierarchy.
Creates a restraint between two Gaussian Mixture Models, "model" and "density".
Extract metadata from an EM density GMM file.
Fit Gaussian-decorated particles to an EM map (also represented with a set of Gaussians) ...
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Object used to hold a set of restraints.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
def center_on_target_density
align the model on target GMM
A decorator for a particle with x,y,z coordinates.
static Copy setup_particle(Model *m, ParticleIndex pi, Int number)
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
def write_target_gmm_to_mrc
Writes target GMM file to MRC.
The general base class for IMP exceptions.
IMP::core::RigidBody create_rigid_body(Hierarchy h)
Class to handle individual particles of a Model object.
Stores a list of Molecules all with the same State index.
def set_weight
set the weight of the restraint
A decorator for a rigid body.
def center_target_density_on_model
aligns the center of mass of the target GMM on the center of mass of the model
Functionality for loading, creating, manipulating and scoring atomic structures.
Hierarchies get_leaves(const Selection &h)
def get_density_as_hierarchy
returns a hierarchy whose leaves are the gaussian particles of the target GMM
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
A decorator for a particle with x,y,z coordinates and a radius.