[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[IMP-dev] Fix behavior of DistanceRestraint at near-zero distance



Attached is a patch which corrects the math in DistanceRestraint. Currently if the interparticle distance is zero (or very near zero) DistanceRestraint pretends that the particles are randomly placed in near proximity, so that a random distance and derivatives are returned. This is because the derivatives of the distance at zero are undefined (the math is 0/0).

IMHO, this behavior is flawed, because the score of a given state should always be the same - and the scoring function should be deterministic. Any necessary randomization of the system should be introduced by perturbing the system state before an optimization, not by having the scoring function return random values. (The current behavior also makes it impossible to introduce a distance restraint with zero mean, to force two particles to be colocated.)

This patch returns a zero derivative for zero distance, which is at least consistent with every molecular mechanics package I've seen.

	Ben
--
                      http://salilab.org/~ben/
"It is a capital mistake to theorize before one has data."
	- Sir Arthur Conan Doyle
Index: test/sd_optimizer/test_sd_optimizer.py
===================================================================
--- test/sd_optimizer/test_sd_optimizer.py	(revision 629)
+++ test/sd_optimizer/test_sd_optimizer.py	(working copy)
@@ -1,6 +1,7 @@
 import unittest
 import IMP.utils
 import IMP.test, IMP
+import random
 
 class SteepestDescentTests(IMP.test.IMPTestCase):
     """Test steepest descent optimizer"""
@@ -64,18 +65,13 @@
     def test_sd_optimizer2(self):
         """Test that optimizer spreads particles apart"""
 
-        self.particles[0].set_x(0.0)
-        self.particles[0].set_y(0.0)
-        self.particles[0].set_z(0.0)
+        # Start off with all particles in close proximity (but not actually
+        # colocated, as the derivative of zero distance is zero):
+        for p in self.particles:
+            p.set_x(random.uniform(-0.01, 0.01))
+            p.set_y(random.uniform(-0.01, 0.01))
+            p.set_z(random.uniform(-0.01, 0.01))
 
-        self.particles[1].set_x(0.0)
-        self.particles[1].set_y(0.0)
-        self.particles[1].set_z(0.0)
-
-        self.particles[2].set_x(0.0)
-        self.particles[2].set_y(0.0)
-        self.particles[2].set_z(0.0)
-
         self.steepest_descent.optimize(50)
 
         for i in range(0, 2):
Index: src/restraints/DistanceRestraint.cpp
===================================================================
--- src/restraints/DistanceRestraint.cpp	(revision 629)
+++ src/restraints/DistanceRestraint.cpp	(working copy)
@@ -98,20 +98,6 @@
   // calculate the distance feature
   distance = sqrt(delta_x * delta_x + delta_y * delta_y + delta_z * delta_z);
 
-  // if distance is too close to zero, set it to some non-zero value
-  if (distance < DistanceRestraint::MIN_DISTANCE) {
-    delta_x = std::rand(); // arbitrary move
-    delta_y = std::rand(); // arbitrary move
-    delta_z = std::rand(); // arbitrary move
-    distance = sqrt(delta_x * delta_x + delta_y * delta_y + delta_z * delta_z);
-
-    // normalize the random move, to the min disance
-    delta_x = DistanceRestraint::MIN_DISTANCE * delta_x / distance;
-    delta_y = DistanceRestraint::MIN_DISTANCE * delta_y / distance;
-    delta_z = DistanceRestraint::MIN_DISTANCE * delta_z / distance;
-    distance = DistanceRestraint::MIN_DISTANCE;
-  }
-
   // if needed, calculate the partial derivatives of the scores with respect
   // to the particle attributes
   if (accum) {
@@ -121,15 +107,21 @@
     // calculate the score and feature derivative based on the distance feature
     score = (*score_func_)(distance, deriv);
 
-    dx = delta_x / distance * deriv;
+    if (distance >= DistanceRestraint::MIN_DISTANCE) {
+      dx = delta_x / distance * deriv;
+      dy = delta_y / distance * deriv;
+      dz = delta_z / distance * deriv;
+    } else {
+      // avoid division by zero
+      dx = dy = dz = 0.;
+    }
+
     get_particle(0)->add_to_derivative(x_, dx, *accum);
     get_particle(1)->add_to_derivative(x_, -dx, *accum);
 
-    dy = delta_y / distance * deriv;
     get_particle(0)->add_to_derivative(y_, dy, *accum);
     get_particle(1)->add_to_derivative(y_, -dy, *accum);
 
-    dz = delta_z / distance * deriv;
     get_particle(0)->add_to_derivative(z_, dz, *accum);
     get_particle(1)->add_to_derivative(z_, -dz, *accum);
   }