[Numpy-svn] r6000 - in trunk/numpy/core: code_generators src tests

numpy-svn@scip... numpy-svn@scip...
Mon Nov 10 22:42:34 CST 2008


Author: charris
Date: 2008-11-10 22:42:25 -0600 (Mon, 10 Nov 2008)
New Revision: 6000

Modified:
   trunk/numpy/core/code_generators/generate_umath.py
   trunk/numpy/core/src/umathmodule.c.src
   trunk/numpy/core/tests/test_umath.py
Log:
Add logaddexp2.
Add tests for log, exp, logaddexp2.

Modified: trunk/numpy/core/code_generators/generate_umath.py
===================================================================
--- trunk/numpy/core/code_generators/generate_umath.py	2008-11-11 03:26:53 UTC (rev 5999)
+++ trunk/numpy/core/code_generators/generate_umath.py	2008-11-11 04:42:25 UTC (rev 6000)
@@ -333,6 +333,11 @@
           "",
           TD(flts, f="logaddexp")
           ),
+'logaddexp2' :
+    Ufunc(2, 1, None,
+          "",
+          TD(flts, f="logaddexp2")
+          ),
 'bitwise_and' :
     Ufunc(2, 1, One,
           docstrings.get('numpy.core.umath.bitwise_and'),

Modified: trunk/numpy/core/src/umathmodule.c.src
===================================================================
--- trunk/numpy/core/src/umathmodule.c.src	2008-11-11 03:26:53 UTC (rev 5999)
+++ trunk/numpy/core/src/umathmodule.c.src	2008-11-11 04:42:25 UTC (rev 6000)
@@ -35,7 +35,11 @@
  * #C = F, ,L#
  */
 
+/* fixme: need more precision for LOG2 and INVLOG2 */
+
 #define PI 3.14159265358979323846264338328@c@
+#define LOG2 0.69314718055994530943@c@
+#define INVLOG2 1.4426950408889634074@c@
 #define degrees@c@ rad2deg@c@
 #define radians@c@ deg2rad@c@
 
@@ -50,6 +54,30 @@
 }
 
 static @type@
+log2_1p@c@(@type@ x)
+{
+    @type@ u = 1 + x;
+    if (u == 1) {
+        return INVLOG2*x;
+    } else {
+        return log2@c@(u) * x / (u - 1);
+    }
+}
+
+static @type@
+exp2_1m@c@(@type@ x)
+{
+    @type@ u = exp@c@(x);
+    if (u == 1.0) {
+        return LOG2*x;
+    } else if (u - 1 == -1) {
+        return -LOG2;
+    } else {
+        return (u - 1) * x/log2@c@(u);
+    }
+}
+
+static @type@
 logaddexp@c@(@type@ x, @type@ y)
 {
     const @type@ tmp = x - y;
@@ -61,7 +89,21 @@
     }
 }
 
+static @type@
+logaddexp2@c@(@type@ x, @type@ y)
+{
+    const @type@ tmp = x - y;
+    if (tmp > 0) {
+        return x + log2_1p@c@(exp2@c@(-tmp));
+    }
+    else {
+        return y + log2_1p@c@(exp2@c@(tmp));
+    }
+}
+
 #undef PI
+#undef LOG2
+#undef INVLOG2
 
 /**end repeat**/
 

Modified: trunk/numpy/core/tests/test_umath.py
===================================================================
--- trunk/numpy/core/tests/test_umath.py	2008-11-11 03:26:53 UTC (rev 5999)
+++ trunk/numpy/core/tests/test_umath.py	2008-11-11 04:42:25 UTC (rev 6000)
@@ -41,28 +41,31 @@
         y = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
         for dt in ['f','d','g'] :
             xf = np.array(x, dtype=dt)
-            assert_almost_equal(np.log2(xf), y)
+            yf = np.array(y, dtype=dt)
+            assert_almost_equal(np.log2(xf), yf)
 
 class TestExp2(TestCase):
     def test_exp2_values(self) :
         x = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]
         y = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
         for dt in ['f','d','g'] :
