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))
81 """write density map from GMM. input can be either particles or gaussians"""
85 ps=[g.get_particle()
for g
in to_draw]
87 print(
'ps must be Particles or Gaussians')
89 print(
'will write GMM map to',out_fn)
90 if bounding_box
is None:
105 print(
'creating map')
108 if origin
is not None:
109 d1.set_origin(origin)
114 """write density map directly from sklearn GMM (kinda slow) """
116 if not dmap_model
is None:
122 print(
'getting coords')
123 nvox=d1.get_number_of_voxels()
124 apos=[list(d1.get_location_by_voxel(nv))
for nv
in range(nvox)]
127 scores=gmm.score(apos)
130 for nv, score
in enumerate(scores):
131 d1.set_value(nv,exp(score))
132 print(
'will write GMM map to',out_fn)
135 def draw_points(pts,out_fn,trans=IMP.algebra.get_identity_transformation_3d(),
137 """ given some points (and optional transform), write them to chimera 'bild' format
138 colors flag only applies to ellipses, otherwise it'll be weird"""
139 with open(out_fn,
'w')
as outf:
145 outf.write(
'.color 1 0 0\n.dotat %.2f %.2f %.2f\n'
146 %(pt[0],pt[1],pt[2]))
151 outf.write(
'.color 0 1 0\n')
152 colors=[
'0 1 0',
'0 0 1',
'0 1 1']
153 for nt,t
in enumerate(pts[start:]):
154 if use_colors
and nt%2==0:
155 outf.write(
'.color %s\n' % colors[nt/2])
157 outf.write(
'.dotat %.2f %.2f %.2f\n' %(pt[0],pt[1],pt[2]))
161 def create_gmm_for_bead(mdl,
164 sampled_points=100000,
166 print(
'fitting bead with',n_components,
'gaussians')
169 IMP.em.write_map(dmap,
'test_intermed.mrc')
170 pts=IMP.isd.sample_points_from_density(dmap,sampled_points)
177 return density_particles,dmap
179 def sample_and_fit_to_particles(model,
182 sampled_points=1000000,
186 covariance_type=
'full',
187 multiply_by_total_mass=
True,
191 if multiply_by_total_mass:
193 print(
'add_component_density: will multiply by mass',mass_multiplier)
196 print(
'add_component_density: sampling points')
202 pts=IMP.isd.sample_points_from_density(dmap,sampled_points)
205 print(
'add_component_density: fitting GMM to',len(pts),
'points')
207 n_components=num_components,
209 ps=density_particles,
211 covariance_type=covariance_type,
212 mass_multiplier=mass_multiplier)
214 if not output_txt
is None:
216 if not output_map
is None:
219 voxel_size=voxel_size,
222 return density_particles
229 covariance_type=
'full',
234 mass_multiplier=1.0):
235 """fit a GMM to some points. Will return the score and the Akaike score.
236 Akaike information criterion for the current model fit. It is a measure
237 of the relative quality of the GMM that takes into account the
238 parsimony and the goodness of the fit.
239 if no particles are provided, they will be created
241 points: list of coordinates (python)
242 n_components: number of gaussians to create
244 ps: list of particles to be decorated. if empty, will add
245 num_iter: number of EM iterations
246 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'spherical'
247 min_covar: assign a minimum value to covariance term. That is used to have more spherical
249 init_centers: initial coordinates of the GMM
250 force_radii: fix the radii (spheres only)
251 force_weight: fix the weights
252 mass_multiplier: multiply the weights of all the gaussians by this value
253 dirichlet: use the DGMM fitting (can reduce number of components, takes longer)
257 import sklearn.mixture
261 if force_radii==-1.0:
265 covariance_type=
'spherical'
266 print(
'forcing spherical with radii',force_radii)
268 if force_weight==-1.0:
272 print(
'forcing weights to be',force_weight)
274 print(
'creating GMM with params',params,
'init params',init_params,
'n_components',n_components,
'n_iter',num_iter,
'covar type',covariance_type)
275 gmm=sklearn.mixture.GMM(n_components=n_components,
277 covariance_type=covariance_type,
280 init_params=init_params)
282 if force_weight!=-1.0:
283 gmm.weights_=np.array([force_weight]*n_components)
284 if force_radii!=-1.0:
285 gmm.covars_=np.array([[force_radii]*3
for i
in range(n_components)])
287 gmm.means_=init_centers
289 model=gmm.fit(points)
290 score=gmm.score(points)
291 akaikescore=model.aic(points)
295 for ng
in range(n_components):
296 covar=gmm.covars_[ng]
298 covar=np.diag(covar).tolist()
301 center=list(gmm.means_[ng])
302 weight=mass_multiplier*gmm.weights_[ng]
310 return (score,akaikescore)
317 covariance_type=
'full',
318 mass_multiplier=1.0):
319 """fit a GMM to some points. Will return core::Gaussians.
320 if no particles are provided, they will be created
322 points: list of coordinates (python)
323 n_components: number of gaussians to create
325 ps: list of particles to be decorated. if empty, will add
326 num_iter: number of EM iterations
327 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'spherical'
328 init_centers: initial coordinates of the GMM
329 force_radii: fix the radii (spheres only)
330 force_weight: fix the weights
331 mass_multiplier: multiply the weights of all the gaussians by this value
335 import sklearn.mixture
338 print(
'using dirichlet prior')
339 gmm=sklearn.mixture.DPGMM(n_components=n_components,
341 covariance_type=covariance_type)
351 for ng
in range(n_components):
352 invcovar=gmm.precs_[ng]
353 covar=np.linalg.inv(invcovar)
355 covar=np.diag(covar).tolist()
358 center=list(gmm.means_[ng])
359 weight=mass_multiplier*gmm.weights_[ng]
static Gaussian setup_particle(Model *m, ParticleIndex pi)
static FloatKey get_mass_key()
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.
DenseGrid3D< float > get_rasterized(const Gaussian3Ds &gmm, const Floats &weights, double cell_width, const BoundingBox3D &bb)
Rasterize the Gaussians to a grid.
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 model particles.
Sphere3D get_enclosing_sphere(const Vector3Ds &ss)
Return a sphere containing the listed vectors.
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.