1 """@namespace IMP.pmi.restraints.saxs
2 Restraints for handling small angle x-ray (SAXS) data.
5 from __future__
import print_function
24 """Basic SAXS restraint."""
26 _include_in_rmf =
True
28 def __init__(self, input_objects, saxs_datafile, weight=1.0,
29 ff_type=IMP.saxs.HEAVY_ATOMS, label=
None, maxq=
"standard"):
30 """Builds the restraint.
31 @param input_objects A list of hierarchies or PMI objects that the
32 SAXS restraint will be applied to. This hierarchy must be
33 atomic unless ff_type=IMP.saxs.RESIDUES is used.
34 @param saxs_datafile the SAXS .dat file.
35 @param weight Restraint score coefficient
36 @param ff_type the form factor to use, of the following types:
37 - IMP.saxs.HEAVY_ATOMS: use form factors with implicit
39 - IMP.saxs.ALL_ATOMS: use individual form factors for all
40 atoms. Does not build missing hydrogens.
41 - IMP.saxs.CA_ATOMS: use residue based form factors
43 - IMP.saxs.RESIDUES: use residue based form factors
44 using per-residue beads
45 @param label Label for the restraint in outputs
46 @param maxq Maximum q value that the restraint will be evaluated at.
47 If set to 'standard' (the default), the following values will
48 be used (these values were eyeballed by comparing ALL_ATOM
49 calculated SAXS profiles to those calculated with the reduced
50 representations, so could be improved):
51 - For ff_type = ALL_ATOMS: 0.5
53 - CA_ATOMS and RESIDUES: 0.15
58 model = list(hiers)[0].get_model()
59 super(SAXSRestraint, self).
__init__(model, label=label, weight=weight)
62 if maxq ==
"standard":
63 if ff_type == IMP.saxs.CA_ATOMS
or ff_type == IMP.saxs.RESIDUES:
65 elif ff_type == IMP.saxs.HEAVY_ATOMS:
69 elif isinstance(maxq, float):
70 if maxq < 0.01
or maxq > 4.0:
72 "SAXSRestraint: maxq must be set between 0.01 and 4.0")
73 if (ff_type == IMP.saxs.CA_ATOMS
or ff_type == IMP.saxs.RESIDUES) \
75 warnings.warn(
"SAXSRestraint: for residue-resolved form "
76 "factors, a maxq > 0.15 is not recommended!",
80 "SAXSRestraint: maxq must be set to a number between 0.01 "
86 if ff_type == IMP.saxs.RESIDUES:
88 hiers, resolution=1).get_selected_particles()
90 elif ff_type == IMP.saxs.CA_ATOMS:
92 hiers, atom_type=IMP.atom.AT_CA).get_selected_particles()
94 elif ff_type == IMP.saxs.HEAVY_ATOMS
or ff_type == IMP.saxs.ALL_ATOMS:
96 hiers, resolution=0).get_selected_particles()
100 "SAXSRestraint: Must provide an IMP.saxs atom type: "
101 "RESIDUES, CA_ATOMS, HEAVY_ATOMS or ALL_ATOMS")
103 if len(self.particles) == 0:
104 raise Exception(
"SAXSRestraint: There are no selected particles")
108 self.rs.add_restraint(self.restraint)
113 """Basic SAXS restraint using ISD."""
115 def __init__(self, representation, profile, resolution=0, weight=1,
116 ff_type=IMP.saxs.HEAVY_ATOMS, label=
None):
118 if not hasattr(IMP,
'isd2'):
119 raise ImportError(
"Module isd2 not installed. "
120 "Cannot use SAXSISDRestraint")
122 model = representation.prot.get_model()
123 super(SAXSISDRestraint, self).__init__(model, label=label,
126 self.taumaxtrans = 0.05
129 self.atoms = IMP.pmi.tools.select(
131 resolution=resolution)
134 self.gamma = IMP.pmi.tools.SetupNuisance(
135 self.model, 1., 0.,
None,
False).get_particle()
138 self.sigma = IMP.pmi.tools.SetupNuisance(
139 self.model, 10.0, 0.,
None,
False).get_particle()
142 self.tau = IMP.pmi.tools.SetupNuisance(self.model, 1., 0.,
None,
False,
146 self.c1 = IMP.pmi.tools.SetupNuisance(self.model, 1.0, 0.95, 1.05,
148 self.c2 = IMP.pmi.tools.SetupNuisance(self.model, 0.0, -2., 4.,
152 self.w = IMP.pmi.tools.SetupWeight(self.model).get_particle()
156 self.cov = [[1
if i == j
else 0
for j
in range(self.prof.size())]
157 for i
in range(self.prof.size())]
159 print(
"create saxs restraint")
160 self.saxs = IMP.isd2.SAXSRestraint(self.prof, self.sigma, self.tau,
161 self.gamma, self.w, self.c1,
163 self.saxs.add_scatterer(self.atoms, self.cov, ff_type)
165 self.rs.add_restraint(self.saxs)
170 self.rs2 = self._create_restraint_set(
'Prior')
173 self.rs2.add_restraint(j1)
175 self.rs2.add_restraint(j2)
177 self.rs2.add_restraint(j3)
180 pw.set_weights(pw.get_unit_simplex().get_barycenter())
181 pw.set_weights_are_optimized(
True)
184 """Set sigma to the value that maximizes its conditional likelihood"""
186 sigma2hat = self.saxs.get_sigmasq_scale_parameter() \
187 / (self.saxs.get_sigmasq_shape_parameter() + 1)
191 """Set gamma to the value that maximizes its conditional likelihood"""
193 gammahat = math.exp(self.saxs.get_loggamma_variance_parameter() *
194 self.saxs.get_loggamma_jOg_parameter())
197 def optimize_tau(self, ltaumin=-2, ltaumax=3, npoints=100):
201 fl = open(
'tauvals.txt',
'w')
202 for tauval
in self._logspace(ltaumin, ltaumax, npoints):
205 values.append((self.model.evaluate(
False), tauval))
208 fl.write(
'%G %G\n' % (values[-1][1], values[-1][0]))
210 ltcenter = math.log(values[0][1]) / math.log(10)
211 spacing = (ltaumax - ltaumin) / float(npoints)
213 for tauval
in self._logspace(
214 ltcenter - 2 * spacing, ltcenter + 2 * spacing,
217 values.append((self.model.evaluate(
False), tauval))
218 fl.write(
'%G %G\n' % (values[-1][1], values[-1][0]))
223 """Get value of gamma."""
224 return self.gamma.get_scale()
226 def set_taumaxtrans(self, taumaxtrans):
227 self.taumaxtrans = taumaxtrans
230 """Draw 1/sigma2 from gamma distribution."""
232 self.saxs.draw_sigma()
235 """Draw gamma from lognormal distribution."""
237 self.saxs.draw_gamma()
239 def update_covariance_matrix(self):
244 self.cov = IMP.isd2.compute_relative_covariance(self.atoms, c1, c2,
249 self.saxs.set_cov(0, self.cov)
251 def write_covariance_matrix(self, fname):
252 fl = open(fname,
'w')
253 for line
in self.cov:
258 def get_output(self):
259 output = super(SAXSISDRestraint, self).get_output()
260 suffix = self._get_label_suffix()
261 output[
"SAXSISDRestraint_Sigma" +
262 suffix] = str(self.sigma.get_scale())
263 output[
"SAXSISDRestraint_Tau" + suffix] = str(self.tau.get_scale())
264 output[
"SAXSISDRestraint_Gamma" +
265 suffix] = str(self.gamma.get_scale())
269 def _logspace(a, b, num=100):
270 """Mimic numpy's logspace function"""
272 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.
Support for small angle X-ray scattering (SAXS) data.
Warning for probably incorrect input parameters.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...