1 ;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp -*-
4 Floating tests and assertions for LISP-UNIT
6 Copyright
(c) 2009-
2012, Thomas M. Hermann
8 Permission is hereby granted
, free of charge
, to any person obtaining
9 a copy of this software and associated documentation files
(the "Software"),
10 to deal in the Software without restriction
, including without limitation
11 the rights to use
, copy
, modify
, merge
, publish
, distribute
, sublicense
,
12 and
/or sell copies of the Software
, and to permit persons to whom the
13 Software is furnished to do so
, subject to the following conditions
:
15 The above copyright notice and this permission notice shall be included
16 in all copies or substantial portions of the Software.
18 THE SOFTWARE IS PROVIDED
"AS IS", WITHOUT WARRANTY OF ANY KIND
, EXPRESS
19 OR IMPLIED
, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY
,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
21 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM
, DAMAGES OR
22 OTHER LIABILITY
, WHETHER IN AN ACTION OF CONTRACT
, TORT OR OTHERWISE
,
23 ARISING FROM
, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24 OTHER DEALINGS IN THE SOFTWARE.
27 [NumAlgoC] Gisela Engeln-Mullges and Frank Uhlig "Numerical
28 Algorithms with C", Springer, 1996
33 (in-package :lisp-unit)
35 ;;; Symbols exported from the floating point extension
39 '(*measure* *epsilon* *significant-figures*))
45 relative-error relative-error-norm
48 ;;; Predicates and assertions
50 '(float-equal assert-float-equal
51 sigfig-equal assert-sigfig-equal
52 norm-equal assert-norm-equal
53 number-equal assert-number-equal
54 numerical-equal assert-numerical-equal))
62 make-random-2d-array))
64 ;;; Floating point extensions
69 "Set the error epsilon if the defaults are not acceptable.")
71 (defvar *significant-figures* 4
72 "Default to 4 significant figures.")
74 (defgeneric default-epsilon (value)
76 "Return the default epsilon for the value."))
78 (defgeneric relative-error (exact approximate)
80 "Return the relative-error between the 2 quantities."))
82 (defgeneric float-equal (data1 data2 &optional epsilon)
84 "Return true if the floating point data is equal."))
86 (defgeneric sumsq (data)
88 "Return the scaling parameter and the sum of the squares of the ~
91 (defgeneric sump (data p)
93 "Return the scaling parameter and the sum of the powers of p of the ~
96 (defgeneric norm (data &optional measure)
98 "Return the element-wise norm of the data."))
100 (defgeneric relative-error-norm (exact approximate &optional measure)
102 "Return the relative error norm "))
104 (defgeneric norm-equal (data1 data2 &optional epsilon measure)
106 "Return true if the norm of the data is equal."))
108 (defgeneric sigfig-equal (data1 data2 &optional significant-figures)
110 "Return true if the data have equal significant figures."))
112 (defgeneric numerical-equal (result1 result2 &key test)
114 "Return true if the results are numerically equal according to :TEST."))
116 ;;; (DEFAULT-EPSILON value) => epsilon
117 (defmethod default-epsilon ((value float))
118 "Return a default epsilon value based on the floating point type."
120 (short-float (* 2S0 short-float-epsilon))
121 (single-float (* 2F0 single-float-epsilon))
122 (double-float (* 2D0 double-float-epsilon))
123 (long-float (* 2L0 long-float-epsilon))))
125 (defmethod default-epsilon ((value complex))
126 "Return a default epsilon value based on the complex type."
128 ((complex short-float) (* 2S0 short-float-epsilon))
129 ((complex single-float) (* 2F0 single-float-epsilon))
130 ((complex double-float) (* 2D0 double-float-epsilon))
131 ((complex long-float) (* 2L0 long-float-epsilon))
134 (defmethod default-epsilon ((value list))
135 "Return the default epsilon based on contents of the list."
136 (loop for val in value maximize (default-epsilon val)))
138 (defmethod default-epsilon ((value vector))
139 "Return the default epsilon based on the contents of the vector."
140 (loop for val across value maximize (default-epsilon val)))
142 (defmethod default-epsilon ((value array))
143 "Return the default epsilon based on the contents of the array."
146 (array-total-size value)
147 :element-type (array-element-type value)
149 maximize (default-epsilon val)))
152 (RELATIVE-ERROR x y) => float
153 [NumAlgoC] : Definition
1.3, pg.
2
154 modified with Definition
1.1, pg.
1
156 The definition of relative error in this routine is modified from
157 the Definition
1.3 in
[NumAlgoC] for cases when either the exact
158 or the approximate value equals zero. According to Definition 1.3,
159 the relative error is identically equal to 1 in those cases. This
160 function returns the absolue error in those cases. This is more
163 (defun %relative-error (exact approximate)
164 "Return the relative error of the numbers."
165 (abs (if (or (zerop exact) (zerop approximate))
166 (- exact approximate)
167 (/ (- exact approximate) exact))))
169 (defmethod relative-error ((exact float) (approximate float))
170 "Return the error delta between the exact and approximate floating
172 (%relative-error exact approximate))
174 (defmethod relative-error ((exact float) (approximate complex))
175 "Return the relative error between the float and complex number."
176 (%relative-error exact approximate))
178 (defmethod relative-error ((exact complex) (approximate float))
179 "Return the relative error between the float and complex number."
180 (%relative-error exact approximate))
182 (defmethod relative-error ((exact complex) (approximate complex))
183 "Return the relative error of the complex numbers."
184 (if (or (typep exact '(complex float))
185 (typep approximate '(complex float)))
186 (%relative-error exact approximate)
187 (error "Relative error is only applicable to complex values with ~
188 floating point parts.")))
190 ;;; (FLOAT-EQUAL data1 data2 epsilon) => true or false
191 (defun %float-equal (data1 data2 epsilon)
192 "Return true if the relative error between the data is less than
195 (and (zerop data1) (zerop data2))
196 (< (%relative-error data1 data2) epsilon)))
198 (defmethod float-equal ((data1 float) (data2 float)
199 &optional (epsilon *epsilon*))
200 "Return true if the relative error between data1 and data2 is less
202 (%float-equal data1 data2
203 (or epsilon (max (default-epsilon data1)
204 (default-epsilon data2)))))
206 (defmethod float-equal ((data1 float) (data2 rational)
207 &optional (epsilon *epsilon*))
208 "Return true if the relative error between data1 and data2 is less
210 (%float-equal data1 (float data2 data1)
211 (or epsilon (default-epsilon data1))))
213 (defmethod float-equal ((data1 rational) (data2 float)
214 &optional (epsilon *epsilon*))
215 "Return true if the relative error between data1 and data2 is less
217 (%float-equal (float data1 data2) data2
218 (or epsilon (default-epsilon data2))))
220 (defmethod float-equal ((data1 float) (data2 complex)
221 &optional (epsilon *epsilon*))
222 "Return true if the relative error between data1 and data2 is less
224 (%float-equal data1 data2
225 (or epsilon (max (default-epsilon data1)
226 (default-epsilon data2)))))
228 (defmethod float-equal ((data1 complex) (data2 float)
229 &optional (epsilon *epsilon*))
230 "Return true if the relative error between data1 and data2 is less
232 (%float-equal data1 data2
233 (or epsilon (max (default-epsilon data1)
234 (default-epsilon data2)))))
236 (defmethod float-equal ((data1 complex) (data2 complex)
237 &optional (epsilon *epsilon*))
238 "Return true if the relative error between data1 and data2 is less
240 (< (relative-error data1 data2)
241 (or epsilon (max (default-epsilon data1)
242 (default-epsilon data2)))))
244 (defun %seq-float-equal (seq1 seq2 epsilon)
245 "Return true if the element-wise comparison of relative error is
248 (and (null seq1) (null seq2))
249 (when (= (length seq1) (length seq2))
251 (lambda (d1 d2) (float-equal d1 d2 epsilon)) seq1 seq2))))
253 (defmethod float-equal ((data1 list) (data2 list)
254 &optional (epsilon *epsilon*))
255 "Return true if the lists are equal in length and element-wise
256 comparison of the relative error is less than epsilon."
257 (%seq-float-equal data1 data2 epsilon))
259 (defmethod float-equal ((data1 list) (data2 vector)
260 &optional (epsilon *epsilon*))
261 "Return true if the vector and the list are equal in length and
262 element-wise comparison of the relative error is less than epsilon."
263 (%seq-float-equal data1 data2 epsilon))
265 (defmethod float-equal ((data1 vector) (data2 list)
266 &optional (epsilon *epsilon*))
267 "Return true if the vector and the list are equal in length and
268 element-wise comparison of the relative error is less than epsilon."
269 (%seq-float-equal data1 data2 epsilon))
271 (defmethod float-equal ((data1 vector) (data2 vector)
272 &optional (epsilon *epsilon*))
273 "Return true if the vectors are equal in length and element-wise
274 comparison of the relative error is less than epsilon."
275 (%seq-float-equal data1 data2 epsilon))
277 (defmethod float-equal ((data1 array) (data2 array)
278 &optional (epsilon *epsilon*))
279 "Return true if the arrays are equal in length and element-wise
280 comparison of the relative error is less than epsilon."
281 (when (equal (array-dimensions data1)
282 (array-dimensions data2))
284 (make-array (array-total-size data1)
285 :element-type (array-element-type data1)
287 (make-array (array-total-size data2)
288 :element-type (array-element-type data2)
292 (defmacro assert-float-equal (expected form &rest extras)
293 `(expand-assert :equal ,form ,form ,expected ,extras :test #'float-equal))
295 ;;; (SUMSQ data) => scale, sumsq
296 (defmethod sumsq ((data list))
297 "Return the scaling parameter and the sum of the squares of the ~
302 (dolist (elm data (values scale sumsq))
303 (when (plusp (setq abs-val (abs elm)))
304 (if (< scale abs-val)
306 sumsq (1+ (* sumsq (expt (/ scale abs-val) 2)))
308 (setq sumsq (+ sumsq (expt (/ elm scale) 2))))))))
310 (defmethod sumsq ((data vector))
311 "Return the scaling parameter and the sum of the squares of the ~
316 (size (length data)))
317 (dotimes (index size (values scale sumsq))
318 (when (plusp (setq abs-val (abs (svref data index))))
319 (if (< scale abs-val)
321 sumsq (1+ (* sumsq (expt (/ scale abs-val) 2)))
325 (+ sumsq (expt (/ (svref data index) scale) 2))))))))
327 (defmethod sumsq ((data array))
328 "Return the scaling parameter and the sum of the squares of the ~
330 (sumsq (make-array (array-total-size data)
331 :element-type (array-element-type data)
332 :displaced-to data)))
334 ;;; (SUMP data) => scale, sump
335 (defmethod sump ((data list) (p real))
336 "Return the scaling parameter and the sum of the powers of p of the ~
341 (dolist (elm data (values scale sump))
342 (when (plusp (setq abs-val (abs elm)))
343 (if (< scale abs-val)
345 sump (1+ (* sump (expt (/ scale abs-val) p)))
347 (setq sump (+ sump (expt (/ elm scale) p))))))))
349 (defmethod sump ((data vector) (p real))
350 "Return the scaling parameter and the sum of the powers of p of the ~
355 (size (length data)))
356 (dotimes (index size (values scale sump))
357 (when (plusp (setq abs-val (abs (svref data index))))
358 (if (< scale abs-val)
360 sump (1+ (* sump (expt (/ scale abs-val) p)))
364 (+ sump (expt (/ (svref data index) scale) p))))))))
366 (defmethod sump ((data array) (p real))
367 "Return the scaling parameter and the sum of the powers of p of the ~
369 (sump (make-array (array-total-size data)
370 :element-type (array-element-type data)
374 ;;; (NORM data) => float
376 (defgeneric %norm (data measure)
378 "Return the norm of the data according to measure."))
380 (defmethod %norm ((data list) (measure (eql 1)))
381 "Return the Taxicab norm of the list."
382 (loop for item in data sum (abs item)))
384 (defmethod %norm ((data vector) (measure (eql 1)))
385 "Return the Taxicab norm of the vector."
386 (loop for item across data sum (abs item)))
388 (defmethod %norm ((data list) (measure (eql 2)))
389 "Return the Euclidean norm of the list."
390 (multiple-value-bind (scale sumsq)
391 (sumsq (map-into (make-array (length data)) #'abs data))
392 (* scale (sqrt sumsq))))
394 (defmethod %norm ((data vector) (measure (eql 2)))
395 "Return the Euclidean norm of the vector."
396 (multiple-value-bind (scale sumsq)
397 (sumsq (map-into (make-array (length data)) #'abs data))
398 (* scale (sqrt sumsq))))
400 (defmethod %norm ((data list) (measure integer))
401 "Return the Euclidean norm of the list."
402 (multiple-value-bind (scale sump)
403 (sump (map-into (make-array (length data)) #'abs data)
405 (* scale (expt sump (/ measure)))))
407 (defmethod %norm ((data vector) (measure integer))
408 "Return the Euclidean norm of the vector."
409 (multiple-value-bind (scale sump)
410 (sump (map-into (make-array (length data)) #'abs data)
412 (* scale (expt sump (/ measure)))))
414 (defmethod %norm ((data list) (measure (eql :infinity)))
415 "Return the infinity, or maximum, norm of the list."
416 (loop for item in data maximize (abs item)))
418 (defmethod %norm ((data vector) (measure (eql :infinity)))
419 "Return the infinity, or maximum, norm of the vector."
420 (loop for item across data maximize (abs item)))
422 (defmethod norm ((data list) &optional (measure *measure*))
423 "Return the norm of the list according to the measure."
424 (%norm data measure))
426 (defmethod norm ((data vector) &optional (measure *measure*))
427 "Return the norm of the vector according to the measure."
428 (%norm data measure))
430 (defmethod norm ((data array) &optional (measure *measure*))
431 "Return the entrywise norm of the array according to the measure."
434 (array-total-size data)
435 :element-type (array-element-type data)
439 ;;; (RELATIVE-ERROR-NORM exact approximate measure) => float
440 (defun %relative-error-norm (exact approximate measure)
441 "Return the relative error norm of the sequences."
442 (/ (norm (map-into (make-array (length exact))
443 (lambda (x1 x2) (abs (- x1 x2)))
444 exact approximate) measure)
445 (norm exact measure)))
447 (defmethod relative-error-norm ((exact list) (approximate list)
448 &optional (measure *measure*))
449 "Return the relative error norm of the lists."
450 (if (= (length exact) (length approximate))
451 (%relative-error-norm exact approximate measure)
452 (error "Lists are not equal in length.")))
454 (defmethod relative-error-norm ((exact list) (approximate vector)
455 &optional (measure *measure*))
456 "Return the relative error norm of the list and the vector."
457 (if (= (length exact) (length approximate))
458 (%relative-error-norm exact approximate measure)
459 (error "The list and vector are not equal in length.")))
461 (defmethod relative-error-norm ((exact vector) (approximate list)
462 &optional (measure *measure*))
463 "Return the relative error norm of the list and the vector."
464 (if (= (length exact) (length approximate))
465 (%relative-error-norm exact approximate measure)
466 (error "The list and vector are not equal in length.")))
468 (defmethod relative-error-norm ((exact vector) (approximate vector)
469 &optional (measure *measure*))
470 "Return the relative error norm of the vectors."
471 (if (= (length exact) (length approximate))
472 (%relative-error-norm exact approximate measure)
473 (error "Vectors are not equal in length.")))
475 (defmethod relative-error-norm ((exact array) (approximate vector)
476 &optional (measure *measure*))
477 "Return the relative error norm of the arrays."
478 (if (equal (array-dimensions exact)
479 (array-dimensions approximate))
480 (%relative-error-norm
481 (make-array (array-total-size exact)
482 :element-type (array-element-type exact)
484 (make-array (array-total-size approximate)
485 :element-type (array-element-type approximate)
486 :displaced-to approximate)
488 (error "Arrays are not equal dimensions.")))
490 ;;; (NORM-EQUAL data1 data2 epsilon measure) => boolean
491 (defun %norm-equal (seq1 seq2 epsilon measure)
492 "Return true if the relative error norm is less than epsilon."
494 (and (null seq1) (null seq2))
495 (< (%relative-error-norm seq1 seq2 measure) epsilon)))
497 (defmethod norm-equal ((data1 list) (data2 list) &optional
498 (epsilon *epsilon*) (measure *measure*))
499 "Return true if the lists are equal in length and the relative error
500 norm is less than epsilon."
501 (%norm-equal data1 data2 epsilon measure))
503 (defmethod norm-equal ((data1 list) (data2 vector) &optional
504 (epsilon *epsilon*) (measure *measure*))
505 "Return true if the vector and the list are equal in length and the
506 relative error norm is less than epsilon."
507 (%norm-equal data1 data2 epsilon measure))
509 (defmethod norm-equal ((data1 vector) (data2 list) &optional
510 (epsilon *epsilon*) (measure *measure*))
511 "Return true if the vector and the list are equal in length and the
512 relative error norm is less than epsilon."
513 (%norm-equal data1 data2 epsilon measure))
515 (defmethod norm-equal ((data1 vector) (data2 vector) &optional
516 (epsilon *epsilon*) (measure *measure*))
517 "Return true if the vectors are equal in length and the relative
518 error norm is less than epsilon."
519 (%norm-equal data1 data2 epsilon measure))
521 (defmethod norm-equal ((data1 array) (data2 array) &optional
522 (epsilon *epsilon*) (measure *measure*))
523 "Return true if the arrays are equal in length and the relative
524 error norm is less than epsilon."
525 (when (equal (array-dimensions data1)
526 (array-dimensions data2))
528 (make-array (array-total-size data1)
529 :element-type (array-element-type data1)
531 (make-array (array-total-size data2)
532 :element-type (array-element-type data2)
536 (defmacro assert-norm-equal (expected form &rest extras)
537 `(expand-assert :equal ,form ,form ,expected ,extras :test #'norm-equal))
539 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
540 ;;; [NumAlgoC] : Definition
1.7, pg.
4
542 ;;; To avoid using 0.1, first 1.0 <= significand < 10. On the final
543 ;;; return, scale 0.1 <= significand < 1.
544 (defun %normalize-float
(significand &optional
(exponent 0))
545 "Return the normalized floating point number and exponent."
546 ;;; FIXME : Use the LOOP.
549 (values significand
0))
550 ((>= (abs significand
) 10)
551 (%normalize-float
(/ significand
10.0) (1+ exponent
)))
552 ((< (abs significand
) 1)
553 (%normalize-float
(* significand
10.0) (1- exponent
)))
554 (t (values (/ significand
10.0) (1+ exponent
)))))
556 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
557 (defun %sigfig-equal
(float1 float2 significant-figures
)
558 "Return true if the floating point numbers have equal significant
560 (if (or (zerop float1
) (zerop float2
))
561 (< (abs (+ float1 float2
)) (* 5D-1
(expt 1D1
(- significant-figures
))))
562 (multiple-value-bind (sig1 exp1
) (%normalize-float float1
)
563 (multiple-value-bind (sig2 exp2
) (%normalize-float float2
)
564 (= (round (* sig1
(expt 1D1 significant-figures
)))
565 (round (* sig2
(expt 1D1
(- significant-figures
(- exp1 exp2
))))))))))
567 (defmethod sigfig-equal ((data1 float
) (data2 float
) &optional
568 (significant-figures *significant-figures
*))
569 "Return true if the floating point numbers have equal significant
571 (%sigfig-equal data1 data2 significant-figures
))
573 (defun %seq-sigfig-equal
(seq1 seq2 significant-figures
)
574 "Return true if the element-wise comparison is equal to the
575 specified significant figures."
577 (and (null seq1
) (null seq2
))
578 (when (= (length seq1
) (length seq2
))
580 (lambda (d1 d2
) (sigfig-equal d1 d2 significant-figures
))
583 (defmethod sigfig-equal ((data1 list
) (data2 list
) &optional
584 (significant-figures *significant-figures
*))
585 "Return true if the lists are equal in length and the element-wise
586 comparison is equal to significant figures."
587 (%seq-sigfig-equal data1 data2 significant-figures
))
589 (defmethod sigfig-equal ((data1 vector
) (data2 list
) &optional
590 (significant-figures *significant-figures
*))
591 "Return true if the vector and the list are equal in length and the
592 element-wise comparison is equal to significant figures."
593 (%seq-sigfig-equal data1 data2 significant-figures
))
595 (defmethod sigfig-equal ((data1 list
) (data2 vector
) &optional
596 (significant-figures *significant-figures
*))
597 "Return true if the list and the vector are equal in length and the
598 element-wise comparison is equal to significant figures."
599 (%seq-sigfig-equal data1 data2 significant-figures
))
601 (defmethod sigfig-equal ((data1 vector
) (data2 vector
) &optional
602 (significant-figures *significant-figures
*))
603 "Return true if the vectors are equal in length and the element-wise
604 comparison is equal to significant figures."
605 (%seq-sigfig-equal data1 data2 significant-figures
))
607 (defmethod sigfig-equal ((data1 array
) (data2 array
) &optional
608 (significant-figures *significant-figures
*))
609 "Return true if the arrays are equal in length and the element-wise
610 comparison is equal to significant figures."
611 (when (equal (array-dimensions data1
)
612 (array-dimensions data2
))
614 (make-array (array-total-size data1
)
615 :element-type
(array-element-type data1
)
617 (make-array (array-total-size data2
)
618 :element-type
(array-element-type data2
)
620 significant-figures
)))
622 (defmacro assert-sigfig-equal
(expected form
&rest extras
)
623 `(expand-assert :equal
,form
,form
,expected
,extras
:test
#'sigfig-equal
))
625 ;;; (NUMBER-EQUAL number1 number2) => true or false
626 (defun number-equal (number1 number2
&optional
(epsilon *epsilon
*) type-eq-p
)
627 "Return true if the numbers are equal within some epsilon,
628 optionally requiring the types to be identical."
630 (or (not type-eq-p
) (eq (type-of number1
) (type-of number2
)))
631 (float-equal (coerce number1
'(complex double-float
))
632 (coerce number2
'(complex double-float
))
635 (defmacro assert-number-equal
(expected form
&rest extras
)
636 `(expand-assert :equal
,form
,form
,expected
,extras
:test
#'number-equal
))
638 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
640 ;;; This is a universal wrapper created by Liam Healy. It is
641 ;;; implemented to support testing in GSLL. The interface is expanded,
642 ;;; but backwards compatible with previous versions.
644 (defmethod numerical-equal ((result1 number
) (result2 number
)
645 &key
(test #'number-equal
))
646 "Return true if the the numbers are equal according to :TEST."
647 (funcall test result1 result2
))
649 (defun %sequence-equal
(seq1 seq2 test
)
650 "Return true if the sequences are equal in length and each element
651 is equal according to :TEST."
652 (when (= (length seq1
) (length seq2
))
653 (every (lambda (s1 s2
) (numerical-equal s1 s2
:test test
))
656 (defmethod numerical-equal ((result1 list
) (result2 list
)
657 &key
(test #'number-equal
))
658 "Return true if the lists are equal in length and each element is
659 equal according to :TEST."
660 (%sequence-equal result1 result2 test
))
662 (defmethod numerical-equal ((result1 vector
) (result2 vector
)
663 &key
(test #'number-equal
))
664 "Return true if the vectors are equal in length and each element is
665 equal according to :TEST."
666 (%sequence-equal result1 result2 test
))
668 (defmethod numerical-equal ((result1 list
) (result2 vector
)
669 &key
(test #'number-equal
))
670 "Return true if every element of the list is equal to the
671 corresponding element of the vector."
672 (%sequence-equal result1 result2 test
))
674 (defmethod numerical-equal ((result1 vector
) (result2 list
)
675 &key
(test #'number-equal
))
676 "Return true if every element of the list is equla to the
677 corresponding element of the vector."
678 (%sequence-equal result1 result2 test
))
680 (defmethod numerical-equal ((result1 array
) (result2 array
)
681 &key
(test #'number-equal
))
682 "Return true if the arrays are equal in dimension and each element
683 is equal according to :TEST."
684 (when (equal (array-dimensions result1
) (array-dimensions result2
))
686 (make-array (array-total-size result1
)
687 :element-type
(array-element-type result1
)
688 :displaced-to result1
)
689 (make-array (array-total-size result2
)
690 :element-type
(array-element-type result2
)
691 :displaced-to result2
))))
693 (defmacro assert-numerical-equal
(expected form
&rest extras
)
694 `(expand-assert :equal
,form
,form
,expected
,extras
:test
#'numerical-equal
))
696 ;;; FIXME : Audit and move the diagnostic functions to a separate
699 ;;; Diagnostic functions
700 ;;; Failing a unit test is only half the problem.
702 ;;; Sequence errors and the indices.
703 (defun %sequence-error
(sequence1 sequence2 test error-function
)
704 "Return a sequence of the indice and error between the sequences."
705 (let ((n1 nil
) (n2 nil
)
707 (dotimes (index (length sequence1
) errseq
)
708 (setf n1
(elt sequence1 index
)
709 n2
(elt sequence2 index
))
710 (unless (funcall test n1 n2
)
711 (push (list (1- index
) n1 n2
(funcall error-function n1 n2
))
714 (defun sequence-error (sequence1 sequence2
&key
715 (test #'number-equal
)
716 (error-function #'relative-error
))
717 "Return a sequence of the indice and error between the sequence elements."
719 ((not (typep sequence1
'sequence
))
720 (error "SEQUENCE1 is not a valid sequence."))
721 ((not (typep sequence2
'sequence
))
722 (error "SEQUENCE2 is not a valid sequence."))
723 ((not (= (length sequence1
) (length sequence2
)))
724 (error "Lengths not equal. SEQUENCE1(~D) /= SEQUENCE2(~D)."
725 (length sequence1
) (length sequence2
)))
726 (t (%sequence-error sequence1 sequence2 test error-function
))))
728 ;;; Array errors and the indices.
729 (defun %array-indices
(row-major-index dimensions
)
730 "Recursively calculate the indices from the row major index."
731 (let ((remaining (rest dimensions
)))
733 (multiple-value-bind (index remainder
)
734 (floor row-major-index
(reduce #'* remaining
))
735 (cons index
(%array-indices remainder remaining
)))
736 (cons row-major-index nil
))))
738 (defun %array-error
(array1 array2 test errfun
)
739 "Return a list of the indices, values and error of the elements that
741 (let ((dimensions (array-dimensions array1
))
745 (dotimes (index (array-total-size array1
) errseq
)
746 (setf indices
(%array-indices index dimensions
)
747 n1
(apply #'aref array1 indices
)
748 n2
(apply #'aref array2 indices
))
749 (unless (funcall test n1 n2
)
750 (push (list indices n1 n2
(funcall errfun n1 n2
))
753 (defun array-error (array1 array2
&key
754 (test #'number-equal
)
755 (error-function #'relative-error
))
756 "Return a list of the indices and error between the array elements."
758 ((not (arrayp array1
))
759 (error "ARRAY1 is not an array."))
760 ((not (arrayp array2
))
761 (error "ARRAY2 is not an array."))
762 ((not (equal (array-dimensions array1
) (array-dimensions array2
)))
763 (error "Arrays are not equal dimensions."))
764 (t (%array-error array1 array2 test error-function
))))
766 ;;; Floating point data functions
767 (defun make-2d-list (rows columns
&key
(initial-element 0))
768 "Return a nested list with INITIAL-ELEMENT."
769 (mapcar (lambda (x) (make-list columns
:initial-element x
))
770 (make-list rows
:initial-element initial-element
)))
772 (defun %complex-float-random
(limit &optional
(state *random-state
*))
773 "Return a random complex float number."
775 (random (realpart limit
) state
)
776 (random (imagpart limit
) state
)))
778 (defun %complex-rational-random
(limit &optional
(state *random-state
*))
779 "Return a random complex rational number."
780 (let ((imaglimit (imagpart limit
)))
783 (random (realpart limit
) state
)
784 ;; Ensure that the imaginary part is not zero.
785 (do ((im (random imaglimit state
)
786 (random imaglimit state
)))
788 (error "Imaginary part must be greater than 1."))))
790 (defun complex-random (limit &optional
(state *random-state
*))
791 "Return a random complex number. "
792 (check-type limit complex
)
793 (if (typep limit
'(complex rational
))
794 (%complex-rational-random limit state
)
795 (%complex-float-random limit state
)))
797 (defun make-random-list (size &optional
(limit 1.0))
798 "Return a list of random numbers."
799 (mapcar (if (complexp limit
) #'complex-random
#'random
)
800 (make-list size
:initial-element limit
)))
802 (defun make-random-2d-list (rows columns
&optional
(limit 1.0))
803 "Return a nested list of random numbers."
804 (mapcar (lambda (x) (make-random-list columns x
))
805 (make-list rows
:initial-element limit
)))
807 (defun make-random-2d-array (rows columns
&optional
(limit 1.0))
808 "Return a 2D array of random numbers."
809 (let ((new-array (make-array (list rows columns
)
810 :element-type
(type-of limit
)))
811 (random-func (if (complexp limit
)
814 (dotimes (i0 rows new-array
)
815 (dotimes (i1 columns
)
816 (setf (aref new-array i0 i1
)
817 (funcall random-func limit
))))))