1 """@namespace IMP.isd.gmm_tools
2 Tools for handling Gaussian Mixture Models.
5 from __future__
import print_function
16 import sklearn.mixture
20 from math
import exp,sqrt,copysign
28 """ read the output from write_gmm_to_text, decorate as Gaussian and Mass"""
30 with open(in_fn,
'r') as inf:
37 weight=float(fields[2])
38 center=list(map(float,fields[3].split()))
39 covar=np.array(list(map(float,
40 fields[4].split()))).reshape((3,3))
53 rmax = sqrt(max(g.get_variances()))*radius_scale
58 if not transform
is None:
63 """write a list of gaussians to text. must be decorated as Gaussian and Mass"""
64 print(
'will write GMM text to',out_fn)
65 with open(out_fn,
'w')
as outf:
66 outf.write(
'#|num|weight|mean|covariance matrix|\n')
67 for ng,g
in enumerate(ps):
71 mean=list(shape.get_center())
72 fm=[ng,weight]+mean+covar
75 outf.write(
'|{}|{}|{} {} {}|{} {} {} {} {} {} {} {} {}|\n'.format(*fm))
78 outf.write(
'|{0}|{1}|{2} {3} {4}|{5} {6} {7} {8} {9} {10} {11} {12} {13}|\n'.format(*fm))
80 def gmm2map(to_draw,voxel_size,bounding_box=None,origin=None, fast=False, factor=2.5):
84 ps=[g.get_particle()
for g
in to_draw]
86 print(
'ps must be Particles or Gaussians')
88 if bounding_box
is None:
106 print(
'creating map')
108 if origin
is not None:
109 d1.set_origin(origin)
111 def write_gmm_to_map(to_draw,out_fn,voxel_size,bounding_box=None,origin=None, fast=False, factor=2.5):
112 """write density map from GMM. input can be either particles or gaussians"""
113 d1 = gmm2map(to_draw,voxel_size,bounding_box,origin, fast)
114 print(
'will write GMM map to',out_fn)
119 """write density map directly from sklearn GMM (kinda slow) """
121 if not dmap_model
is None:
127 print(
'getting coords')
128 nvox=d1.get_number_of_voxels()
129 apos=[list(d1.get_location_by_voxel(nv))
for nv
in range(nvox)]
132 scores=gmm.score(apos)
135 for nv, score
in enumerate(scores):
136 d1.set_value(nv,exp(score))
137 print(
'will write GMM map to',out_fn)
140 def draw_points(pts,out_fn,trans=IMP.algebra.get_identity_transformation_3d(),
142 """ given some points (and optional transform), write them to chimera 'bild' format
143 colors flag only applies to ellipses, otherwise it'll be weird"""
144 with open(out_fn,
'w')
as outf:
150 outf.write(
'.color 1 0 0\n.dotat %.2f %.2f %.2f\n'
151 %(pt[0],pt[1],pt[2]))
156 outf.write(
'.color 0 1 0\n')
157 colors=[
'0 1 0',
'0 0 1',
'0 1 1']
158 for nt,t
in enumerate(pts[start:]):
159 if use_colors
and nt%2==0:
160 outf.write(
'.color %s\n' % colors[nt/2])
162 outf.write(
'.dotat %.2f %.2f %.2f\n' %(pt[0],pt[1],pt[2]))
166 def create_gmm_for_bead(mdl,
169 sampled_points=100000,
171 print(
'fitting bead with',n_components,
'gaussians')
174 IMP.em.write_map(dmap,
'test_intermed.mrc')
175 pts=IMP.isd.sample_points_from_density(dmap,sampled_points)
182 return density_particles,dmap
184 def sample_and_fit_to_particles(model,
187 sampled_points=1000000,
191 covariance_type=
'full',
192 multiply_by_total_mass=
True,
196 if multiply_by_total_mass:
198 print(
'add_component_density: will multiply by mass',mass_multiplier)
201 print(
'add_component_density: sampling points')
207 pts=IMP.isd.sample_points_from_density(dmap,sampled_points)
210 print(
'add_component_density: fitting GMM to',len(pts),
'points')
212 n_components=num_components,
214 ps=density_particles,
216 covariance_type=covariance_type,
217 mass_multiplier=mass_multiplier)
219 if not output_txt
is None:
221 if not output_map
is None:
224 voxel_size=voxel_size,
227 return density_particles
234 covariance_type=
'full',
239 mass_multiplier=1.0):
240 """fit a GMM to some points. Will return the score and the Akaike score.
241 Akaike information criterion for the current model fit. It is a measure
242 of the relative quality of the GMM that takes into account the
243 parsimony and the goodness of the fit.
244 if no particles are provided, they will be created
246 points: list of coordinates (python)
247 n_components: number of gaussians to create
249 ps: list of particles to be decorated. if empty, will add
250 num_iter: number of EM iterations
251 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'spherical'
252 min_covar: assign a minimum value to covariance term. That is used to have more spherical
254 init_centers: initial coordinates of the GMM
255 force_radii: fix the radii (spheres only)
256 force_weight: fix the weights
257 mass_multiplier: multiply the weights of all the gaussians by this value
258 dirichlet: use the DGMM fitting (can reduce number of components, takes longer)
262 import sklearn.mixture
266 if force_radii==-1.0:
270 covariance_type=
'spherical'
271 print(
'forcing spherical with radii',force_radii)
273 if force_weight==-1.0:
277 print(
'forcing weights to be',force_weight)
279 print(
'creating GMM with params',params,
'init params',init_params,
'n_components',n_components,
'n_iter',num_iter,
'covar type',covariance_type)
280 gmm=sklearn.mixture.GMM(n_components=n_components,
282 covariance_type=covariance_type,
285 init_params=init_params)
287 if force_weight!=-1.0:
288 gmm.weights_=np.array([force_weight]*n_components)
289 if force_radii!=-1.0:
290 gmm.covars_=np.array([[force_radii]*3
for i
in range(n_components)])
292 gmm.means_=init_centers
294 model=gmm.fit(points)
295 score=gmm.score(points)
296 akaikescore=model.aic(points)
300 for ng
in range(n_components):
301 covar=gmm.covars_[ng]
303 covar=np.diag(covar).tolist()
306 center=list(gmm.means_[ng])
307 weight=mass_multiplier*gmm.weights_[ng]
315 return (score,akaikescore)
322 covariance_type=
'full',
323 mass_multiplier=1.0):
324 """fit a GMM to some points. Will return core::Gaussians.
325 if no particles are provided, they will be created
327 points: list of coordinates (python)
328 n_components: number of gaussians to create
330 ps: list of particles to be decorated. if empty, will add
331 num_iter: number of EM iterations
332 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'spherical'
333 init_centers: initial coordinates of the GMM
334 force_radii: fix the radii (spheres only)
335 force_weight: fix the weights
336 mass_multiplier: multiply the weights of all the gaussians by this value
340 import sklearn.mixture
343 print(
'using dirichlet prior')
344 gmm=sklearn.mixture.DPGMM(n_components=n_components,
346 covariance_type=covariance_type)
356 for ng
in range(n_components):
357 invcovar=gmm.precs_[ng]
358 covar=np.linalg.inv(invcovar)
360 covar=np.diag(covar).tolist()
363 center=list(gmm.means_[ng])
364 weight=mass_multiplier*gmm.weights_[ng]
static Gaussian setup_particle(Model *m, ParticleIndex pi)
static FloatKey get_mass_key()
DenseGrid3D< double > get_rasterized(const Gaussian3Ds &gmm, const Floats &weights, double cell_width, const BoundingBox3D &bb)
Rasterize the Gaussians to a grid.
BoundingBoxD< 3 > get_bounding_box(const Geometry &)
Compute the bounding box of any geometric object.
static XYZR setup_particle(Model *m, ParticleIndex pi)
double get_mass(ResidueType c)
Get the mass from the residue type.
IMP_Eigen::Matrix3d get_covariance(const Gaussian3D &g)
static bool get_is_setup(Model *m, ParticleIndex pi)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Class for sampling a density map from particles.
The standard decorator for manipulating molecular structures.
DensityMap * create_density_map(const IMP::algebra::GridD< 3, S, V, E > &arg)
Create a density map from an arbitrary IMP::algebra::GridD.
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
Basic utilities for handling cryo-electron microscopy 3D density maps.
A decorator for a particle with x,y,z coordinates.
algebra::BoundingBoxD< 3 > get_bounding_box(const DensityMap *m)
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...
Gaussian3D get_gaussian_from_covariance(const IMP_Eigen::Matrix3d &covariance, const Vector3D ¢er)
Return a Gaussian centered at the origin from a covariance matrix.
Class to handle individual particles of a Model object.
Sphere3D get_enclosing_sphere(const Vector3Ds &ss)
Return a sphere containing the listed vectors.
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.
A decorator for helping deal with a generalized hierarchy.
A decorator for a rigid body.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Functionality for loading, creating, manipulating and scoring atomic structures.
A decorator for a particle with x,y,z coordinates and a radius.