[Guile-commits] 08/08: Support general arrays in random:hollow-sphere!

[Daniel Llorens]

lloda pushed a commit to branch wip-lloda
in repository guile.

commit 6311ea3a384d8fc0d956f7f503ab043a4a2bc4b1
Author: Daniel Llorens <address@hidden>
Date:   Thu Mar 30 14:29:59 2017 +0200

    Support general arrays in random:hollow-sphere!
    
    * libguile/random.c (vector_scale_x, vector_sum_squares): Handle general
      rank-1 #t or 'f64 arrays.
    * test-suite/tests/random.test: Add tests for random:hollow-sphere!.
---
 libguile/random.c            | 161 +++++++++++++++++++++++--------------------
 test-suite/tests/random.test |  49 ++++++++++++-
 2 files changed, 132 insertions(+), 78 deletions(-)

diff --git a/libguile/random.c b/libguile/random.c
index a8ad075..76d214e 100644
--- a/libguile/random.c
+++ b/libguile/random.c
@@ -140,7 +140,7 @@ scm_i_rstate_from_datum (scm_t_rstate *state, SCM value)
   scm_t_i_rstate *istate = (scm_t_i_rstate*) state;
   scm_t_uint32 w, c;
   long length;
-  
+
   SCM_VALIDATE_LIST_COPYLEN (SCM_ARG1, value, length);
   SCM_ASSERT (length == 3, value, SCM_ARG1, FUNC_NAME);
   SCM_ASSERT (scm_is_eq (SCM_CAR (value), scm_i_rstate_tag),
@@ -227,13 +227,13 @@ scm_c_normal01 (scm_t_rstate *state)
   else
     {
       double r, a, n;
-      
+
       r = sqrt (-2.0 * log (scm_c_uniform01 (state)));
       a = 2.0 * M_PI * scm_c_uniform01 (state);
-      
+
       n = r * sin (a);
       state->normal_next = r * cos (a);
-      
+
       return n;
     }
 }
@@ -274,7 +274,7 @@ scm_c_random64 (scm_t_rstate *state, scm_t_uint64 m)
 
   if (m <= SCM_T_UINT32_MAX)
     return scm_c_random (state, (scm_t_uint32) m);
-  
+
   mask = scm_i_mask32 (m >> 32);
   while ((r = ((scm_t_uint64) (state->rng->random_bits (state) & mask) << 32)
           | state->rng->random_bits (state)) >= m)
@@ -322,7 +322,7 @@ scm_c_random_bignum (scm_t_rstate *state, SCM m)
       scm_t_uint32 chunks_left = num_chunks;
 
       mpz_set_ui (SCM_I_BIG_MPZ (result), 0);
-      
+
       if (end_bits)
         {
           /* generate a mask with ones in the end_bits position, i.e. if
@@ -334,7 +334,7 @@ scm_c_random_bignum (scm_t_rstate *state, SCM m)
           *current_chunk-- = highest_bits;
           chunks_left--;
         }
-      
+
       while (chunks_left)
         {
           /* now fill in the remaining scm_t_uint32 sized chunks */
@@ -360,7 +360,7 @@ scm_c_random_bignum (scm_t_rstate *state, SCM m)
 /*
  * Scheme level representation of random states.
  */
- 
+
 scm_t_bits scm_tc16_rstate;
 
 static SCM
@@ -376,7 +376,7 @@ make_rstate (scm_t_rstate *state)
 
 SCM_GLOBAL_VARIABLE_INIT (scm_var_random_state, "*random-state*", 
scm_seed_to_random_state (scm_from_locale_string 
("URL:http://stat.fsu.edu/~geo/diehard.html";)));
 
-SCM_DEFINE (scm_random, "random", 1, 1, 0, 
+SCM_DEFINE (scm_random, "random", 1, 1, 0,
             (SCM n, SCM state),
             "Return a number in [0, N).\n"
             "\n"
@@ -420,7 +420,7 @@ SCM_DEFINE (scm_random, "random", 1, 1, 0,
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, 
+SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
             (SCM state),
             "Return a copy of the random state @var{state}.")
 #define FUNC_NAME s_scm_copy_random_state
@@ -432,7 +432,7 @@ SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 
1, 0,
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, 
+SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
             (SCM seed),
             "Return a new random state using @var{seed}.")
 #define FUNC_NAME s_scm_seed_to_random_state
@@ -445,11 +445,11 @@ SCM_DEFINE (scm_seed_to_random_state, 
"seed->random-state", 1, 0, 0,
                                        scm_i_string_length (seed)));
   scm_remember_upto_here_1 (seed);
   return res;
-  
+
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0, 
+SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0,
             (SCM datum),
             "Return a new random state using @var{datum}, which should have\n"
             "been obtained from @code{random-state->datum}.")
