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
375 (defun %seq-1-norm (data)
376 "Return the Taxicab norm of the sequence."
377 ;; FIXME : Use the LOOP.
378 (reduce (lambda (x y) (+ x (abs y)))
379 data :initial-value 0))
381 (defun %seq-2-norm (data)
382 "Return the Euclidean norm of the sequence."
383 (multiple-value-bind (scale sumsq)
384 (sumsq (map-into (make-array (length data)) #'abs data))
385 (* scale (sqrt sumsq))))
387 (defun %seq-p-norm (data p)
388 "Return the p norm of the sequence."
389 (multiple-value-bind (scale sump)
390 (sump (map-into (make-array (length data)) #'abs data) p)
391 (* scale (expt sump (/ p)))))
393 (defun %seq-inf-norm (data)
394 "Return the infinity, or maximum, norm of the sequence."
395 ;; FIXME : Use the LOOP.
396 (reduce (lambda (x y) (max x (abs y)))
397 data :initial-value 0))
399 (defun %seq-norm (data measure)
400 "Return the norm of the sequence according to the measure."
407 (%seq-p-norm data measure))
408 ((equalp measure :infinity)
409 (%seq-inf-norm data))
410 (t (error "Unrecognized norm, ~A." measure))))
412 (defmethod norm ((data list) &optional (measure *measure*))
413 "Return the norm of the list according to the measure."
414 (%seq-norm data measure))
416 (defmethod norm ((data vector) &optional (measure *measure*))
417 "Return the norm of the vector according to the measure."
418 (%seq-norm data measure))
420 (defmethod norm ((data array) &optional (measure *measure*))
421 "Return the entrywise norm of the array according to the measure."
422 (let ((flat-data (make-array (array-total-size data)
423 :element-type (array-element-type data)
424 :displaced-to data)))
426 ((and (numberp measure) (< 0 measure))
427 (warn "Measure ~D results in an entrywise p-norm." measure)
428 (%seq-p-norm flat-data measure))
429 ((equalp measure :frobenius)
430 (%seq-2-norm flat-data))
431 ((equalp measure :max)
432 (%seq-inf-norm flat-data))
433 (t (error "Unrecognized norm, ~A." measure)))))
435 ;;; (RELATIVE-ERROR-NORM exact approximate measure) => float
436 (defun %relative-error-norm (exact approximate measure)
437 "Return the relative error norm of the sequences."
438 (/ (norm (map-into (make-array (length exact))
439 (lambda (x1 x2) (abs (- x1 x2)))
440 exact approximate) measure)
441 (norm exact measure)))
443 (defmethod relative-error-norm ((exact list) (approximate list)
444 &optional (measure *measure*))
445 "Return the relative error norm of the lists."
446 (if (= (length exact) (length approximate))
447 (%relative-error-norm exact approximate measure)
448 (error "Lists are not equal in length.")))
450 (defmethod relative-error-norm ((exact list) (approximate vector)
451 &optional (measure *measure*))
452 "Return the relative error norm of the list and the vector."
453 (if (= (length exact) (length approximate))
454 (%relative-error-norm exact approximate measure)
455 (error "The list and vector are not equal in length.")))
457 (defmethod relative-error-norm ((exact vector) (approximate list)
458 &optional (measure *measure*))
459 "Return the relative error norm of the list and the vector."
460 (if (= (length exact) (length approximate))
461 (%relative-error-norm exact approximate measure)
462 (error "The list and vector are not equal in length.")))
464 (defmethod relative-error-norm ((exact vector) (approximate vector)
465 &optional (measure *measure*))
466 "Return the relative error norm of the vectors."
467 (if (= (length exact) (length approximate))
468 (%relative-error-norm exact approximate measure)
469 (error "Vectors are not equal in length.")))
471 (defmethod relative-error-norm ((exact array) (approximate vector)
472 &optional (measure *measure*))
473 "Return the relative error norm of the arrays."
474 (if (equal (array-dimensions exact)
475 (array-dimensions approximate))
476 (%relative-error-norm
477 (make-array (array-total-size exact)
478 :element-type (array-element-type exact)
480 (make-array (array-total-size approximate)
481 :element-type (array-element-type approximate)
482 :displaced-to approximate)
484 (error "Arrays are not equal dimensions.")))
486 ;;; (NORM-EQUAL data1 data2 epsilon measure) => boolean
487 (defun %norm-equal (seq1 seq2 epsilon measure)
488 "Return true if the relative error norm is less than epsilon."
490 (and (null seq1) (null seq2))
491 (< (%relative-error-norm seq1 seq2 measure) epsilon)))
493 (defmethod norm-equal ((data1 list) (data2 list) &optional
494 (epsilon *epsilon*) (measure *measure*))
495 "Return true if the lists are equal in length and the relative error
496 norm is less than epsilon."
497 (%norm-equal data1 data2 epsilon measure))
499 (defmethod norm-equal ((data1 list) (data2 vector) &optional
500 (epsilon *epsilon*) (measure *measure*))
501 "Return true if the vector and the list are equal in length and the
502 relative error norm is less than epsilon."
503 (%norm-equal data1 data2 epsilon measure))
505 (defmethod norm-equal ((data1 vector) (data2 list) &optional
506 (epsilon *epsilon*) (measure *measure*))
507 "Return true if the vector and the list are equal in length and the
508 relative error norm is less than epsilon."
509 (%norm-equal data1 data2 epsilon measure))
511 (defmethod norm-equal ((data1 vector) (data2 vector) &optional
512 (epsilon *epsilon*) (measure *measure*))
513 "Return true if the vectors are equal in length and the relative
514 error norm is less than epsilon."
515 (%norm-equal data1 data2 epsilon measure))
517 (defmethod norm-equal ((data1 array) (data2 array) &optional
518 (epsilon *epsilon*) (measure *measure*))
519 "Return true if the arrays are equal in length and the relative
520 error norm is less than epsilon."
521 (when (equal (array-dimensions data1)
522 (array-dimensions data2))
524 (make-array (array-total-size data1)
525 :element-type (array-element-type data1)
527 (make-array (array-total-size data2)
528 :element-type (array-element-type data2)
532 (defmacro assert-norm-equal (expected form &rest extras)
533 `(expand-assert :equal ,form ,form ,expected ,extras :test #'norm-equal))
535 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
536 ;;; [NumAlgoC] : Definition
1.7, pg.
4
538 ;;; To avoid using 0.1, first 1.0 <= significand < 10. On the final
539 ;;; return, scale 0.1 <= significand < 1.
540 (defun %normalize-float
(significand &optional
(exponent 0))
541 "Return the normalized floating point number and exponent."
542 ;;; FIXME : Use the LOOP.
545 (values significand
0))
546 ((>= (abs significand
) 10)
547 (%normalize-float
(/ significand
10.0) (1+ exponent
)))
548 ((< (abs significand
) 1)
549 (%normalize-float
(* significand
10.0) (1- exponent
)))
550 (t (values (/ significand
10.0) (1+ exponent
)))))
552 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
553 (defun %sigfig-equal
(float1 float2 significant-figures
)
554 "Return true if the floating point numbers have equal significant
556 (if (or (zerop float1
) (zerop float2
))
557 (< (abs (+ float1 float2
)) (* 5D-1
(expt 1D1
(- significant-figures
))))
558 (multiple-value-bind (sig1 exp1
) (%normalize-float float1
)
559 (multiple-value-bind (sig2 exp2
) (%normalize-float float2
)
560 (= (round (* sig1
(expt 1D1 significant-figures
)))
561 (round (* sig2
(expt 1D1
(- significant-figures
(- exp1 exp2
))))))))))
563 (defmethod sigfig-equal ((data1 float
) (data2 float
) &optional
564 (significant-figures *significant-figures
*))
565 "Return true if the floating point numbers have equal significant
567 (%sigfig-equal data1 data2 significant-figures
))
569 (defun %seq-sigfig-equal
(seq1 seq2 significant-figures
)
570 "Return true if the element-wise comparison is equal to the
571 specified significant figures."
573 (and (null seq1
) (null seq2
))
574 (when (= (length seq1
) (length seq2
))
576 (lambda (d1 d2
) (sigfig-equal d1 d2 significant-figures
))
579 (defmethod sigfig-equal ((data1 list
) (data2 list
) &optional
580 (significant-figures *significant-figures
*))
581 "Return true if the lists are equal in length and the element-wise
582 comparison is equal to significant figures."
583 (%seq-sigfig-equal data1 data2 significant-figures
))
585 (defmethod sigfig-equal ((data1 vector
) (data2 list
) &optional
586 (significant-figures *significant-figures
*))
587 "Return true if the vector and the list are equal in length and the
588 element-wise comparison is equal to significant figures."
589 (%seq-sigfig-equal data1 data2 significant-figures
))
591 (defmethod sigfig-equal ((data1 list
) (data2 vector
) &optional
592 (significant-figures *significant-figures
*))
593 "Return true if the list and the vector are equal in length and the
594 element-wise comparison is equal to significant figures."
595 (%seq-sigfig-equal data1 data2 significant-figures
))
597 (defmethod sigfig-equal ((data1 vector
) (data2 vector
) &optional
598 (significant-figures *significant-figures
*))
599 "Return true if the vectors are equal in length and the element-wise
600 comparison is equal to significant figures."
601 (%seq-sigfig-equal data1 data2 significant-figures
))
603 (defmethod sigfig-equal ((data1 array
) (data2 array
) &optional
604 (significant-figures *significant-figures
*))
605 "Return true if the arrays are equal in length and the element-wise
606 comparison is equal to significant figures."
607 (when (equal (array-dimensions data1
)
608 (array-dimensions data2
))
610 (make-array (array-total-size data1
)
611 :element-type
(array-element-type data1
)
613 (make-array (array-total-size data2
)
614 :element-type
(array-element-type data2
)
616 significant-figures
)))
618 (defmacro assert-sigfig-equal
(expected form
&rest extras
)
619 `(expand-assert :equal
,form
,form
,expected
,extras
:test
#'sigfig-equal
))
621 ;;; (NUMBER-EQUAL number1 number2) => true or false
622 (defun number-equal (number1 number2
&optional
(epsilon *epsilon
*) type-eq-p
)
623 "Return true if the numbers are equal within some epsilon,
624 optionally requiring the types to be identical."
626 (or (not type-eq-p
) (eq (type-of number1
) (type-of number2
)))
627 (float-equal (coerce number1
'(complex double-float
))
628 (coerce number2
'(complex double-float
))
631 (defmacro assert-number-equal
(expected form
&rest extras
)
632 `(expand-assert :equal
,form
,form
,expected
,extras
:test
#'number-equal
))
634 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
636 ;;; This is a universal wrapper created by Liam Healy. It is
637 ;;; implemented to support testing in GSLL. The interface is expanded,
638 ;;; but backwards compatible with previous versions.
640 (defmethod numerical-equal ((result1 number
) (result2 number
)
641 &key
(test #'number-equal
))
642 "Return true if the the numbers are equal according to :TEST."
643 (funcall test result1 result2
))
645 (defun %sequence-equal
(seq1 seq2 test
)
646 "Return true if the sequences are equal in length and each element
647 is equal according to :TEST."
648 (when (= (length seq1
) (length seq2
))
649 (every (lambda (s1 s2
) (numerical-equal s1 s2
:test test
))
652 (defmethod numerical-equal ((result1 list
) (result2 list
)
653 &key
(test #'number-equal
))
654 "Return true if the lists are equal in length and each element is
655 equal according to :TEST."
656 (%sequence-equal result1 result2 test
))
658 (defmethod numerical-equal ((result1 vector
) (result2 vector
)
659 &key
(test #'number-equal
))
660 "Return true if the vectors are equal in length and each element is
661 equal according to :TEST."
662 (%sequence-equal result1 result2 test
))
664 (defmethod numerical-equal ((result1 list
) (result2 vector
)
665 &key
(test #'number-equal
))
666 "Return true if every element of the list is equal to the
667 corresponding element of the vector."
668 (%sequence-equal result1 result2 test
))
670 (defmethod numerical-equal ((result1 vector
) (result2 list
)
671 &key
(test #'number-equal
))
672 "Return true if every element of the list is equla to the
673 corresponding element of the vector."
674 (%sequence-equal result1 result2 test
))
676 (defmethod numerical-equal ((result1 array
) (result2 array
)
677 &key
(test #'number-equal
))
678 "Return true if the arrays are equal in dimension and each element
679 is equal according to :TEST."
680 (when (equal (array-dimensions result1
) (array-dimensions result2
))
682 (make-array (array-total-size result1
)
683 :element-type
(array-element-type result1
)
684 :displaced-to result1
)
685 (make-array (array-total-size result2
)
686 :element-type
(array-element-type result2
)
687 :displaced-to result2
))))
689 (defmacro assert-numerical-equal
(expected form
&rest extras
)
690 `(expand-assert :equal
,form
,form
,expected
,extras
:test
#'numerical-equal
))
692 ;;; FIXME : Audit and move the diagnostic functions to a separate
695 ;;; Diagnostic functions
696 ;;; Failing a unit test is only half the problem.
698 ;;; Sequence errors and the indices.
699 (defun %sequence-error
(sequence1 sequence2 test error-function
)
700 "Return a sequence of the indice and error between the sequences."
701 (let ((n1 nil
) (n2 nil
)
703 (dotimes (index (length sequence1
) errseq
)
704 (setf n1
(elt sequence1 index
)
705 n2
(elt sequence2 index
))
706 (unless (funcall test n1 n2
)
707 (push (list (1- index
) n1 n2
(funcall error-function n1 n2
))
710 (defun sequence-error (sequence1 sequence2
&key
711 (test #'number-equal
)
712 (error-function #'relative-error
))
713 "Return a sequence of the indice and error between the sequence elements."
715 ((not (typep sequence1
'sequence
))
716 (error "SEQUENCE1 is not a valid sequence."))
717 ((not (typep sequence2
'sequence
))
718 (error "SEQUENCE2 is not a valid sequence."))
719 ((not (= (length sequence1
) (length sequence2
)))
720 (error "Lengths not equal. SEQUENCE1(~D) /= SEQUENCE2(~D)."
721 (length sequence1
) (length sequence2
)))
722 (t (%sequence-error sequence1 sequence2 test error-function
))))
724 ;;; Array errors and the indices.
725 (defun %array-indices
(row-major-index dimensions
)
726 "Recursively calculate the indices from the row major index."
727 (let ((remaining (rest dimensions
)))
729 (multiple-value-bind (index remainder
)
730 (floor row-major-index
(reduce #'* remaining
))
731 (cons index
(%array-indices remainder remaining
)))
732 (cons row-major-index nil
))))
734 (defun %array-error
(array1 array2 test errfun
)
735 "Return a list of the indices, values and error of the elements that
737 (let ((dimensions (array-dimensions array1
))
741 (dotimes (index (array-total-size array1
) errseq
)
742 (setf indices
(%array-indices index dimensions
)
743 n1
(apply #'aref array1 indices
)
744 n2
(apply #'aref array2 indices
))
745 (unless (funcall test n1 n2
)
746 (push (list indices n1 n2
(funcall errfun n1 n2
))
749 (defun array-error (array1 array2
&key
750 (test #'number-equal
)
751 (error-function #'relative-error
))
752 "Return a list of the indices and error between the array elements."
754 ((not (arrayp array1
))
755 (error "ARRAY1 is not an array."))
756 ((not (arrayp array2
))
757 (error "ARRAY2 is not an array."))
758 ((not (equal (array-dimensions array1
) (array-dimensions array2
)))
759 (error "Arrays are not equal dimensions."))
760 (t (%array-error array1 array2 test error-function
))))
762 ;;; Floating point data functions
763 (defun make-2d-list (rows columns
&key
(initial-element 0))
764 "Return a nested list with INITIAL-ELEMENT."
765 (mapcar (lambda (x) (make-list columns
:initial-element x
))
766 (make-list rows
:initial-element initial-element
)))
768 (defun %complex-float-random
(limit &optional
(state *random-state
*))
769 "Return a random complex float number."
771 (random (realpart limit
) state
)
772 (random (imagpart limit
) state
)))
774 (defun %complex-rational-random
(limit &optional
(state *random-state
*))
775 "Return a random complex rational number."
776 (let ((imaglimit (imagpart limit
)))
779 (random (realpart limit
) state
)
780 ;; Ensure that the imaginary part is not zero.
781 (do ((im (random imaglimit state
)
782 (random imaglimit state
)))
784 (error "Imaginary part must be greater than 1."))))
786 (defun complex-random (limit &optional
(state *random-state
*))
787 "Return a random complex number. "
788 (check-type limit complex
)
789 (if (typep limit
'(complex rational
))
790 (%complex-rational-random limit state
)
791 (%complex-float-random limit state
)))
793 (defun make-random-list (size &optional
(limit 1.0))
794 "Return a list of random numbers."
795 (mapcar (if (complexp limit
) #'complex-random
#'random
)
796 (make-list size
:initial-element limit
)))
798 (defun make-random-2d-list (rows columns
&optional
(limit 1.0))
799 "Return a nested list of random numbers."
800 (mapcar (lambda (x) (make-random-list columns x
))
801 (make-list rows
:initial-element limit
)))
803 (defun make-random-2d-array (rows columns
&optional
(limit 1.0))
804 "Return a 2D array of random numbers."
805 (let ((new-array (make-array (list rows columns
)
806 :element-type
(type-of limit
)))
807 (random-func (if (complexp limit
)
810 (dotimes (i0 rows new-array
)
811 (dotimes (i1 columns
)
812 (setf (aref new-array i0 i1
)
813 (funcall random-func limit
))))))