IMP logo
IMP Reference Guide  2.6.1
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 IMP
7 import IMP.core
8 import IMP.algebra
9 import IMP.atom
10 import IMP.container
11 import IMP.pmi.tools
12 
13 class SAXSISDRestraint(object):
14 
15  import IMP.saxs
16  import IMP.isd
17  import IMP.isd2
18  import IMP.pmi.tools
19 
20  def __init__(self, representation, profile, resolution=0, weight=1,
21  ff_type=IMP.saxs.HEAVY_ATOMS):
22 
23  self.m = representation.prot.get_model()
24  self.label = "None"
25  self.rs = IMP.RestraintSet(self.m, 'saxs')
26 
27  self.taumaxtrans = 0.05
28  self.prof = IMP.saxs.Profile(profile)
29 
30  self.atoms = IMP.pmi.tools.select(
31  representation,
32  resolution=resolution)
33 
34  # gamma nuisance
35  self.gamma = IMP.pmi.tools.SetupNuisance(
36  self.m, 1., 0., None, False).get_particle()
37 
38  # sigma nuisance
39  self.sigma = IMP.pmi.tools.SetupNuisance(self.m, 10.0, 0., None, False
40  ).get_particle()
41 
42  # tau nuisance, optimized
43  self.tau = IMP.pmi.tools.SetupNuisance(self.m, 1., 0., None, False,
44  ).get_particle()
45 
46  #c1 and c2, optimized
47  self.c1 = IMP.pmi.tools.SetupNuisance(self.m, 1.0, 0.95, 1.05,
48  True).get_particle()
49  self.c2 = IMP.pmi.tools.SetupNuisance(self.m, 0.0, -2., 4.,
50  True).get_particle()
51 
52  #weight, optimized
53  self.w = IMP.pmi.tools.SetupWeight(self.m).get_particle()
54  IMP.isd.Weight(self.w).set_weights_are_optimized(True)
55 
56  # take identity covariance matrix for the start
57  self.cov = [[1 if i == j else 0 for j in range(self.prof.size())]
58  for i in range(self.prof.size())]
59 
60  print("create saxs restraint")
61  self.saxs = IMP.isd2.SAXSRestraint(self.prof, self.sigma, self.tau,
62  self.gamma, self.w, self.c1, self.c2)
63  self.saxs.add_scatterer(self.atoms, self.cov, ff_type)
64 
65  self.rs.add_restraint(self.saxs)
66  self.rs.set_weight(weight)
67 
68  # self.saxs_stuff={'nuis':(sigma,gamma),'cov':cov,
69  # 'exp':prof,'th':tmp}
70 
71  self.rs2 = IMP.RestraintSet(self.m, 'jeffreys')
72  # jeffreys restraints for nuisances
73  j1 = IMP.isd.JeffreysRestraint(self.m, self.sigma)
74  self.rs2.add_restraint(j1)
75  j2 = IMP.isd.JeffreysRestraint(self.m, self.tau)
76  self.rs2.add_restraint(j2)
77  j3 = IMP.isd.JeffreysRestraint(self.m, self.gamma)
78  self.rs2.add_restraint(j3)
79 
80  def optimize_sigma(self):
81  """set sigma to the value that maximizes its conditional likelihood"""
82  from math import sqrt
83  self.m.update()
84  sigma2hat = self.saxs.get_sigmasq_scale_parameter() \
85  / (self.saxs.get_sigmasq_shape_parameter() + 1)
86  IMP.isd.Scale(self.sigma).set_scale(sqrt(sigma2hat))
87 
88  def optimize_gamma(self):
89  """set gamma to the value that maximizes its conditional likelihood"""
90  from math import exp
91  self.m.update()
92  gammahat = exp(self.saxs.get_loggamma_variance_parameter()
93  * self.saxs.get_loggamma_jOg_parameter())
94  IMP.isd.Scale(self.gamma).set_scale(gammahat)
95 
96  def logspace(self, a, b, num=100):
97  """mimick numpy's logspace function"""
98  for i in range(num):
99  val = a + float(b - a) / float(num - 1) * i
100  yield 10 ** val
101 
102  def optimize_tau(self, ltaumin=-2, ltaumax=3, npoints=100):
103  from math import log
104  import IMP.atom
105  values = []
106  self.m.update()
107  IMP.atom.write_pdb(self.atoms, 'tauvals.pdb')
108  fl = open('tauvals.txt', 'w')
109  for tauval in self.logspace(ltaumin, ltaumax, npoints):
110  IMP.isd.Scale(self.tau).set_scale(tauval)
111  try:
112  values.append((self.m.evaluate(False), tauval))
113  except:
114  pass
115  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
116  values.sort()
117  ltcenter = log(values[0][1]) / log(10)
118  spacing = (ltaumax - ltaumin) / float(npoints)
119  values = []
120  for tauval in self.logspace(
121  ltcenter - 2 * spacing, ltcenter + 2 * spacing,
122  npoints):
123  IMP.isd.Scale(self.tau).set_scale(tauval)
124  values.append((self.m.evaluate(False), tauval))
125  fl.write('%G %G\n' % (values[-1][1], values[-1][0]))
126  values.sort()
127  IMP.isd.Scale(self.tau).set_scale(values[0][1])
128 
129  def draw_sigma(self):
130  """draw 1/sigma2 from gamma distribution"""
131  self.m.update()
132  self.saxs.draw_sigma()
133 
134  def draw_gamma(self):
135  """draw gamma from lognormal distribution"""
136  self.m.update()
137  self.saxs.draw_gamma()
138 
139  def update_covariance_matrix(self):
140  c1 = IMP.isd.Nuisance(self.c1).get_nuisance()
141  c2 = IMP.isd.Nuisance(self.c2).get_nuisance()
142  #tau = IMP.isd.Nuisance(self.tau).get_nuisance()
143  tau = 1.0
144  self.cov = IMP.isd2.compute_relative_covariance(self.atoms, c1, c2,
145  tau, self.prof)
146  # for i in xrange(len(self.cov)):
147  # for j in xrange(len(self.cov)):
148  # self.cov[i][j] = self.cov[i][j]/tau**2
149  self.saxs.set_cov(0, self.cov)
150 
151  def write_covariance_matrix(self, fname):
152  fl = open(fname, 'w')
153  for line in self.cov:
154  for i in line:
155  fl.write('%G ' % i)
156  fl.write('\n')
157 
158  def get_gamma_value(self):
159  return self.gamma.get_scale()
160 
161  def set_label(self, label):
162  self.label = label
163 
164  def add_to_model(self):
165  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
166  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs2)
167 
168  def set_taumaxtrans(self, taumaxtrans):
169  self.taumaxtrans = taumaxtrans
170 
171  def get_particles_to_sample(self):
172  ps = {}
173  # ps["Nuisances_SAXSISDRestraint_Tau_" +
174  # self.label] = ([self.tau], self.taumaxtrans)
175  return ps
176 
177  def get_output(self):
178  self.m.update()
179  output = {}
180  score = self.rs.unprotected_evaluate(None)
181  score2 = self.rs2.unprotected_evaluate(None)
182  output["_TotalScore"] = str(score + score2)
183 
184  output["SAXSISDRestraint_Likelihood_" + self.label] = str(score)
185  output["SAXSISDRestraint_Prior_" + self.label] = str(score2)
186  output["SAXSISDRestraint_Sigma_" +
187  self.label] = str(self.sigma.get_scale())
188  output["SAXSISDRestraint_Tau_" +
189  self.label] = str(self.tau.get_scale())
190  output["SAXSISDRestraint_Gamma_" +
191  self.label] = str(self.gamma.get_scale())
192  return output
Add weights for a set of states to a particle.
Definition: Weight.h:24
Various classes to hold sets of particles.
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
log
Definition: log.py:1
Object used to hold a set of restraints.
Definition: RestraintSet.h:36
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:37
Add nuisance parameter to particle.
Definition: Nuisance.h:25
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...
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:702
Functionality for loading, creating, manipulating and scoring atomic structures.
Support for small angle X-ray scattering (SAXS) data.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...