@@ -459,7 +459,7 @@ SCM_DEFINE (scm_datum_to_random_state, 
"datum->random-state", 1, 0, 0,
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0, 
+SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0,
             (SCM state),
             "Return a datum representation of @var{state} that may be\n"
             "written out and read back with the Scheme reader.")
@@ -470,7 +470,7 @@ SCM_DEFINE (scm_random_state_to_datum, 
"random-state->datum", 1, 0, 0,
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, 
+SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
             (SCM state),
            "Return a uniformly distributed inexact real random number in\n"
            "[0,1).")
@@ -483,7 +483,7 @@ SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, 
+SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
             (SCM state),
            "Return an inexact real in a normal distribution.  The\n"
            "distribution used has mean 0 and standard deviation 1.  For a\n"
@@ -501,63 +501,73 @@ SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
 static void
 vector_scale_x (SCM v, double c)
 {
-  size_t n;
-  if (scm_is_vector (v))
-    {
-      n = SCM_SIMPLE_VECTOR_LENGTH (v);
-      while (n-- > 0)
-       SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)) *= c;
-    }
-  else
-    {
-      /* must be a f64vector. */
-      scm_t_array_handle handle;
-      size_t i, len;
-      ssize_t inc;
-      double *elts;
+  scm_t_array_handle handle;
+  scm_t_array_dim const * dims;
+  ssize_t i, inc, ubnd;
 
-      elts = scm_f64vector_writable_elements (v, &handle, &len, &inc);
+  scm_array_get_handle (v, &handle);
+  dims = scm_array_handle_dims (&handle);
+  if (1 == scm_array_handle_rank (&handle))
+    {
+      ubnd = dims[0].ubnd;
+      inc = dims[0].inc;
 
-      for (i = 0; i < len; i++, elts += inc)
-       *elts *= c;
-      
-      scm_array_handle_release (&handle);
+      if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_F64)
+        {
+          double *elts = (double *)(handle.writable_elements) + handle.base;
+          for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc)
+            *elts *= c;
+          return;
+        }
+      else if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
+        {
+          SCM *elts = (SCM *)(handle.writable_elements) + handle.base;
+          for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc)
+            SCM_REAL_VALUE (*elts) *= c;
+          return;
+        }
     }
+  scm_array_handle_release (&handle);
+  scm_misc_error (NULL, "must be a rank-1 array of type #t or 'f64", 
scm_list_1 (v));
 }
 
 static double
 vector_sum_squares (SCM v)
 {
   double x, sum = 0.0;
-  size_t n;
-  if (scm_is_vector (v))
-    {
-      n = SCM_SIMPLE_VECTOR_LENGTH (v);
-      while (n-- > 0)
-       {
-         x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n));
-         sum += x * x;
-       }
-    }
-  else
-    {
-      /* must be a f64vector. */
-      scm_t_array_handle handle;
-      size_t i, len;
-      ssize_t inc;
-      const double *elts;
-
-      elts = scm_f64vector_elements (v, &handle, &len, &inc);
-
-      for (i = 0; i < len; i++, elts += inc)
-       {
-         x = *elts;
-         sum += x * x;
-       }
+  scm_t_array_handle handle;
+  scm_t_array_dim const * dims;
+  ssize_t i, inc, ubnd;
 
-      scm_array_handle_release (&handle);
+  scm_array_get_handle (v, &handle);
+  dims = scm_array_handle_dims (&handle);
+  if (1 == scm_array_handle_rank (&handle))
+    {
+      ubnd = dims[0].ubnd;
+      inc = dims[0].inc;
+      if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_F64)
+        {
+          const double *elts = (const double *)(handle.elements) + handle.base;
+          for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc)
+            {
+              x = *elts;
+              sum += x * x;
+            }
+          return sum;
+        }
+      else if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
+        {
+          const SCM *elts = (const SCM *)(handle.elements) + handle.base;
+          for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc)
+            {
+              x = SCM_REAL_VALUE (*elts);
+              sum += x * x;
+            }
+          return sum;
+        }
     }
-  return sum;
+  scm_array_handle_release (&handle);
+  scm_misc_error (NULL, "must be an array of type #t or 'f64", scm_list_1 (v));
 }
 
 /* For the uniform distribution on the solid sphere, note that in
@@ -565,7 +575,7 @@ vector_sum_squares (SCM v)
  * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
  * generated as r=u^(1/n).
  */
-SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, 
+SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
             (SCM v, SCM state),
            "Fills @var{vect} with inexact real random numbers the sum of\n"
            "whose squares is less than 1.0.  Thinking of @var{vect} as\n"
