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)
274 points = np.array(points)
275 weights_init = precisions_init =
None
276 if force_radii != -1.0:
277 print(
'warning: radii can no longer be forced, but setting '
278 'initial values to ', force_radii)
279 precisions_init = np.array([[1./force_radii]*3
280 for i
in range(n_components)])
281 if force_weight != -1.0:
282 print(
'warning: weights can no longer be forced, but setting '
283 'initial values to ', force_weight)
284 weights_init = np.array([force_weight]*n_components)
286 gmm = GaussianMixture(n_components=n_components,
288 covariance_type=covariance_type,
289 weights_init=weights_init,
290 precisions_init=precisions_init,
291 means_init=
None if init_centers==[]
296 if force_radii==-1.0:
300 covariance_type=
'spherical'
301 print(
'forcing spherical with radii',force_radii)
303 if force_weight==-1.0:
307 print(
'forcing weights to be',force_weight)
309 gmm = GMM(n_components=n_components, n_iter=num_iter,
310 covariance_type=covariance_type, min_covar=min_covar,
311 params=params, init_params=init_params)
312 if force_weight!=-1.0:
313 gmm.weights_=np.array([force_weight]*n_components)
314 if force_radii!=-1.0:
315 gmm.covars_=np.array([[force_radii]*3
for i
in range(n_components)])
317 gmm.means_=init_centers
319 model=gmm.fit(points)
320 score=gmm.score(points)
321 akaikescore=model.aic(points)
326 covars = gmm.covariances_
329 for ng
in range(n_components):
332 covar=np.diag(covar).tolist()
335 center=list(gmm.means_[ng])
336 weight=mass_multiplier*gmm.weights_[ng]
344 return (score,akaikescore)
351 covariance_type=
'full',
352 mass_multiplier=1.0):
353 """fit a GMM to some points. Will return core::Gaussians.
354 if no particles are provided, they will be created
356 points: list of coordinates (python)
357 n_components: number of gaussians to create
359 ps: list of particles to be decorated. if empty, will add
360 num_iter: number of EM iterations
361 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'spherical'
362 init_centers: initial coordinates of the GMM
363 force_radii: fix the radii (spheres only)
364 force_weight: fix the weights
365 mass_multiplier: multiply the weights of all the gaussians by this value
369 import sklearn.mixture
372 print(
'using dirichlet prior')
373 gmm=sklearn.mixture.DPGMM(n_components=n_components,
375 covariance_type=covariance_type)
385 for ng
in range(n_components):
386 invcovar=gmm.precs_[ng]
387 covar=np.linalg.inv(invcovar)
389 covar=np.diag(covar).tolist()
392 center=list(gmm.means_[ng])
393 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.