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
115 self.atoms = IMP.pmi.tools.select(
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 pw.set_weights(pw.get_unit_simplex().get_barycenter())
167 pw.set_weights_are_optimized(
True)
170 """Set sigma to the value that maximizes its conditional likelihood"""
172 sigma2hat = self.saxs.get_sigmasq_scale_parameter() \
173 / (self.saxs.get_sigmasq_shape_parameter() + 1)
177 """Set gamma to the value that maximizes its conditional likelihood"""
179 gammahat = math.exp(self.saxs.get_loggamma_variance_parameter() *
180 self.saxs.get_loggamma_jOg_parameter())
183 def optimize_tau(self, ltaumin=-2, ltaumax=3, npoints=100):
187 fl = open(
'tauvals.txt',
'w')
188 for tauval
in self._logspace(ltaumin, ltaumax, npoints):
191 values.append((self.model.evaluate(
False), tauval))
194 fl.write(
'%G %G\n' % (values[-1][1], values[-1][0]))
196 ltcenter = math.log(values[0][1]) / math.log(10)
197 spacing = (ltaumax - ltaumin) / float(npoints)
199 for tauval
in self._logspace(
200 ltcenter - 2 * spacing, ltcenter + 2 * spacing,
203 values.append((self.model.evaluate(
False), tauval))
204 fl.write(
'%G %G\n' % (values[-1][1], values[-1][0]))
209 """Get value of gamma."""
210 return self.gamma.get_scale()
212 def set_taumaxtrans(self, taumaxtrans):
213 self.taumaxtrans = taumaxtrans
216 """Draw 1/sigma2 from gamma distribution."""
218 self.saxs.draw_sigma()
221 """Draw gamma from lognormal distribution."""
223 self.saxs.draw_gamma()
225 def update_covariance_matrix(self):
230 self.cov = IMP.isd2.compute_relative_covariance(self.atoms, c1, c2,
235 self.saxs.set_cov(0, self.cov)
237 def write_covariance_matrix(self, fname):
238 fl = open(fname,
'w')
239 for line
in self.cov:
245 output = super(SAXSISDRestraint, self).
get_output()
246 suffix = self._get_label_suffix()
247 output[
"SAXSISDRestraint_Sigma" +
248 suffix] = str(self.sigma.get_scale())
249 output[
"SAXSISDRestraint_Tau" + suffix] = str(self.tau.get_scale())
250 output[
"SAXSISDRestraint_Gamma" +
251 suffix] = str(self.gamma.get_scale())
255 def _logspace(a, b, num=100):
256 """Mimick numpy's logspace function"""
258 val = a + float(b - a) / float(num - 1) * i
Add weights 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 ...