From 04dfa37f4c55cd45a3d86a1bbae5a874433cc13d Mon Sep 17 00:00:00 2001
From: Sumant Oemrawsingh <soemraws@xs4all.nl>
Date: Wed, 1 Sep 2010 23:11:21 +0200
Subject: [PATCH] Fix for random errors in tests

The test functions did not properly initialize newly created vectors to
contain unique values. This caused comparison of elements that are out of the
stride to randomly succeed. The new test functions properly initialize the
vectors just like is done in the GSL test functions. The result is that
comparison of elements that are out of stride should fail consistently.

This automatically means that the comparison should only check elements which
are on stride, just like is done in GSL. This is not fixed yet, so the
fast-fourier-transform unit tests are still left disabled for now.
---
 fast-fourier-transforms/example.lisp | 94 ++++++++++++++++++++++------
 1 file changed, 75 insertions(+), 19 deletions(-)

diff --git a/fast-fourier-transforms/example.lisp b/fast-fourier-transforms/example.lisp
index f7ccaede..c5b1592f 100644
--- a/fast-fourier-transforms/example.lisp
+++ b/fast-fourier-transforms/example.lisp
@@ -73,11 +73,28 @@
     (setf urand-seed 1)
     (values)))
 
+(defun make-and-init-vector (element-type size &key init-offset)
+  "Make a vector of the given element type and size. If init-offset is given,
+   it is assumed to be a valid number, with which the vector is initialised;
+   each element is set to a unique, predetermined value. See also test.c in
+   GSL's fft directory."
+  (let ((vec (grid:make-foreign-array element-type :dimensions (list size))))
+    (when init-offset
+      (loop for i from 0 below size
+            do
+            (setf (grid:gref vec i)
+                  (if (subtypep element-type 'complex)
+                    (coerce (complex (+ (* 2 i) init-offset)
+                                     (+ (* 2 i) init-offset 1))
+                            element-type)
+                    (coerce (+ i init-offset) element-type))))
+    vec)))
+
 ;; (make-urand-vector '(complex double-float) 5)
-(defun make-urand-vector (element-type dimension &key (stride 1))
+(defun make-urand-vector (element-type dimension &key (stride 1) init-offset)
   "Make a vector with random elements."
-  (let ((vec (grid:make-foreign-array `(complex ,(grid:component-float-type element-type))
-			  :dimensions (list (* stride dimension)))))
+  (let ((vec (make-and-init-vector `(complex ,(grid:component-float-type element-type))
+                                   (* stride dimension) :init-offset init-offset)))
     (loop for i from 0 below (* stride dimension) by stride
        do
        (setf (grid:gref vec i)
@@ -86,17 +103,30 @@
 		 (complex (coerce (urand) element-type)))))
     vec))
 
-(defun realpart-vector (complex-vector)
+(defun realpart-vector (complex-vector &key (stride 1) init-offset)
   "The real vector consisting of the real part of the complex vector."
   (let ((real-vector
-	 (grid:make-foreign-array
-	  (grid:component-float-type (element-type complex-vector))
-	  :dimensions (dimensions complex-vector))))
-    (loop for i below (size complex-vector) do
-	 (setf (grid:gref real-vector i)
-	       (realpart (grid:gref complex-vector i))))
+          (make-and-init-vector
+            (grid:component-float-type (element-type complex-vector))
+            (size complex-vector)
+            :init-offset init-offset)))
+    (loop for i below (size complex-vector) by stride
+          do
+          (setf (grid:gref real-vector i)
+                (realpart (grid:gref complex-vector i))))
     real-vector))
 
+(defun copy-with-stride (vector &key (stride 1) init-offset)
+  "Copy a vector and initialize it."
+  (let ((vec
+          (make-and-init-vector (element-type vector)
+                                (size vector)
+                                :init-offset init-offset)))
+    (loop for i below (size vector) by stride
+          do
+          (setf (grid:gref vec i) (grid:gref vector i)))
+    vec))
+
 (defun size-vector-scalar (vector &key (stride 1))
   "Return the size of a vector while taking the stride into account."
   (coerce (floor (size vector) stride)
@@ -104,35 +134,58 @@
 	      (element-type vector)
 	      'double-float)))
 
+#+nil
 (defun vector/length (vector stride)
   (elt/ vector (size-vector-scalar vector :stride stride)))
 
+(defun vector/length (vector &key (stride 1))
+  (let ((element-type (element-type vector)))
+    (loop with length = (size-vector-scalar vector :stride stride)
+          for i from 0 below (size vector) by stride
+          do
+          (setf (grid:gref vector i) (coerce (/ (grid:gref vector i) length) element-type)))
+  vector))
+
 (defun test-real-fft-noise (vector &key (stride 1))
   "Test forward and inverse FFT for a real vector, and return both results in unpacked form."
   (let* ((forward
-	  (forward-fourier-transform (realpart-vector vector) :stride stride))
+	  (forward-fourier-transform (realpart-vector vector :stride stride :init-offset 0)
+                                     :stride stride))
          (inverse
 	  (forward-fourier-transform (copy forward) :half-complex t :stride stride)))
     (values (unpack forward :unpack-type 'complex :stride stride)
-            (unpack (vector/length inverse stride) :stride stride))))
+            (unpack (vector/length inverse :stride stride) :stride stride))))
 
 (defun test-complex-fft-noise (vector &key (stride 1))
   "Test forward, inverse and backward FFT for a complex vector and return all three results."
-  (let ((forward (forward-fourier-transform (copy vector) :stride stride)))
+  (let ((forward
+          (forward-fourier-transform
+            (copy-with-stride vector :stride stride :init-offset 2000)
+            :stride stride)))
     (values forward
-	    (inverse-fourier-transform (copy forward) :stride stride)
-	    (backward-fourier-transform (copy forward) :stride stride))))
+	    (inverse-fourier-transform (copy-with-stride forward :init-offset 0)
+                                       :stride stride)
+            ;; in backward-fourier-transform, we could use copy instead of ;;
+            ;; copy-with-stride
+            (backward-fourier-transform (copy-with-stride forward :init-offset
+                                                          2000)
+                                        :stride stride))))
 
 (defun test-fft-noise (element-type size &key (stride 1))
-  "A test of real forward and complex forward, revese, and inverse FFT
+  "A test of real forward and complex forward, reverse, and inverse FFT
    on random noise.  Returns the result of the DFT forward Fourier transformation
    and the forward FFT, which should be the same, and the original vector and the
    inverse FFT, which should also be the same.  In addition,
    the backward Fourier transform is returned for complex vectors, which
    should be the same as the last two."
-  (let* ((random-vector (make-urand-vector element-type size :stride stride))
-	 (dft-random-vector
-	  (forward-discrete-fourier-transform random-vector :stride stride)))
+  (let* ((random-vector (make-urand-vector element-type size :stride stride
+                                           :init-offset 1000))
+         (dft-random-vector
+           (forward-discrete-fourier-transform random-vector :stride stride
+                                               :result (make-and-init-vector
+                                                         `(complex ,(grid:component-float-type element-type))
+                                                         (* size stride)
+                                                         :init-offset 3000))))
     (if (subtypep element-type 'complex)
 	(multiple-value-bind (forward inverse backward)
 	    (test-complex-fft-noise random-vector :stride stride)
@@ -155,6 +208,9 @@
 	    (test-real-fft-noise random-vector :stride stride)
 	  (values dft-random-vector forward random-vector inverse)))))
 
+(test-fft-noise 'double-float 10 :stride 2)
+(test-fft-noise '(complex double-float) 10 :stride 2)
+
 ;; (test-fft-noise 'double-float 10 :stride 1)
 ;; (test-fft-noise '(complex double-float) 10 :stride 1)
 
-- 
GitLab