IMP logo
IMP Reference Guide  2.12.0
The Integrative Modeling Platform
restraints/saxs.py
1 """@namespace IMP.pmi.restraints.saxs
2 Restraints for handling small angle x-ray (SAXS) data.
3 """
4 
5 from __future__ import print_function
6 import math
7 import IMP
8 import IMP.core
9 import IMP.algebra
10 import IMP.atom
11 import IMP.container
12 import IMP.pmi.tools
13 import IMP.pmi.restraints
14 import IMP.saxs
15 import warnings
16 
18 
19  """Basic SAXS restraint."""
20 
21  def __init__(self, input_objects, saxs_datafile, weight=1.0,
22  ff_type=IMP.saxs.HEAVY_ATOMS, label=None, maxq="standard"):
23  """Builds the restraint.
24  @param input_objects A list of hierarchies or PMI objects that the
25  SAXS restraint will be applied to. This hierarchy MUST be
26  atomic. You can pass a list of CA atom particles to evaluate
27  at residue resolution
28  @param saxs_datafile the SAXS .dat file.
29  @param weight Restraint score coefficient
30  @param ff_type the form factor to use, of the following types:
31  - IMP.saxs.HEAVY_ATOMS: use form factors with implicit
32  hydrogens
33  - IMP.saxs.ALL_ATOMS: use individual form factors for all
34  atoms. Does not build missing hydrogens.
35  - IMP.saxs.CA_ATOMS: use residue based form factors
36  centered at CA atoms
37  @param label Label for the restraint in outputs
38 
39  @param maxq - maximum q value that the restraint will be evaluated at
40  Default vaules for ff_type = ALL_ATOMS : 0.5. HEAVY_ATOMS : 0.4,
41  CA_ATOMS and RESIDUES = 0.15. These values were eyeballed
42  by comparing ALL_ATOM calculated SAXS profiles to those calculated
43  with the reduced representations, so could be improved.
44  """
45  # Get all hierarchies.
46  hiers = IMP.pmi.tools.input_adaptor(input_objects,
47  flatten=True)
48  model = list(hiers)[0].get_model()
49  super(SAXSRestraint, self).__init__(model, label=label, weight=weight)
50 
51  # Determine maxq to compare computed and experimental profiles
52  if maxq == "standard":
53  if ff_type == IMP.saxs.CA_ATOMS or ff_type == IMP.saxs.RESIDUES:
54  maxq = 0.15
55  elif ff_type == IMP.saxs.HEAVY_ATOMS:
56  maxq = 0.4
57  else:
58  maxq = 0.5
59  elif type(maxq) == float:
60  if maxq < 0.01 or maxq > 4.0:
61  raise Exception("SAXSRestraint: maxq must be set between 0.01 and 4.0")
62  if (ff_type == IMP.saxs.CA_ATOMS or ff_type == IMP.saxs.RESIDUES) and maxq > 0.15:
63  warnings.warn("SAXSRestraint: for residue-resolved form "
64  "factors, a maxq > 0.15 is not recommended!",
66  else:
67  raise Exception("SAXSRestraint: maxq must be set to a number between 0.01 and 4.0")
68 
69  # Depending on the type of FF used, get the correct particles
70  # from the hierarchies list and create an IMP::saxs::Profile()
71  # at the appropriate maxq.
72  if ff_type == IMP.saxs.RESIDUES:
73  self.particles = IMP.atom.Selection(
74  hiers, resolution=1).get_selected_particles()
75  self.profile = IMP.saxs.Profile(saxs_datafile, False, maxq)
76  elif ff_type == IMP.saxs.CA_ATOMS:
77  self.particles = IMP.atom.Selection(
78  hiers, atom_type=IMP.atom.AT_CA).get_selected_particles()
79  self.profile = IMP.saxs.Profile(saxs_datafile, False, maxq)
80  elif ff_type == IMP.saxs.HEAVY_ATOMS or ff_type == IMP.saxs.ALL_ATOMS:
81  self.particles = IMP.atom.Selection(
82  hiers, resolution=0).get_selected_particles()
83  self.profile = IMP.saxs.Profile(saxs_datafile, False, maxq)
84  else:
85  raise Exception("SAXSRestraint: Must provide an IMP.saxs atom type: RESIDUES, CA_ATOMS, HEAVY_ATOMS or ALL_ATOMS")
86 
87  if len(self.particles) == 0:
88  raise Exception("SAXSRestraint: There are no selected particles")
89 
90  self.restraint = IMP.saxs.Restraint(self.particles, self.profile,
91  ff_type)
92  self.rs.add_restraint(self.restraint)
93 
94 
96 
97  """Basic SAXS restraint using ISD."""
98 
99  import IMP.isd
100  try:
101  import IMP.isd2
102  except:
103  print("Module isd2 not installed. Cannot use SAXSISDRestraint")
104 
105  def __init__(self, representation, profile, resolution=0, weight=1,
106  ff_type=IMP.saxs.HEAVY_ATOMS, label=None):
107 
108  model = representation.prot.get_model()
109  super(SAXSISDRestraint, self).__init__(model, label=label,
110  weight=weight)
111 
112  self.taumaxtrans = 0.05
113  self.prof = IMP.saxs.Profile(profile)
114 
115  self.atoms = IMP.pmi.tools.select(
116  representation,
117  resolution=resolution)
118 
119  # gamma nuisance
120  self.gamma = IMP.pmi.tools.SetupNuisance(
121  self.model, 1., 0., None, False).get_particle()
122 
123  # sigma nuisance
124  self.sigma = IMP.pmi.tools.SetupNuisance(self.model, 10.0, 0., None, False
125  ).get_particle()
126 
127  # tau nuisance, optimized
128  self.tau = IMP.pmi.tools.SetupNuisance(self.model, 1., 0., None, False,
129  ).get_particle()
130 
131  # c1 and c2, optimized
132  self.c1 = IMP.pmi.tools.SetupNuisance(self.model, 1.0, 0.95, 1.05,
133  True).get_particle()
134  self.c2 = IMP.pmi.tools.SetupNuisance(self.model, 0.0, -2., 4.,
135  True).get_particle()
136 
137  # weight, optimized
138  self.w = IMP.pmi.tools.SetupWeight(self.model).get_particle()
139  IMP.isd.Weight(self.w).set_weights_are_optimized(False)
140 
141  # take identity covariance matrix for the start
142  self.cov = [[1 if i == j else 0 for j in range(self.prof.size())]
143  for i in range(self.prof.size())]
144 
145  print("create saxs restraint")
146  self.saxs = IMP.isd2.SAXSRestraint(self.prof, self.sigma, self.tau,
147  self.gamma, self.w, self.c1,
148  self.c2)
149  self.saxs.add_scatterer(self.atoms, self.cov, ff_type)
150 
151  self.rs.add_restraint(self.saxs)
152 
153  # self.saxs_stuff={'nuis':(sigma,gamma),'cov':cov,
154  # 'exp':prof,'th':tmp}
155 
156  self.rs2 = self._create_restraint_set('Prior')
157  # jeffreys restraints for nuisances
158  j1 = IMP.isd.JeffreysRestraint(self.model, self.sigma)
159  self.rs2.add_restraint(j1)
160  j2 = IMP.isd.JeffreysRestraint(self.model, self.tau)
161  self.rs2.add_restraint(j2)
162  j3 = IMP.isd.JeffreysRestraint(self.model, self.gamma)
163  self.rs2.add_restraint(j3)
164 
165  pw = IMP.isd.Weight(self.w)
166  pw.set_weights(pw.get_unit_simplex().get_barycenter())
167  pw.set_weights_are_optimized(True)
168 
169  def optimize_sigma(self):
170  """Set sigma to the value that maximizes its conditional likelihood"""
171  self.model.update()
172  sigma2hat = self.saxs.get_sigmasq_scale_parameter() \
173  / (self.saxs.get_sigmasq_shape_parameter() + 1)
174  IMP.isd.Scale(self.sigma).set_scale(math.sqrt(sigma2hat))
175 
176  def optimize_gamma(self):
177  """Set gamma to the value that maximizes its conditional likelihood"""
178  self.model.update()
179  gammahat = math.exp(self.saxs.get_loggamma_variance_parameter() *
180  self.saxs.get_loggamma_jOg_parameter())
181  IMP.isd.Scale(self.gamma).set_scale(gammahat)
182 
183  def optimize_tau(self, ltaumin=-2, ltaumax=3, npoints=100):
184  values = []
185  self.model.update()
186  IMP.atom.write_pdb(self.atoms, 'tauvals.pdb')
187  fl = open('tauvals.txt', 'w')
188  for tauval in self._logspace(ltaumin, ltaumax, npoints):
189  IMP.isd.Scale(self.tau).set_scale(tauval)
190  try:
191  values.append((self.model.evaluate(False), tauval))
192  except:
193  pass
194  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
195  values.sort()
196  ltcenter = math.log(values[0][1]) / math.log(10)
197  spacing = (ltaumax - ltaumin) / float(npoints)
198  values = []
199  for tauval in self._logspace(
200  ltcenter - 2 * spacing, ltcenter + 2 * spacing,
201  npoints):
202  IMP.isd.Scale(self.tau).set_scale(tauval)
203  values.append((self.model.evaluate(False), tauval))
204  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
205  values.sort()
206  IMP.isd.Scale(self.tau).set_scale(values[0][1])
207 
208  def get_gamma_value(self):
209  """Get value of gamma."""
210  return self.gamma.get_scale()
211 
212  def set_taumaxtrans(self, taumaxtrans):
213  self.taumaxtrans = taumaxtrans
214 
215  def draw_sigma(self):
216  """Draw 1/sigma2 from gamma distribution."""
217  self.model.update()
218  self.saxs.draw_sigma()
219 
220  def draw_gamma(self):
221  """Draw gamma from lognormal distribution."""
222  self.model.update()
223  self.saxs.draw_gamma()
224 
225  def update_covariance_matrix(self):
226  c1 = IMP.isd.Nuisance(self.c1).get_nuisance()
227  c2 = IMP.isd.Nuisance(self.c2).get_nuisance()
228  # tau = IMP.isd.Nuisance(self.tau).get_nuisance()
229  tau = 1.0
230  self.cov = IMP.isd2.compute_relative_covariance(self.atoms, c1, c2,
231  tau, self.prof)
232  # for i in xrange(len(self.cov)):
233  # for j in xrange(len(self.cov)):
234  # self.cov[i][j] = self.cov[i][j]/tau**2
235  self.saxs.set_cov(0, self.cov)
236 
237  def write_covariance_matrix(self, fname):
238  fl = open(fname, 'w')
239  for line in self.cov:
240  for i in line:
241  fl.write('%G ' % i)
242  fl.write('\n')
243 
244  def get_output(self):
245  output = super(SAXSISDRestraint, self).get_output()
246  suffix = self._get_label_suffix()
247  output["SAXSISDRestraint_Sigma" +
248  suffix] = str(self.sigma.get_scale())
249  output["SAXSISDRestraint_Tau" + suffix] = str(self.tau.get_scale())
250  output["SAXSISDRestraint_Gamma" +
251  suffix] = str(self.gamma.get_scale())
252  return output
253 
254  @staticmethod
255  def _logspace(a, b, num=100):
256  """Mimick numpy's logspace function"""
257  for i in range(num):
258  val = a + float(b - a) / float(num - 1) * i
259  yield 10 ** val
Add weights to a particle.
Definition: Weight.h:29
Various classes to hold sets of particles.
def optimize_gamma
Set gamma to the value that maximizes its conditional likelihood.
Calculate score based on fit to SAXS profile.
Miscellaneous utilities.
Definition: tools.py:1
void write_pdb(const Selection &mhd, TextOutput out, unsigned int model=1)
Add scale parameter to particle.
Definition: Scale.h:24
Classes to handle different kinds of restraints.
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: tools.py:918
def optimize_sigma
Set sigma to the value that maximizes its conditional likelihood.
def draw_gamma
Draw gamma from lognormal distribution.
Add nuisance parameter to particle.
Definition: Nuisance.h:25
Basic SAXS restraint using ISD.
def __init__
Builds the restraint.
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...
The general base class for IMP exceptions.
Definition: exception.h:49
def draw_sigma
Draw 1/sigma2 from gamma distribution.
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
def get_output
Get outputs to write to stat files.
Support for small angle X-ray scattering (SAXS) data.
Warning for probably incorrect input parameters.
Base class for PMI restraints, which wrap IMP.Restraint(s).
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...