+            xf = np.array(x, dtype=dt)
             yf = np.array(y, dtype=dt)
-            assert_almost_equal(np.exp2(yf), x)
+            assert_almost_equal(np.exp2(yf), xf)
 
-class TestLogAddExp(TestCase):
-    def test_logaddexp_values(self) :
+class TestLogAddExp2(TestCase):
+    # Need test for intermediate precisions
+    def test_logaddexp2_values(self) :
         x = [1, 2, 3, 4, 5]
         y = [5, 4, 3, 2, 1]
         z = [6, 6, 6, 6, 6]
-        for dt in ['f','d','g'] :
-            logxf = np.log(np.array(x, dtype=dt))
-            logyf = np.log(np.array(y, dtype=dt))
-            logzf = np.log(np.array(z, dtype=dt))
-            assert_almost_equal(np.logaddexp(logxf, logyf), logzf)
+        for dt, dec in zip(['f','d','g'],[6, 15, 15]) :
+            xf = np.log2(np.array(x, dtype=dt))
+            yf = np.log2(np.array(y, dtype=dt))
+            zf = np.log2(np.array(z, dtype=dt))
+            assert_almost_equal(np.logaddexp2(xf, yf), zf, decimal=dec)
 
-    def test_logaddexp_range(self) :
+    def test_logaddexp2_range(self) :
         x = [1000000, -1000000, 2000000, -2000000]
         y = [2000000, -2000000, 1000000, -1000000]
         z = [2000000, -1000000, 2000000, -1000000]
@@ -72,16 +75,46 @@
             logzf = np.array(z, dtype=dt)
             assert_almost_equal(np.logaddexp(logxf, logyf), logzf)
 
+class TestLog(TestCase):
+    def test_log_values(self) :
+        x = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]
+        y = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+        for dt in ['f','d','g'] :
+            log2_ = 0.69314718055994530943
+            xf = np.array(x, dtype=dt)
+            yf = np.array(y, dtype=dt)*log2_
+            assert_almost_equal(np.log(xf), yf)
+
+class TestExp(TestCase):
+    def test_exp_values(self) :
+        x = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]
+        y = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+        for dt in ['f','d','g'] :
+            log2_ = 0.69314718055994530943
+            xf = np.array(x, dtype=dt)
+            yf = np.array(y, dtype=dt)*log2_
+            assert_almost_equal(np.exp(yf), xf)
+
 class TestLogAddExp(TestCase):
-    def test_logaddexp(self) :
+    def test_logaddexp_values(self) :
         x = [1, 2, 3, 4, 5]
         y = [5, 4, 3, 2, 1]
         z = [6, 6, 6, 6, 6]
+        for dt, dec in zip(['f','d','g'],[6, 15, 15]) :
+            xf = np.log(np.array(x, dtype=dt))
+            yf = np.log(np.array(y, dtype=dt))
+            zf = np.log(np.array(z, dtype=dt))
+            assert_almost_equal(np.logaddexp(xf, yf), zf, decimal=dec)
+
+    def test_logaddexp_range(self) :
+        x = [1000000, -1000000, 2000000, -2000000]
+        y = [2000000, -2000000, 1000000, -1000000]
+        z = [2000000, -1000000, 2000000, -1000000]
         for dt in ['f','d','g'] :
-            logxf = np.log2(np.array(x, dtype=dt))
-            logyf = np.log2(np.array(y, dtype=dt))
-            logzf = np.log2(np.array(z, dtype=dt))
-            assert_almost_equal(np.logaddexp2(logxf, logyf), logzf)
+            logxf = np.array(x, dtype=dt)
+            logyf = np.array(y, dtype=dt)
+            logzf = np.array(z, dtype=dt)
+            assert_almost_equal(np.logaddexp(logxf, logyf), logzf)
 
 class TestLog1p(TestCase):
     def test_log1p(self):



More information about the Numpy-svn mailing list