@@ -606,16 +616,16 @@ SCM_DEFINE (scm_random_hollow_sphere_x, 
"random:hollow-sphere!", 1, 1, 0,
 #undef FUNC_NAME
 
 
-SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, 
+SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0,
             (SCM v, SCM state),
             "Fills vect with inexact real random numbers that are\n"
             "independent and standard normally distributed\n"
             "(i.e., with mean 0 and variance 1).")
 #define FUNC_NAME s_scm_random_normal_vector_x
 {
-  long i;
   scm_t_array_handle handle;
-  scm_t_array_dim *dim;
+  scm_t_array_dim const * dims;
+  ssize_t i;
 
   if (SCM_UNBNDP (state))
     state = SCM_VARIABLE_REF (scm_var_random_state);
@@ -627,30 +637,29 @@ SCM_DEFINE (scm_random_normal_vector_x, 
"random:normal-vector!", 1, 1, 0,
       scm_array_handle_release (&handle);
       scm_wrong_type_arg_msg (NULL, 0, v, "rank 1 array");
     }
-  
-  dim = scm_array_handle_dims (&handle);
+
+  dims = scm_array_handle_dims (&handle);
 
   if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
     {
       SCM *elts = scm_array_handle_writable_elements (&handle);
-      for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc)
-       *elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
+      for (i = dims->lbnd; i <= dims->ubnd; i++, elts += dims->inc)
+        *elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
     }
   else
     {
       /* must be a f64vector. */
       double *elts = scm_array_handle_f64_writable_elements (&handle);
-      for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc)
-       *elts = scm_c_normal01 (SCM_RSTATE (state));
+      for (i = dims->lbnd; i <= dims->ubnd; i++, elts += dims->inc)
+        *elts = scm_c_normal01 (SCM_RSTATE (state));
     }
 
   scm_array_handle_release (&handle);
-
   return SCM_UNSPECIFIED;
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, 
+SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
             (SCM state),
            "Return an inexact real in an exponential distribution with mean\n"
            "1.  For an exponential distribution with mean u use (* u\n"
@@ -781,13 +790,13 @@ scm_init_random ()
     scm_i_rstate_to_datum
   };
   scm_the_rng = rng;
-  
+
   scm_tc16_rstate = scm_make_smob_type ("random-state", 0);
 
   for (m = 1; m <= 0x100; m <<= 1)
     for (i = m >> 1; i < m; ++i)
       scm_masktab[i] = m - 1;
-  
+
 #include "libguile/random.x"
 
   scm_add_feature ("random");
diff --git a/test-suite/tests/random.test b/test-suite/tests/random.test
index ab20b58..b5c9692 100644
--- a/test-suite/tests/random.test
+++ b/test-suite/tests/random.test
@@ -20,7 +20,8 @@
   #:use-module ((system base compile) #:select (compile))
   #:use-module (test-suite lib)
   #:use-module (srfi srfi-4)
-  #:use-module (srfi srfi-4 gnu))
+  #:use-module (srfi srfi-4 gnu)
+  #:use-module ((ice-9 control) #:select (let/ec)))
 
 ; see strings.test, arrays.test.
 (define exception:wrong-type-arg
@@ -52,4 +53,48 @@
       (begin
         (random:normal-vector! b (random-state-from-platform))
         (random:normal-vector! c (random-state-from-platform))
-        (and (not (equal? a b)) (not (equal? a c)))))))
+        (and (not (equal? a b)) (not (equal? a c))))))
+
+  (pass-if "empty argument"
+    (random:normal-vector! (vector) (random-state-from-platform))
+    (random:normal-vector! (f64vector) (random-state-from-platform))
+    #t))
+
+;;;
+;;; random:hollow-sphere!
+;;;
+
+(with-test-prefix "random:hollow-sphere!"
+
+  (define (sqr a)
+    (* a a))
+  (define (norm a)
+    (sqrt (+ (sqr (array-ref a 0)) (sqr (array-ref a 1)) (sqr (array-ref a 
2)))))
+  (define double-eps 1e-15)
+  
+  (pass-if "non uniform"
+    (let ((a (transpose-array (make-array 0. 3 10) 1 0)))
+      (let/ec exit
+        (array-slice-for-each 1
+          (lambda (a)
+            (random:hollow-sphere! a)
+            (if (> (magnitude (- 1 (norm a))) double-eps) (exit #f)))
+          a)
+        #t)))
+
+  (pass-if "uniform (f64)"
+    (let ((a (transpose-array (make-array 0. 3 10) 1 0)))
+      (let/ec exit
+        (array-slice-for-each 1
+          (lambda (a)
+            (random:hollow-sphere! a)
+            (if (> (magnitude (- 1 (norm a))) double-eps) (exit #f)))
+          a)
+        #t)))
+
+  (pass-if "empty argument"
+    (random:hollow-sphere! (vector))
+    (random:hollow-sphere! (f64vector))
+    #t))
+
+