Improvements to random number routines used in constant-pH MD 37/4337/1
authorradakb <brian.radak@gmail.com>
Tue, 10 Jul 2018 01:32:09 +0000 (21:32 -0400)
committerradakb <brian.radak@gmail.com>
Tue, 10 Jul 2018 01:32:09 +0000 (21:32 -0400)
- Some users reported a bug when sampling from a PMF with a
  dominant mode. Underflow during subtraction meant that no
  choice was made. A catch has been added that finds and returns
  the modal value in this case.

- The normal() code has been changed to use the iterative
  Box-Muller-like transformation that is implemented in Random.h
  The new code also stores the second normal variable and returns
  it later rather than just discarding it.

Change-Id: I7d66b93b7135665db91208412d482b5a47ed268d

lib/namdcph/namdcph/numtcl.tcl

index 166bc12..ba43de6 100644 (file)
@@ -150,17 +150,32 @@ proc lmultiply {list args} {
 # Generate a Gaussian (normal) random variable with the given mean and standard
 # deviation (default to zero and one, respectively). 
 #
-# This uses the standard Box-Muller transformation of two uniform random 
-# variables on (0, 1).
+# This uses a Box-Muller-like transformation by generating two uniform variates
+# within the unit circle. The cost of rejecting some random numbers generally
+# outweighs the cost of the cos() and sin() needed for the standard Box-Muller
+# transformation.
 #
+variable SecondGaussianWaiting 0
+variable SecondGaussian 0.0
 proc tcl::mathfunc::normal {{mu 0.0} {sigma 1.0}} {
-    set two_pi [expr {2*3.141592653589793}]
-    set u1 [expr {rand()}]
-    set u2 [expr {rand()}]
-    # The transformation produces two random variates - either of these can be
-    # used (or both).
-    expr {$mu + $sigma*sqrt(-2*log($u1))*cos($two_pi*$u2)}
-    # expr {$mu + $sigma*sqrt(-2*log($u1))*sin($two_pi*$u2)}
+    if {$::SecondGaussianWaiting} {
+        set ::SecondGaussianWaiting 0
+        expr {$mu + $sigma*$::SecondGaussian}
+    } else {
+        # Pick two variates within the unit circle. The lower bound is from
+        # NAMD Random.h and purportedly ensures "abs result < 6".
+        set r 2.0
+        while {$r >=1. || $r < 1.523e-8} {
+            set u1 [expr {2*rand() - 1}]
+            set u2 [expr {2*rand() - 1}]
+            set r [expr {$u1*$u1 + $u2*$u2}]
+        }
+        # The transformation produces two random variates - store one.
+        set fac [expr sqrt(-2*log($r)/$r)]
+        set ::SecondGaussianWaiting 1
+        set ::SecondGaussian [expr {$fac*$u1}]
+        expr {$mu + $sigma*$fac*$u2}
+    }
 }
 
 # choice
@@ -205,7 +220,18 @@ proc choice {list {weights {}}} {
         }
         incr j
     }
-    # This should never be reached.
-    error "Something bad happened in choice!"
+    # This may be reached if one of the weights is essentially one within
+    # rounding precision. In that case, returning the most probable choice is
+    # likely a good approximation and really the only stable thing to do.
+    set MaxWeight 0.0
+    set i 0
+    foreach Weight $weights {
+        if {$Weight >= $MaxWeight} {
+            set MaxWeight $Weight
+            set j $i
+        }
+        incr i
+    }
+    return [list [lindex $list $j] $j]
 }