[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Guile-commits] 05/06: Support general arrays in random:hollow-sphere!
From: |
Daniel Llorens |
Subject: |
[Guile-commits] 05/06: Support general arrays in random:hollow-sphere! |
Date: |
Mon, 11 Sep 2017 07:37:39 -0400 (EDT) |
lloda pushed a commit to branch wip-exception-truncate
in repository guile.
commit 61049cbce337872a50ee2bcacfdf22af6674e45c
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 | 138 ++++++++++++++++++++++++-------------------
test-suite/tests/random.test | 49 ++++++++++++++-
2 files changed, 123 insertions(+), 64 deletions(-)
diff --git a/libguile/random.c b/libguile/random.c
index a8ad075..2992ace 100644
--- a/libguile/random.c
+++ b/libguile/random.c
@@ -501,63 +501,76 @@ 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;
+ if (dims[0].lbnd > ubnd) /* avoid _writable_elements on empty vector */
+ return;
- 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
@@ -613,10 +626,10 @@ SCM_DEFINE (scm_random_normal_vector_x,
"random:normal-vector!", 1, 1, 0,
"(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);
SCM_VALIDATE_RSTATE (2, state);
@@ -628,24 +641,25 @@ SCM_DEFINE (scm_random_normal_vector_x,
"random:normal-vector!", 1, 1, 0,
scm_wrong_type_arg_msg (NULL, 0, v, "rank 1 array");
}
- dim = scm_array_handle_dims (&handle);
-
- if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
+ dims = scm_array_handle_dims (&handle);
+ if (dims->lbnd <= dims->ubnd) /* avoid _writable_elements on empty vector */
{
- 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)));
- }
- 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));
+ if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
+ {
+ SCM *elts = scm_array_handle_writable_elements (&handle);
+ 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 = 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
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))
+
+
- [Guile-commits] branch wip-exception-truncate updated (129d631 -> 4376871), Daniel Llorens, 2017/09/11
- [Guile-commits] 06/06: (wip) give a handle into format used in exceptions, Daniel Llorens, 2017/09/11
- [Guile-commits] 04/06: Remove scm_generalized_vector_get_handle, Daniel Llorens, 2017/09/11
- [Guile-commits] 05/06: Support general arrays in random:hollow-sphere!,
Daniel Llorens <=
- [Guile-commits] 01/06: Fix sort, sort! for arrays with nonzero lower bound, Daniel Llorens, 2017/09/11
- [Guile-commits] 02/06: Support non-zero lower bounds in array-slice-for-each, Daniel Llorens, 2017/09/11
- [Guile-commits] 03/06: Fix bitvectors and non-zero lower bound arrays in truncated-print, Daniel Llorens, 2017/09/11