IMP logo
IMP Reference Guide  2.11.1
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(True)
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  def optimize_sigma(self):
166  """Set sigma to the value that maximizes its conditional likelihood"""
167  self.model.update()
168  sigma2hat = self.saxs.get_sigmasq_scale_parameter() \
169  / (self.saxs.get_sigmasq_shape_parameter() + 1)
170  IMP.isd.Scale(self.sigma).set_scale(math.sqrt(sigma2hat))
171 
172  def optimize_gamma(self):
173  """Set gamma to the value that maximizes its conditional likelihood"""
174  self.model.update()
175  gammahat = math.exp(self.saxs.get_loggamma_variance_parameter() *
176  self.saxs.get_loggamma_jOg_parameter())
177  IMP.isd.Scale(self.gamma).set_scale(gammahat)
178 
179  def optimize_tau(self, ltaumin=-2, ltaumax=3, npoints=100):
180  values = []
181  self.model.update()
182  IMP.atom.write_pdb(self.atoms, 'tauvals.pdb')
183  fl = open('tauvals.txt', 'w')
184  for tauval in self._logspace(ltaumin, ltaumax, npoints):
185  IMP.isd.Scale(self.tau).set_scale(tauval)
186  try:
187  values.append((self.model.evaluate(False), tauval))
188  except:
189  pass
190  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
191  values.sort()
192  ltcenter = math.log(values[0][1]) / math.log(10)
193  spacing = (ltaumax - ltaumin) / float(npoints)
194  values = []
195  for tauval in self._logspace(
196  ltcenter - 2 * spacing, ltcenter + 2 * spacing,
197  npoints):
198  IMP.isd.Scale(self.tau).set_scale(tauval)
199  values.append((self.model.evaluate(False), tauval))
200  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
201  values.sort()
202  IMP.isd.Scale(self.tau).set_scale(values[0][1])
203 
204  def get_gamma_value(self):
205  """Get value of gamma."""
206  return self.gamma.get_scale()
207 
208  def set_taumaxtrans(self, taumaxtrans):
209  self.taumaxtrans = taumaxtrans
210 
211  def draw_sigma(self):
212  """Draw 1/sigma2 from gamma distribution."""
213  self.model.update()
214  self.saxs.draw_sigma()
215 
216  def draw_gamma(self):
217  """Draw gamma from lognormal distribution."""
218  self.model.update()
219  self.saxs.draw_gamma()
220 
221  def update_covariance_matrix(self):
222  c1 = IMP.isd.Nuisance(self.c1).get_nuisance()
223  c2 = IMP.isd.Nuisance(self.c2).get_nuisance()
224  # tau = IMP.isd.Nuisance(self.tau).get_nuisance()
225  tau = 1.0
226  self.cov = IMP.isd2.compute_relative_covariance(self.atoms, c1, c2,
227  tau, self.prof)
228  # for i in xrange(len(self.cov)):
229  # for j in xrange(len(self.cov)):
230  # self.cov[i][j] = self.cov[i][j]/tau**2
231  self.saxs.set_cov(0, self.cov)
232 
233  def write_covariance_matrix(self, fname):
234  fl = open(fname, 'w')
235  for line in self.cov:
236  for i in line:
237  fl.write('%G ' % i)
238  fl.write('\n')
239 
240  def get_output(self):
241  output = super(SAXSISDRestraint, self).get_output()
242  suffix = self._get_label_suffix()
243  output["SAXSISDRestraint_Sigma" +
244  suffix] = str(self.sigma.get_scale())
245  output["SAXSISDRestraint_Tau" + suffix] = str(self.tau.get_scale())
246  output["SAXSISDRestraint_Gamma" +
247  suffix] = str(self.gamma.get_scale())
248  return output
249 
250  @staticmethod
251  def _logspace(a, b, num=100):
252  """Mimick numpy's logspace function"""
253  for i in range(num):
254  val = a + float(b - a) / float(num - 1) * i
255  yield 10 ** val
Add weights for a set of states to a particle.
Definition: Weight.h:24
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:1275
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.
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:545
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 ...