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 for comment
in comments:
67 outf.write(
'# ' + comment +
'\n')
68 outf.write(
'#|num|weight|mean|covariance matrix|\n')
69 for ng,g
in enumerate(ps):
73 mean=list(shape.get_center())
74 fm=[ng,weight]+mean+covar
77 outf.write(
'|{}|{}|{} {} {}|{} {} {} {} {} {} {} {} {}|\n'.format(*fm))
80 outf.write(
'|{0}|{1}|{2} {3} {4}|{5} {6} {7} {8} {9} {10} {11} {12} {13}|\n'.format(*fm))
82 def gmm2map(to_draw,voxel_size,bounding_box=None,origin=None, fast=False, factor=2.5):
86 ps=[g.get_particle()
for g
in to_draw]
88 print(
'ps must be Particles or Gaussians')
90 if bounding_box
is None:
108 print(
'creating map')
110 if origin
is not None:
111 d1.set_origin(origin)
113 def write_gmm_to_map(to_draw,out_fn,voxel_size,bounding_box=None,origin=None, fast=False, factor=2.5):
114 """write density map from GMM. input can be either particles or gaussians"""
115 d1 = gmm2map(to_draw,voxel_size,bounding_box,origin, fast)
116 print(
'will write GMM map to',out_fn)
121 """write density map directly from sklearn GMM (kinda slow) """
123 if not dmap_model
is None:
129 print(
'getting coords')
130 nvox=d1.get_number_of_voxels()
131 apos=[list(d1.get_location_by_voxel(nv))
for nv
in range(nvox)]
134 scores=gmm.score(apos)
137 for nv, score
in enumerate(scores):
138 d1.set_value(nv,exp(score))
139 print(
'will write GMM map to',out_fn)
142 def draw_points(pts,out_fn,trans=IMP.algebra.get_identity_transformation_3d(),
144 """ given some points (and optional transform), write them to chimera 'bild' format
145 colors flag only applies to ellipses, otherwise it'll be weird"""
146 with open(out_fn,
'w')
as outf:
152 outf.write(
'.color 1 0 0\n.dotat %.2f %.2f %.2f\n'
153 %(pt[0],pt[1],pt[2]))
158 outf.write(
'.color 0 1 0\n')
159 colors=[
'0 1 0',
'0 0 1',
'0 1 1']
160 for nt,t
in enumerate(pts[start:]):
161 if use_colors
and nt%2==0:
162 outf.write(
'.color %s\n' % colors[nt/2])
164 outf.write(
'.dotat %.2f %.2f %.2f\n' %(pt[0],pt[1],pt[2]))
168 def create_gmm_for_bead(mdl,
171 sampled_points=100000,
173 print(
'fitting bead with',n_components,
'gaussians')
176 IMP.em.write_map(dmap,
'test_intermed.mrc')
177 pts=IMP.isd.sample_points_from_density(dmap,sampled_points)
184 return density_particles,dmap
186 def sample_and_fit_to_particles(model,
189 sampled_points=1000000,
193 covariance_type=
'full',
194 multiply_by_total_mass=
True,
198 if multiply_by_total_mass:
200 print(
'add_component_density: will multiply by mass',mass_multiplier)
203 print(
'add_component_density: sampling points')
209 pts=IMP.isd.sample_points_from_density(dmap,sampled_points)
212 print(
'add_component_density: fitting GMM to',len(pts),
'points')
214 n_components=num_components,
216 ps=density_particles,
218 covariance_type=covariance_type,
219 mass_multiplier=mass_multiplier)
221 if not output_txt
is None:
223 if not output_map
is None:
226 voxel_size=voxel_size,
229 return density_particles
236 covariance_type=
'full',
241 mass_multiplier=1.0):
242 """fit a GMM to some points. Will return the score and the Akaike score.
243 Akaike information criterion for the current model fit. It is a measure
244 of the relative quality of the GMM that takes into account the
245 parsimony and the goodness of the fit.
246 if no particles are provided, they will be created
248 points: list of coordinates (python)
249 n_components: number of gaussians to create
251 ps: list of particles to be decorated. if empty, will add
252 num_iter: number of EM iterations
253 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'spherical'
254 min_covar: assign a minimum value to covariance term. That is used to have more spherical
256 init_centers: initial coordinates of the GMM
257 force_radii: fix the radii (spheres only)
258 force_weight: fix the weights
259 mass_multiplier: multiply the weights of all the gaussians by this value
260 dirichlet: use the DGMM fitting (can reduce number of components, takes longer)
266 from sklearn.mixture
import GMM
268 from sklearn.mixture
import GaussianMixture
271 print(
'creating GMM with n_components',n_components,
'n_iter',num_iter,
'covar type',covariance_type)
275 points = np.array(points)
276 weights_init = precisions_init =
None
277 if force_radii != -1.0:
278 print(
'warning: radii can no longer be forced, but setting '
279 'initial values to ', force_radii)
280 precisions_init = np.array([[1./force_radii]*3
281 for i
in range(n_components)])
282 if force_weight != -1.0:
283 print(
'warning: weights can no longer be forced, but setting '
284 'initial values to ', force_weight)
285 weights_init = np.array([force_weight]*n_components)
287 gmm = GaussianMixture(n_components=n_components,
289 covariance_type=covariance_type,
290 weights_init=weights_init,
291 precisions_init=precisions_init,
292 means_init=
None if init_centers==[]
297 if force_radii==-1.0:
301 covariance_type=
'spherical'
302 print(
'forcing spherical with radii',force_radii)
304 if force_weight==-1.0:
308 print(
'forcing weights to be',force_weight)
310 gmm = GMM(n_components=n_components, n_iter=num_iter,
311 covariance_type=covariance_type, min_covar=min_covar,
312 params=params, init_params=init_params)
313 if force_weight!=-1.0:
314 gmm.weights_=np.array([force_weight]*n_components)
315 if force_radii!=-1.0:
316 gmm.covars_=np.array([[force_radii]*3
for i
in range(n_components)])
318 gmm.means_=init_centers
320 model=gmm.fit(points)
321 score=gmm.score(points)
322 akaikescore=model.aic(points)
327 covars = gmm.covariances_
330 for ng
in range(n_components):
333 covar=np.diag(covar).tolist()
336 center=list(gmm.means_[ng])
337 weight=mass_multiplier*gmm.weights_[ng]
345 return (score,akaikescore)
352 covariance_type=
'full',
353 mass_multiplier=1.0):
354 """fit a GMM to some points. Will return core::Gaussians.
355 if no particles are provided, they will be created
357 points: list of coordinates (python)
358 n_components: number of gaussians to create
360 ps: list of particles to be decorated. if empty, will add
361 num_iter: number of EM iterations
362 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'spherical'
363 init_centers: initial coordinates of the GMM
364 force_radii: fix the radii (spheres only)
365 force_weight: fix the weights
366 mass_multiplier: multiply the weights of all the gaussians by this value
372 from sklearn.mixture
import BayesianGaussianMixture
374 from sklearn.mixture
import DPGMM
378 print(
'using dirichlet prior')
380 gmm = BayesianGaussianMixture(
381 weight_concentration_prior_type=
'dirichlet_process',
382 n_components=n_components, max_iter=num_iter,
383 covariance_type=covariance_type)
385 gmm = DPGMM(n_components=n_components, n_iter=num_iter,
386 covariance_type=covariance_type)
397 gmm.precisions_ = gmm.precs_
398 for ng
in range(n_components):
399 invcovar=gmm.precisions_[ng]
400 covar=np.linalg.inv(invcovar)
402 covar=np.diag(covar).tolist()
405 center=list(gmm.means_[ng])
406 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.
Gaussian3D get_gaussian_from_covariance(const Eigen::Matrix3d &covariance, const Vector3D ¢er)
Return a Gaussian centered at the origin from a covariance matrix.
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...
Eigen::Matrix3d get_covariance(const Gaussian3D &g)
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.