IMP logo
IMP Reference Guide  2.8.0
The Integrative Modeling Platform
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 
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):
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  hiers = IMP.pmi.tools.input_adaptor(input_objects, pmi_resolution=0,
40  flatten=True)
41  m = list(hiers)[0].get_model()
42  super(SAXSRestraint, self).__init__(m, label=label, weight=weight)
43  self.profile = IMP.saxs.Profile(saxs_datafile)
44 
45  if ff_type == IMP.saxs.CA_ATOMS:
46  self.particles = IMP.atom.Selection(
47  hiers, atom_type=IMP.atom.AT_CA).get_selected_particles()
48  elif ff_type == IMP.saxs.HEAVY_ATOMS:
49  self.particles = IMP.atom.Selection(
50  hiers, resolution=0).get_selected_particles()
51  elif ff_type == IMP.saxs.ALL_ATOMS:
52  self.particles = IMP.atom.Selection(
53  hiers, resolution=0).get_selected_particles()
54  else:
55  raise Exception("SAXSRestraint: Must provide an IMP.saxs atom "
56  "type: CA_ATOMS, HEAVY_ATOMS or ALL_ATOMS")
57  if len(self.particles) == 0:
58  raise Exception("SAXSRestraint: There are no selected particles")
59 
60  self.restraint = IMP.saxs.Restraint(self.particles, self.profile,
61  ff_type)
62  self.rs.add_restraint(self.restraint)
63 
64 
66 
67  """Basic SAXS restraint using ISD."""
68 
69  import IMP.isd
70  try:
71  import IMP.isd2
72  except:
73  print("Module isd2 not installed. Cannot use SAXSISDRestraint")
74 
75  def __init__(self, representation, profile, resolution=0, weight=1,
76  ff_type=IMP.saxs.HEAVY_ATOMS, label=None):
77 
78  m = representation.prot.get_model()
79  super(SAXSISDRestraint, self).__init__(m, label=label, weight=weight)
80 
81  self.taumaxtrans = 0.05
82  self.prof = IMP.saxs.Profile(profile)
83 
84  self.atoms = IMP.pmi.tools.select(
85  representation,
86  resolution=resolution)
87 
88  # gamma nuisance
89  self.gamma = IMP.pmi.tools.SetupNuisance(
90  self.m, 1., 0., None, False).get_particle()
91 
92  # sigma nuisance
93  self.sigma = IMP.pmi.tools.SetupNuisance(self.m, 10.0, 0., None, False
94  ).get_particle()
95 
96  # tau nuisance, optimized
97  self.tau = IMP.pmi.tools.SetupNuisance(self.m, 1., 0., None, False,
98  ).get_particle()
99 
100  # c1 and c2, optimized
101  self.c1 = IMP.pmi.tools.SetupNuisance(self.m, 1.0, 0.95, 1.05,
102  True).get_particle()
103  self.c2 = IMP.pmi.tools.SetupNuisance(self.m, 0.0, -2., 4.,
104  True).get_particle()
105 
106  # weight, optimized
107  self.w = IMP.pmi.tools.SetupWeight(self.m).get_particle()
108  IMP.isd.Weight(self.w).set_weights_are_optimized(True)
109 
110  # take identity covariance matrix for the start
111  self.cov = [[1 if i == j else 0 for j in range(self.prof.size())]
112  for i in range(self.prof.size())]
113 
114  print("create saxs restraint")
115  self.saxs = IMP.isd2.SAXSRestraint(self.prof, self.sigma, self.tau,
116  self.gamma, self.w, self.c1,
117  self.c2)
118  self.saxs.add_scatterer(self.atoms, self.cov, ff_type)
119 
120  self.rs.add_restraint(self.saxs)
121 
122  # self.saxs_stuff={'nuis':(sigma,gamma),'cov':cov,
123  # 'exp':prof,'th':tmp}
124 
125  self.rs2 = self._create_restraint_set('Prior')
126  # jeffreys restraints for nuisances
127  j1 = IMP.isd.JeffreysRestraint(self.m, self.sigma)
128  self.rs2.add_restraint(j1)
129  j2 = IMP.isd.JeffreysRestraint(self.m, self.tau)
130  self.rs2.add_restraint(j2)
131  j3 = IMP.isd.JeffreysRestraint(self.m, self.gamma)
132  self.rs2.add_restraint(j3)
133 
134  def optimize_sigma(self):
135  """Set sigma to the value that maximizes its conditional likelihood"""
136  self.m.update()
137  sigma2hat = self.saxs.get_sigmasq_scale_parameter() \
138  / (self.saxs.get_sigmasq_shape_parameter() + 1)
139  IMP.isd.Scale(self.sigma).set_scale(math.sqrt(sigma2hat))
140 
141  def optimize_gamma(self):
142  """Set gamma to the value that maximizes its conditional likelihood"""
143  self.m.update()
144  gammahat = math.exp(self.saxs.get_loggamma_variance_parameter() *
145  self.saxs.get_loggamma_jOg_parameter())
146  IMP.isd.Scale(self.gamma).set_scale(gammahat)
147 
148  def optimize_tau(self, ltaumin=-2, ltaumax=3, npoints=100):
149  values = []
150  self.m.update()
151  IMP.atom.write_pdb(self.atoms, 'tauvals.pdb')
152  fl = open('tauvals.txt', 'w')
153  for tauval in self._logspace(ltaumin, ltaumax, npoints):
154  IMP.isd.Scale(self.tau).set_scale(tauval)
155  try:
156  values.append((self.m.evaluate(False), tauval))
157  except:
158  pass
159  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
160  values.sort()
161  ltcenter = math.log(values[0][1]) / math.log(10)
162  spacing = (ltaumax - ltaumin) / float(npoints)
163  values = []
164  for tauval in self._logspace(
165  ltcenter - 2 * spacing, ltcenter + 2 * spacing,
166  npoints):
167  IMP.isd.Scale(self.tau).set_scale(tauval)
168  values.append((self.m.evaluate(False), tauval))
169  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
170  values.sort()
171  IMP.isd.Scale(self.tau).set_scale(values[0][1])
172 
173  def get_gamma_value(self):
174  """Get value of gamma."""
175  return self.gamma.get_scale()
176 
177  def set_taumaxtrans(self, taumaxtrans):
178  self.taumaxtrans = taumaxtrans
179 
180  def draw_sigma(self):
181  """Draw 1/sigma2 from gamma distribution."""
182  self.m.update()
183  self.saxs.draw_sigma()
184 
185  def draw_gamma(self):
186  """Draw gamma from lognormal distribution."""
187  self.m.update()
188  self.saxs.draw_gamma()
189 
190  def update_covariance_matrix(self):
191  c1 = IMP.isd.Nuisance(self.c1).get_nuisance()
192  c2 = IMP.isd.Nuisance(self.c2).get_nuisance()
193  # tau = IMP.isd.Nuisance(self.tau).get_nuisance()
194  tau = 1.0
195  self.cov = IMP.isd2.compute_relative_covariance(self.atoms, c1, c2,
196  tau, self.prof)
197  # for i in xrange(len(self.cov)):
198  # for j in xrange(len(self.cov)):
199  # self.cov[i][j] = self.cov[i][j]/tau**2
200  self.saxs.set_cov(0, self.cov)
201 
202  def write_covariance_matrix(self, fname):
203  fl = open(fname, 'w')
204  for line in self.cov:
205  for i in line:
206  fl.write('%G ' % i)
207  fl.write('\n')
208 
209  def get_output(self):
210  output = super(SAXSISDRestraint, self).get_output()
211  suffix = self._get_label_suffix()
212  output["SAXSISDRestraint_Sigma" +
213  suffix] = str(self.sigma.get_scale())
214  output["SAXSISDRestraint_Tau" + suffix] = str(self.tau.get_scale())
215  output["SAXSISDRestraint_Gamma" +
216  suffix] = str(self.gamma.get_scale())
217  return output
218 
219  @staticmethod
220  def _logspace(a, b, num=100):
221  """Mimick numpy's logspace function"""
222  for i in range(num):
223  val = a + float(b - a) / float(num - 1) * i
224  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.
Definition: saxs.py:141
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.
Basic SAXS restraint.
Definition: saxs.py:17
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: tools.py:1579
def optimize_sigma
Set sigma to the value that maximizes its conditional likelihood.
Definition: saxs.py:134
def draw_gamma
Draw gamma from lognormal distribution.
Definition: saxs.py:185
Add nuisance parameter to particle.
Definition: Nuisance.h:25
Basic SAXS restraint using ISD.
Definition: saxs.py:65
def get_gamma_value
Get value of gamma.
Definition: saxs.py:173
def __init__
Builds the restraint.
Definition: saxs.py:33
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.
Definition: saxs.py:180
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:700
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.
Base class for PMI restraints, which wrap IMP.Restraint(s).
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...