1 """@namespace IMP.pmi.restraints.saxs
2 Restraints for handling small angle x-ray (SAXS) data.
5 from __future__
import print_function
19 """Basic SAXS restraint."""
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
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
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
37 @param label Label for the restraint in outputs
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.
48 model = list(hiers)[0].get_model()
49 super(SAXSRestraint, self).
__init__(model, label=label, weight=weight)
52 if maxq ==
"standard":
53 if ff_type == IMP.saxs.CA_ATOMS
or ff_type == IMP.saxs.RESIDUES:
55 elif ff_type == IMP.saxs.HEAVY_ATOMS:
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!",
67 raise Exception(
"SAXSRestraint: maxq must be set to a number between 0.01 and 4.0")
72 if ff_type == IMP.saxs.RESIDUES:
74 hiers, resolution=1).get_selected_particles()
76 elif ff_type == IMP.saxs.CA_ATOMS:
78 hiers, atom_type=IMP.atom.AT_CA).get_selected_particles()
80 elif ff_type == IMP.saxs.HEAVY_ATOMS
or ff_type == IMP.saxs.ALL_ATOMS:
82 hiers, resolution=0).get_selected_particles()
85 raise Exception(
"SAXSRestraint: Must provide an IMP.saxs atom type: RESIDUES, CA_ATOMS, HEAVY_ATOMS or ALL_ATOMS")
87 if len(self.particles) == 0:
88 raise Exception(
"SAXSRestraint: There are no selected particles")
92 self.rs.add_restraint(self.restraint)
97 """Basic SAXS restraint using ISD."""
103 print(
"Module isd2 not installed. Cannot use SAXSISDRestraint")
105 def __init__(self, representation, profile, resolution=0, weight=1,
106 ff_type=IMP.saxs.HEAVY_ATOMS, label=
None):
108 model = representation.prot.get_model()
109 super(SAXSISDRestraint, self).
__init__(model, label=label,
112 self.taumaxtrans = 0.05
117 resolution=resolution)
120 self.gamma = IMP.pmi.tools.SetupNuisance(
121 self.model, 1., 0.,
None,
False).get_particle()
124 self.sigma = IMP.pmi.tools.SetupNuisance(self.model, 10.0, 0.,
None,
False
128 self.tau = IMP.pmi.tools.SetupNuisance(self.model, 1., 0.,
None,
False,
132 self.c1 = IMP.pmi.tools.SetupNuisance(self.model, 1.0, 0.95, 1.05,
134 self.c2 = IMP.pmi.tools.SetupNuisance(self.model, 0.0, -2., 4.,
138 self.w = IMP.pmi.tools.SetupWeight(self.model).get_particle()
142 self.cov = [[1
if i == j
else 0
for j
in range(self.prof.size())]
143 for i
in range(self.prof.size())]
145 print(
"create saxs restraint")
146 self.saxs = IMP.isd2.SAXSRestraint(self.prof, self.sigma, self.tau,
147 self.gamma, self.w, self.c1,
149 self.saxs.add_scatterer(self.atoms, self.cov, ff_type)
151 self.rs.add_restraint(self.saxs)
156 self.rs2 = self._create_restraint_set(
'Prior')
159 self.rs2.add_restraint(j1)
161 self.rs2.add_restraint(j2)
163 self.rs2.add_restraint(j3)
166 """Set sigma to the value that maximizes its conditional likelihood"""
168 sigma2hat = self.saxs.get_sigmasq_scale_parameter() \
169 / (self.saxs.get_sigmasq_shape_parameter() + 1)
173 """Set gamma to the value that maximizes its conditional likelihood"""
175 gammahat = math.exp(self.saxs.get_loggamma_variance_parameter() *
176 self.saxs.get_loggamma_jOg_parameter())
179 def optimize_tau(self, ltaumin=-2, ltaumax=3, npoints=100):
183 fl = open(
'tauvals.txt',
'w')
184 for tauval
in self._logspace(ltaumin, ltaumax, npoints):
187 values.append((self.model.evaluate(
False), tauval))
190 fl.write(
'%G %G\n' % (values[-1][1], values[-1][0]))
192 ltcenter = math.log(values[0][1]) / math.log(10)
193 spacing = (ltaumax - ltaumin) / float(npoints)
195 for tauval
in self._logspace(
196 ltcenter - 2 * spacing, ltcenter + 2 * spacing,
199 values.append((self.model.evaluate(
False), tauval))
200 fl.write(
'%G %G\n' % (values[-1][1], values[-1][0]))
205 """Get value of gamma."""
206 return self.gamma.get_scale()
208 def set_taumaxtrans(self, taumaxtrans):
209 self.taumaxtrans = taumaxtrans
212 """Draw 1/sigma2 from gamma distribution."""
214 self.saxs.draw_sigma()
217 """Draw gamma from lognormal distribution."""
219 self.saxs.draw_gamma()
221 def update_covariance_matrix(self):
226 self.cov = IMP.isd2.compute_relative_covariance(self.atoms, c1, c2,
231 self.saxs.set_cov(0, self.cov)
233 def write_covariance_matrix(self, fname):
234 fl = open(fname,
'w')
235 for line
in self.cov:
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())
251 def _logspace(a, b, num=100):
252 """Mimick numpy's logspace function"""
254 val = a + float(b - a) / float(num - 1) * i
Add weights for a set of states to a particle.
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.
void write_pdb(const Selection &mhd, TextOutput out, unsigned int model=1)
Add scale parameter to particle.
Classes to handle different kinds of restraints.
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.
Basic SAXS restraint using ISD.
def get_gamma_value
Get value of gamma.
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.
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.
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 ...