r15348 - in /branches/frame_order_testing/extern/sobol: test.log test.py -- February 21, 2012 - 14:17

Author: bugman
Date: Tue Feb 21 14:17:03 2012
New Revision: 15348

URL: http://svn.gna.org/viewcvs/relax?rev=15348&view=rev
Log:
Created a simple script to test for bias in the Sobol' sequence 
implementation.


Added:
    branches/frame_order_testing/extern/sobol/test.log
    branches/frame_order_testing/extern/sobol/test.py

Added: branches/frame_order_testing/extern/sobol/test.log
URL: 
http://svn.gna.org/viewcvs/relax/branches/frame_order_testing/extern/sobol/test.log?rev=15348&view=auto
==============================================================================
--- branches/frame_order_testing/extern/sobol/test.log (added)
+++ branches/frame_order_testing/extern/sobol/test.log Tue Feb 21 14:17:03 
2012
@@ -1,0 +1,21 @@
+
+N = 1
+Average vector length: 0.866025403784
+
+N = 10
+Average vector length: 0.0819679815538
+
+N = 100
+Average vector length: 0.0082620650415
+
+N = 1000
+Average vector length: 0.000446544348826
+
+N = 10000
+Average vector length: 0.000127752105696
+
+N = 100000
+Average vector length: 7.09521003856e-06
+
+N = 1000000
+Average vector length: 1.32445797827e-06

Added: branches/frame_order_testing/extern/sobol/test.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/frame_order_testing/extern/sobol/test.py?rev=15348&view=auto
==============================================================================
--- branches/frame_order_testing/extern/sobol/test.py (added)
+++ branches/frame_order_testing/extern/sobol/test.py Tue Feb 21 14:17:03 2012
@@ -1,0 +1,35 @@
+# Python module imports.
+from numpy import float64, ones, zeros
+from numpy.linalg import norm
+
+# relax module imports.
+from sobol_lib import i4_sobol
+
+
+# Some variables.
+DIM = 3
+OFFSET = 0.5 * ones(DIM)
+
+# Loop over different number of points.
+for exponent in range(7): 
+    # The number of points.
+    N = int(10**exponent)
+
+    # Initialise a vector.
+    ave_pos = zeros(DIM, float64)
+
+    # Print out.
+    print("\nN = %s" % N)
+
+    # Loop over the points.
+    for i in range(N):
+        # The raw point.
+        point, seed = i4_sobol(DIM, i)
+
+        # Sum the point, minus the offset.
+        ave_pos += point - OFFSET
+
+    # The average vector length.
+    ave_pos = ave_pos / float(N)
+    r = norm(ave_pos)
+    print("Average vector length: %s" % r)