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
23 ''' read the output from write_gmm_to_text, decorate as Gaussian and Mass'''
29 weight=float(fields[2])
30 center=list(map(float,fields[3].split()))
31 covar=np.array(list(map(float,fields[4].split()))).reshape((3,3))
38 rmax=sqrt(max(g.get_variances()))*radius_scale
40 if not transform
is None:
44 '''write a list of gaussians to text. must be decorated as Gaussian and Mass'''
45 print(
'will write GMM text to',out_fn)
47 outf.write(
'#|num|weight|mean|covariance matrix|\n')
48 for ng,g
in enumerate(ps):
52 mean=list(shape.get_center())
53 fm=[ng,weight]+mean+covar
56 outf.write(
'|{}|{}|{} {} {}|{} {} {} {} {} {} {} {} {}|\n'.format(*fm))
59 outf.write(
'|{0}|{1}|{2} {3} {4}|{5} {6} {7} {8} {9} {10} {11} {12} {13}|\n'.format(*fm))
63 '''write density map from GMM. input can be either particles or gaussians'''
67 ps=[g.get_particle()
for g
in to_draw]
69 print(
'ps must be Particles or Gaussians')
71 print(
'will write GMM map to',out_fn)
72 if bounding_box
is None:
94 '''write density map directly from sklearn GMM (kinda slow) '''
96 if not dmap_model
is None:
102 print(
'getting coords')
103 nvox=d1.get_number_of_voxels()
104 apos=[list(d1.get_location_by_voxel(nv))
for nv
in range(nvox)]
107 scores=gmm.score(apos)
110 for nv, score
in enumerate(scores):
111 d1.set_value(nv,exp(score))
112 print(
'will write GMM map to',out_fn)
115 def draw_points(pts,out_fn,trans=IMP.algebra.get_identity_transformation_3d(),
117 ''' given some points (and optional transform), write them to chimera 'bild' format
118 colors flag only applies to ellipses, otherwise it'll be weird'''
119 outf=open(out_fn,
'w')
125 outf.write(
'.color 1 0 0\n.dotat %.2f %.2f %.2f\n' %(pt[0],pt[1],pt[2]))
130 outf.write(
'.color 0 1 0\n')
131 colors=[
'0 1 0',
'0 0 1',
'0 1 1']
132 for nt,t
in enumerate(pts[start:]):
133 if use_colors
and nt%2==0:
134 outf.write(
'.color %s\n' % colors[nt/2])
136 outf.write(
'.dotat %.2f %.2f %.2f\n' %(pt[0],pt[1],pt[2]))
141 def create_gmm_for_bead(mdl,
144 sampled_points=100000,
146 print(
'fitting bead with',n_components,
'gaussians')
150 pts=IMP.isd.sample_points_from_density(dmap,sampled_points)
157 return density_particles,dmap
159 def sample_and_fit_to_particles(model,
162 sampled_points=1000000,
166 covariance_type=
'full',
167 multiply_by_total_mass=
True,
171 if multiply_by_total_mass:
173 print(
'add_component_density: will multiply by mass',mass_multiplier)
176 print(
'add_component_density: sampling points')
182 pts=IMP.isd.sample_points_from_density(dmap,sampled_points)
185 print(
'add_component_density: fitting GMM to',len(pts),
'points')
187 n_components=num_components,
189 ps=density_particles,
191 covariance_type=covariance_type,
192 mass_multiplier=mass_multiplier)
194 if not output_txt
is None:
196 if not output_map
is None:
199 voxel_size=voxel_size,
202 return density_particles
209 covariance_type=
'full',
213 mass_multiplier=1.0):
214 '''fit a GMM to some points. Will return core::Gaussians.
215 if no particles are provided, they will be created
217 points: list of coordinates (python)
218 n_components: number of gaussians to create
220 ps: list of particles to be decorated. if empty, will add
221 num_iter: number of EM iterations
222 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'sphereical'
223 init_centers: initial coordinates of the GMM
224 force_radii: fix the radii (spheres only)
225 force_weight: fix the weights
226 mass_multiplier: multiply the weights of all the gaussians by this value
227 dirichlet: use the DGMM fitting (can reduce number of components, takes longer)
231 import sklearn.mixture
235 if force_radii==-1.0:
239 covariance_type=
'spherical'
240 print(
'forcing spherical with radii',force_radii)
242 if force_weight==-1.0:
246 print(
'forcing weights to be',force_weight)
248 print(
'creating GMM with params',params,
'init params',init_params,
'n_components',n_components,
'n_iter',num_iter,
'covar type',covariance_type)
249 gmm=sklearn.mixture.GMM(n_components=n_components,
251 covariance_type=covariance_type,
253 init_params=init_params)
255 if force_weight!=-1.0:
256 gmm.weights_=np.array([force_weight]*n_components)
257 if force_radii!=-1.0:
258 gmm.covars_=np.array([[force_radii]*3
for i
in range(n_components)])
260 gmm.means_=init_centers
264 print(
'>>> GMM score',gmm.score(points))
267 for ng
in range(n_components):
268 covar=gmm.covars_[ng]
270 covar=np.diag(covar).tolist()
273 center=list(gmm.means_[ng])
274 weight=mass_multiplier*gmm.weights_[ng]
287 covariance_type=
'full',
288 mass_multiplier=1.0):
289 '''fit a GMM to some points. Will return core::Gaussians.
290 if no particles are provided, they will be created
292 points: list of coordinates (python)
293 n_components: number of gaussians to create
295 ps: list of particles to be decorated. if empty, will add
296 num_iter: number of EM iterations
297 covariance_type: covar type for the gaussians. options: 'full', 'diagonal', 'sphereical'
298 init_centers: initial coordinates of the GMM
299 force_radii: fix the radii (spheres only)
300 force_weight: fix the weights
301 mass_multiplier: multiply the weights of all the gaussians by this value
305 import sklearn.mixture
308 print(
'using dirichlet prior')
309 gmm=sklearn.mixture.DPGMM(n_components=n_components,
311 covariance_type=covariance_type)
315 print(
'>>> GMM score',gmm.score(points))
321 for ng
in range(n_components):
322 invcovar=gmm.precs_[ng]
323 covar=np.linalg.inv(invcovar)
325 covar=np.diag(covar).tolist()
328 center=list(gmm.means_[ng])
329 weight=mass_multiplier*gmm.weights_[ng]
static Gaussian setup_particle(Model *m, ParticleIndex pi)
static FloatKey get_mass_key()
DensityMap * create_density_map(const algebra::GridD< 3, algebra::DenseGridStorageD< 3, float >, float > &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.
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)
void write_map(DensityMap *m, std::string filename)
Write a density map to a file.
Class for sampling a density map from particles.
The standard decorator for manipulating molecular structures.
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 hierarchy.
A decorator for a rigid body.
Functionality for loading, creating, manipulating and scoring atomic structures.