IMP logo
IMP Reference Guide  develop.8b7eb21230,2024/06/24
The Integrative Modeling Platform
/restraints/saxs.py
1 """@namespace IMP.pmi1.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.pmi1.tools
14 import IMP.saxs
15 
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 values 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.pmi1.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  print("SAXSRestraint: WARNING> for residue-resolved form factors, a maxq > 0.15 is not recommended!")
64  else:
65  raise Exception("SAXSRestraint: maxq must be set to a number between 0.01 and 4.0")
66 
67  # Depending on the type of FF used, get the correct particles
68  # from the hierarchies list and create an IMP::saxs::Profile()
69  # at the appropriate maxq.
70  if ff_type == IMP.saxs.RESIDUES:
71  self.particles = IMP.atom.Selection(
72  hiers, resolution=1).get_selected_particles()
73  self.profile = IMP.saxs.Profile(saxs_datafile, False, maxq)
74  elif ff_type == IMP.saxs.CA_ATOMS:
75  self.particles = IMP.atom.Selection(
76  hiers, atom_type=IMP.atom.AT_CA).get_selected_particles()
77  self.profile = IMP.saxs.Profile(saxs_datafile, False, maxq)
78  elif ff_type == IMP.saxs.HEAVY_ATOMS or ff_type == IMP.saxs.ALL_ATOMS:
79  self.particles = IMP.atom.Selection(
80  hiers, resolution=0).get_selected_particles()
81  self.profile = IMP.saxs.Profile(saxs_datafile, False, maxq)
82  else:
83  raise Exception("SAXSRestraint: Must provide an IMP.saxs atom type: RESIDUES, CA_ATOMS, HEAVY_ATOMS or ALL_ATOMS")
84 
85  if len(self.particles) == 0:
86  raise Exception("SAXSRestraint: There are no selected particles")
87 
88  self.restraint = IMP.saxs.Restraint(self.particles, self.profile,
89  ff_type)
90  self.rs.add_restraint(self.restraint)
91 
92 
94 
95  """Basic SAXS restraint using ISD."""
96 
97  import IMP.isd
98  try:
99  import IMP.isd2
100  except:
101  print("Module isd2 not installed. Cannot use SAXSISDRestraint")
102 
103  def __init__(self, representation, profile, resolution=0, weight=1,
104  ff_type=IMP.saxs.HEAVY_ATOMS, label=None):
105 
106  model = representation.prot.get_model()
107  super(SAXSISDRestraint, self).__init__(model, label=label,
108  weight=weight)
109 
110  self.taumaxtrans = 0.05
111  self.prof = IMP.saxs.Profile(profile)
112 
113  self.atoms = IMP.pmi1.tools.select(
114  representation,
115  resolution=resolution)
116 
117  # gamma nuisance
118  self.gamma = IMP.pmi1.tools.SetupNuisance(
119  self.model, 1., 0., None, False).get_particle()
120 
121  # sigma nuisance
122  self.sigma = IMP.pmi1.tools.SetupNuisance(self.model, 10.0, 0., None, False
123  ).get_particle()
124 
125  # tau nuisance, optimized
126  self.tau = IMP.pmi1.tools.SetupNuisance(self.model, 1., 0., None, False,
127  ).get_particle()
128 
129  # c1 and c2, optimized
130  self.c1 = IMP.pmi1.tools.SetupNuisance(self.model, 1.0, 0.95, 1.05,
131  True).get_particle()
132  self.c2 = IMP.pmi1.tools.SetupNuisance(self.model, 0.0, -2., 4.,
133  True).get_particle()
134 
135  # weight, optimized
136  self.w = IMP.pmi1.tools.SetupWeight(self.model).get_particle()
137  IMP.isd.Weight(self.w).set_weights_are_optimized(False)
138 
139  # take identity covariance matrix for the start
140  self.cov = [[1 if i == j else 0 for j in range(self.prof.size())]
141  for i in range(self.prof.size())]
142 
143  print("create saxs restraint")
144  self.saxs = IMP.isd2.SAXSRestraint(self.prof, self.sigma, self.tau,
145  self.gamma, self.w, self.c1,
146  self.c2)
147  self.saxs.add_scatterer(self.atoms, self.cov, ff_type)
148 
149  self.rs.add_restraint(self.saxs)
150 
151  # self.saxs_stuff={'nuis':(sigma,gamma),'cov':cov,
152  # 'exp':prof,'th':tmp}
153 
154  self.rs2 = self._create_restraint_set('Prior')
155  # jeffreys restraints for nuisances
156  j1 = IMP.isd.JeffreysRestraint(self.model, self.sigma)
157  self.rs2.add_restraint(j1)
158  j2 = IMP.isd.JeffreysRestraint(self.model, self.tau)
159  self.rs2.add_restraint(j2)
160  j3 = IMP.isd.JeffreysRestraint(self.model, self.gamma)
161  self.rs2.add_restraint(j3)
162 
163  pw = IMP.isd.Weight(self.w)
164  pw.set_weights(pw.get_unit_simplex().get_barycenter())
165  pw.set_weights_are_optimized(True)
166 
167  def optimize_sigma(self):
168  """Set sigma to the value that maximizes its conditional likelihood"""
169  self.model.update()
170  sigma2hat = self.saxs.get_sigmasq_scale_parameter() \
171  / (self.saxs.get_sigmasq_shape_parameter() + 1)
172  IMP.isd.Scale(self.sigma).set_scale(math.sqrt(sigma2hat))
173 
174  def optimize_gamma(self):
175  """Set gamma to the value that maximizes its conditional likelihood"""
176  self.model.update()
177  gammahat = math.exp(self.saxs.get_loggamma_variance_parameter() *
178  self.saxs.get_loggamma_jOg_parameter())
179  IMP.isd.Scale(self.gamma).set_scale(gammahat)
180 
181  def optimize_tau(self, ltaumin=-2, ltaumax=3, npoints=100):
182  values = []
183  self.model.update()
184  IMP.atom.write_pdb(self.atoms, 'tauvals.pdb')
185  fl = open('tauvals.txt', 'w')
186  for tauval in self._logspace(ltaumin, ltaumax, npoints):
187  IMP.isd.Scale(self.tau).set_scale(tauval)
188  try:
189  values.append((self.model.evaluate(False), tauval))
190  except:
191  pass
192  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
193  values.sort()
194  ltcenter = math.log(values[0][1]) / math.log(10)
195  spacing = (ltaumax - ltaumin) / float(npoints)
196  values = []
197  for tauval in self._logspace(
198  ltcenter - 2 * spacing, ltcenter + 2 * spacing,
199  npoints):
200  IMP.isd.Scale(self.tau).set_scale(tauval)
201  values.append((self.model.evaluate(False), tauval))
202  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
203  values.sort()
204  IMP.isd.Scale(self.tau).set_scale(values[0][1])
205 
206  def get_gamma_value(self):
207  """Get value of gamma."""
208  return self.gamma.get_scale()
209 
210  def set_taumaxtrans(self, taumaxtrans):
211  self.taumaxtrans = taumaxtrans
212 
213  def draw_sigma(self):
214  """Draw 1/sigma2 from gamma distribution."""
215  self.model.update()
216  self.saxs.draw_sigma()
217 
218  def draw_gamma(self):
219  """Draw gamma from lognormal distribution."""
220  self.model.update()
221  self.saxs.draw_gamma()
222 
223  def update_covariance_matrix(self):
224  c1 = IMP.isd.Nuisance(self.c1).get_nuisance()
225  c2 = IMP.isd.Nuisance(self.c2).get_nuisance()
226  # tau = IMP.isd.Nuisance(self.tau).get_nuisance()
227  tau = 1.0
228  self.cov = IMP.isd2.compute_relative_covariance(self.atoms, c1, c2,
229  tau, self.prof)
230  # for i in xrange(len(self.cov)):
231  # for j in xrange(len(self.cov)):
232  # self.cov[i][j] = self.cov[i][j]/tau**2
233  self.saxs.set_cov(0, self.cov)
234 
235  def write_covariance_matrix(self, fname):
236  fl = open(fname, 'w')
237  for line in self.cov:
238  for i in line:
239  fl.write('%G ' % i)
240  fl.write('\n')
241 
242  def get_output(self):
243  output = super(SAXSISDRestraint, self).get_output()
244  suffix = self._get_label_suffix()
245  output["SAXSISDRestraint_Sigma" +
246  suffix] = str(self.sigma.get_scale())
247  output["SAXSISDRestraint_Tau" + suffix] = str(self.tau.get_scale())
248  output["SAXSISDRestraint_Gamma" +
249  suffix] = str(self.gamma.get_scale())
250  return output
251 
252  @staticmethod
253  def _logspace(a, b, num=100):
254  """Mimic numpy's logspace function"""
255  for i in range(num):
256  val = a + float(b - a) / float(num - 1) * i
257  yield 10 ** val
Add weights to a particle.
Definition: Weight.h:29
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: /tools.py:1621
def draw_gamma
Draw gamma from lognormal distribution.
def optimize_sigma
Set sigma to the value that maximizes its conditional likelihood.
Various classes to hold sets of particles.
Calculate score based on fit to SAXS profile.
void write_pdb(const Selection &mhd, TextOutput out, unsigned int model=1)
Classes to handle different kinds of restraints.
Add scale parameter to particle.
Definition: Scale.h:24
def optimize_gamma
Set gamma to the value that maximizes its conditional likelihood.
def get_output
Get outputs to write to stat files.
Miscellaneous utilities.
Definition: /tools.py:1
Add nuisance parameter to particle.
Definition: Nuisance.h:27
Base class for PMI restraints, which wrap IMP.Restraint(s).
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: /tools.py:727
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:48
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:70
Basic SAXS restraint using ISD.
Support for small angle X-ray scattering (SAXS) data.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...