Changeset 4f2478e in sasmodels

Timestamp:

Jan 28, 2016 1:19:03 PM (9 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

a3e78c3

Parents:

bf6631a

Message:

simplify calculation for sph_j1c

Location:

sasmodels

Files:

• compare.py (modified) (4 diffs)
• compare_many.py (modified) (1 diff)
• models/lib/sph_j1c.c (modified) (3 diffs)

sasmodels/compare.py

@@ -117 +117 @@
         return dt.total_seconds()
 
+
+class push_seed(object):
+    """
+    Set the seed value for the random number generator.
+
+    When used in a with statement, the random number generator state is
+    restored after the with statement is complete.
+
+    :Parameters:
+
+    *seed* : int or array_like, optional
+        Seed for RandomState
+
+    :Example:
+
+    Seed can be used directly to set the seed::
+
+        >>> from numpy.random import randint
+        >>> push_seed(24)
+        <...push_seed object at...>
+        >>> print(randint(0,1000000,3))
+        [242082    899 211136]
+
+    Seed can also be used in a with statement, which sets the random
+    number generator state for the enclosed computations and restores
+    it to the previous state on completion::
+
+        >>> with push_seed(24):
+        ...    print(randint(0,1000000,3))
+        [242082    899 211136]
+
+    Using nested contexts, we can demonstrate that state is indeed
+    restored after the block completes::
+
+        >>> with push_seed(24):
+        ...    print(randint(0,1000000))
+        ...    with push_seed(24):
+        ...        print(randint(0,1000000,3))
+        ...    print(randint(0,1000000))
+        242082
+        [242082    899 211136]
+        899
+
+    The restore step is protected against exceptions in the block::
+
+        >>> with push_seed(24):
+        ...    print(randint(0,1000000))
+        ...    try:
+        ...        with push_seed(24):
+        ...            print(randint(0,1000000,3))
+        ...            raise Exception()
+        ...    except:
+        ...        print("Exception raised")
+        ...    print(randint(0,1000000))
+        242082
+        [242082    899 211136]
+        Exception raised
+        899
+    """
+    def __init__(self, seed=None):
+        self._state = np.random.get_state()
+        np.random.seed(seed)
+
+    def __enter__(self):
+        return None
+
+    def __exit__(self, *args):
+        np.random.set_state(self._state)
 
 def tic():
@@ -179 +247 @@
         return [0, (2*v if v > 0 else 1)]
 
+
 def _randomize_one(p, v):
     """
@@ -188 +257 @@
         return np.random.uniform(*parameter_range(p, v))
 
+
 def randomize_pars(pars, seed=None):
     """
@@ -197 +267 @@
     :func:`constrain_pars` needs to be called afterward..
     """
-    np.random.seed(seed)
-    # Note: the sort guarantees order `of calls to random number generator
-    pars = dict((p, _randomize_one(p, v))
-                for p, v in sorted(pars.items()))
+    with push_seed(seed):
+        # Note: the sort guarantees order `of calls to random number generator
+        pars = dict((p, _randomize_one(p, v))
+                    for p, v in sorted(pars.items()))
     return pars
 

sasmodels/compare_many.py

@@ -261 +261 @@
 
 if __name__ == "__main__":
+    #from .compare import push_seed
+    #with push_seed(1): main()
     main()

sasmodels/models/lib/sph_j1c.c

@@ -1 +1 @@
 /**
-* Spherical Bessel function j1(x)/x
+* Spherical Bessel function 3*j1(x)/x
 *
 * Used for low q to avoid cancellation error.
@@ -8 +8 @@
 * requires the first 3 terms.  Double precision requires the 4th term.
 * The fifth term is not needed, and is commented out.
+* Taylor expansion:
+*      1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.))))
+* Expression returned from Herbie (herbie.uwpise.org/demo):
+*      const double t = ((1. + 3.*q2*q2/5600.) - q2/20.);
+*      return t*t;
 */
 
@@ -18 +23 @@
     SINCOS(q, sin_q, cos_q);
 
-    const double bessel = (q < 1.e-1) ?
-        1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))) // + q2*(3./3991680.)))
+    const double bessel = (q < 0.384038453352533)
+        ? (1.0 + q2*(-3./30. + q2*(3./840.))
         : 3.0*(sin_q/q - cos_q)/q2;
 
     return bessel;
+
+ /*
+    // Code to test various expressions
+    if (sizeof(q2) > 4) {
+        return 3.0*(sin_q/q - cos_q)/q2;
+    } else if (q < 0.384038453352533) {
+        //const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); return t*t;
+        return 1.0 + q2*q2*(3./840.) - q2*(3./30.);
+        //return 1.0 + q2*(-3./30. + q2*(3./840.));
+        //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.)));
+        //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.))));
+    } else {
+        return 3.0*(sin_q/q - cos_q)/q2;
+    }
+*/
 }