Add methods to magicl:=-lisp so that magicl:= will compare scalars. by svspire · Pull Request #154 · quil-lang/magicl

svspire · 2021-12-12T17:21:33Z

Add #'test-scalar-equality to test those methods.
Modify #'test-p-norm to use magicl:= rather than common-lisp:= because
it's (potentially) comparing floats.

Add #'test-scalar-equality to test those methods. Modify #'test-p-norm to use magicl:= rather than common-lisp:= because it's (potentially) comparing floats.

svspire · 2021-12-12T17:34:42Z

Some of things to note:

Search for "Use least precise epsilon" in single-float.lisp. Is that correct, or should we be using the most precise epsilon because of FP contagion?
Is #'test-scalar-equality thorough enough?
Should we care about short-floats and long-floats too? In my implementation (CCL) they're synonymous with single-float and double-float, respectively.

ghost · 2022-01-07T22:31:36Z

Re

Should we care about short-floats and long-floats too? In my implementation (CCL) they're synonymous with single-float and double-float, respectively.

we have the same sameness in our primary implementation (SBCL) and in emerging Allegro implementation. Ignoring types short- and long-float is preexisting practice, so I think it's OK to continue.

ghost · 2022-01-08T02:18:10Z

Re

Search for "Use least precise epsilon" in single-float.lisp. Is that correct, or should we be using the most precise epsilon because of FP contagion?

I believe the original code is not correct, and we should be using most precise (double-float) epsilon by default because of floating-point contagion, namely *double-comparison-threshold*.

Consider that (= 1.0s0 1.0d0) should be like (= 1.0d0 1.0d0) by default. And (= 1.0d0 (+ 1.0d0 *double-comparison-threshold*)) is true, whereas (= 1.0d0 (+ 1.0d0 *double-comparison-threshold* *double-comparison-threshold*)) is false. Therefore, (= 1.0s0 (+ 1.0d0 *double-comparison-threshold* *double-comparison-threshold*)) should be false by default.

In the original implementation, specifying epsilon

MAGICL> (= 1.0s0 (+ 1.0d0 *double-comparison-threshold* *double-comparison-threshold*) *double-comparison-threshold*)
NIL
MAGICL> (= 1.0s0 (+ 1.0d0 *double-comparison-threshold* *double-comparison-threshold*) *float-comparison-threshold*)
T

ghost · 2022-01-08T02:22:12Z

Probably worth a review by @stylewarning to be sure. Here's a reference on CL precision, contagion, and coercion: https://www.cs.cmu.edu/Groups/AI/html/cltl/clm/node122.html

svspire · 2022-01-08T04:42:49Z

I believe the original code is not correct, and we should be using most precise (double-float) epsilon by default because of floating-point contagion, namely *double-comparison-threshold*.

That's not the original code; it's mine. I asked the question because I wasn't sure I'd done the right thing. So if it's wrong it's on me :-)

ghost · 2022-01-08T21:29:10Z

That's not the original code; it's mine.

OK, no problem @svspire -- I just meant for someone reading this, "original" = code being reviewed, before (my) recommended changes. Thanks

stylewarning

(Belated) happy holidays! Sorry for the delay on the review.

I think there's room to simplify (which, if I'm right, should be great news?!). This is definitely a welcome addition to MAGICL.

stylewarning · 2022-01-09T08:38:51Z


 (in-package #:magicl-tests)

+(defmacro swapping-arguments-is ((predicate arg1 arg2))


This can be a function if you'll tolerate a FUNCALL. is doesn't need to be lexically visible.

stylewarning · 2022-01-09T08:39:22Z

+       (is (,predicate ,arg1sym ,arg2sym))
+       (is (,predicate ,arg2sym ,arg1sym)))))
+
+(defmacro swapping-arguments-not ((predicate arg1 arg2))


if swapping-arguments-is is made a function, you can use complement on predicate

stylewarning · 2022-01-09T08:44:50Z

+;; Might want to inline this
+(defun %scalar= (s1 s2 epsilon)
+  "For equality checks of inexact scalars."
+  (declare (type number s1 s2) ; you can use this on complex but it might be slower because of the sqrt and two multiplies. Maybe.
+           (type real epsilon))
+  (<= (abs (- s1
+              s2))
+      epsilon))
+
+(defmethod =-lisp ((scalar1 rational) (scalar2 rational) &optional epsilon)
+  (declare (ignore epsilon))
+   "Rationals (integers and ratios) should be compared exactly."
+  (common-lisp:= scalar1 scalar2))
+


Could we move this to a separate file like numerical-equality.lisp? The extensible function defined in abstract-tensor should probably be moved there too. That way we can knock out a bunch of general/generic stuff all in one go, and we don't have it hidden here in single-float.lisp

stylewarning · 2022-01-09T08:46:42Z

+(defmethod =-lisp ((scalar1 single-float) (scalar2 single-float) &optional (epsilon *float-comparison-threshold*))
+  (%scalar= scalar1 scalar2 epsilon))
+
+(defmethod =-lisp ((scalar1 rational) (scalar2 single-float) &optional (epsilon *float-comparison-threshold*))
+  (%scalar= scalar1 scalar2 epsilon))
+
+(defmethod =-lisp ((scalar1 single-float) (scalar2 rational) &optional (epsilon *float-comparison-threshold*))
+  (%scalar= scalar1 scalar2 epsilon))
+
+(defmethod =-lisp ((scalar1 float) (scalar2 single-float) &optional (epsilon *float-comparison-threshold*))
+  "This covers comparing double-float to single-float. Use least precise epsilon."
+  (%scalar= scalar1 scalar2 epsilon))
+
+(defmethod =-lisp ((scalar1 single-float) (scalar2 float) &optional (epsilon *float-comparison-threshold*))
+  "This covers comparing single-float to double-float. Use least precise epsilon."
+  (%scalar= scalar1 scalar2 epsilon))
+


If we make the separate file, could we "just" implement a single defmethod on real and real to cover all these cases not covered by rational?

stylewarning · 2022-01-09T08:47:10Z

+(defmethod =-lisp ((scalar1 double-float) (scalar2 double-float) &optional (epsilon *double-comparison-threshold*))
+  (%scalar= scalar1 scalar2 epsilon))
+
+(defmethod =-lisp ((scalar1 rational) (scalar2 double-float) &optional (epsilon *double-comparison-threshold*))
+ (%scalar= scalar1 scalar2 epsilon))
+
+(defmethod =-lisp ((scalar1 double-float) (scalar2 rational) &optional (epsilon *double-comparison-threshold*))
+  (%scalar= scalar1 scalar2 epsilon))
+


See comments elsewhere. I think these might be able to be eliminated.

stylewarning · 2022-01-09T08:48:10Z

  (loop :for i :below (size vector1)
        :sum (* (tref vector1 i) (conjugate (tref vector2 i)))))

+(defmethod =-lisp ((scalar1 complex) (scalar2 complex) &optional epsilon)


See comment elsewhere about moving these functions. I think this should be moved to the file I mention elsewhere.

stylewarning · 2022-01-09T08:52:48Z

+(defmethod =-lisp ((scalar1 complex) (scalar2 real) &optional epsilon)
+  (let ((imagpart (imagpart scalar1))
+        (zero nil))
+    (typecase imagpart ;; TODO: do we need to care about short-floats and long-floats? On which implementations does it matter?
+      (single-float (setf epsilon (or epsilon *float-comparison-threshold*)
+                          zero 0.0s0))
+      (double-float (setf epsilon (or epsilon *double-comparison-threshold*)
+                          zero 0.0d0))
+      ; If imagpart is rational it's also guaranteed nonzero per ANSI, which means we can just go ahead and exit now
+      (t (return-from =-lisp nil)))
+    (and (%scalar= zero imagpart epsilon)
+         (%scalar= (realpart scalar1) scalar2 epsilon))))
+
+(defmethod =-lisp ((scalar1 real) (scalar2 complex) &optional epsilon)
+  (let ((imagpart (imagpart scalar2))
+        (zero nil))
+    (typecase imagpart
+      (single-float (setf epsilon (or epsilon *float-comparison-threshold*)
+                          zero 0.0s0))
+      (double-float (setf epsilon (or epsilon *double-comparison-threshold*)
+                          zero 0.0d0))
+      ; If imagpart is rational it's also guaranteed nonzero per ANSI, which means we can just go ahead and exit now
+      (t (return-from =-lisp nil)))
+    (and (%scalar= zero imagpart epsilon)
+         (%scalar= (realpart scalar2) scalar1 epsilon))))


I might be missing some nuance but I think these can be implemented as

(and (%scalar= 0 (realpart z) epsilon) (%scalar= x (imagpart z) epsilon))

I can't imagine finding a default for epsilon and using a typecase is needed.

stylewarning · 2022-01-09T08:53:11Z

+       (is (not (,predicate ,arg1sym ,arg2sym)))
+       (is (not (,predicate ,arg2sym ,arg1sym))))))
+
+(deftest test-scalar-equality ()


great and thorough tests! Thanks for writing them.

ghost · 2022-01-27T22:54:59Z

@stylewarning re:

I can't imagine finding a default for epsilon and using a typecase is needed.

I don't really understand your comment. Recalling that there already was defaulting of epsilon in preexisting code, are you objecting to that here?

stylewarning · 2022-01-31T20:28:46Z

@stylewarning re:

I can't imagine finding a default for epsilon and using a typecase is needed.

I don't really understand your comment. Recalling that there already was defaulting of epsilon in preexisting code, are you objecting to that here?

I don't understand why we aren't just recursing on =-lisp instead of baking in the typecase dispatching ourselves.

ghost · 2022-02-04T21:53:06Z

@stylewarning re:

I can't imagine finding a default for epsilon and using a typecase is needed.

I don't really understand your comment. Recalling that there already was defaulting of epsilon in preexisting code, are you objecting to that here?

I don't understand why we aren't just recursing on =-lisp instead of baking in the typecase dispatching ourselves.

Aha, I see the code is figuring out what to pass if a caller does not supply epsilon (or passes in nil). It's a bit of a pain in Lisp to pass along optionally supplied value if and only if it's supplied, and otherwise default, but it can be done using a supplied-p parameter with the optional arg (as described here: http://clhs.lisp.se/Body/03_dab.htm).

But as things stand, code still ends up calling %scalar=, which takes a required epsilon arg. You could change %scalar= to have epsilon be optional and to take care of figuring out the default for epsilon (which I think should be the highest resolution threshold, i.e., *double-comparison-threshold*), so then, at least, you'd have it all in just one place.

But taking @stylewarning's suggestion to heart, you could recursively call =-lisp in each case instead. The issue is we did not have a version of that for mixed types. So I think this would mean systematically eliminating calls to %scalar= and replacing them with recursive calls to =-lisp with args coerced to the "appropriate" type. For a single- and double-float mixed case, coerce the single to double. For a rational and <float-type> mixed case, coerce the rational to <float-type>. This is just a sketch. What do folks think?

Add methods to magicl:=-lisp so that magicl:= will compare scalars.

4355013

Add #'test-scalar-equality to test those methods. Modify #'test-p-norm to use magicl:= rather than common-lisp:= because it's (potentially) comparing floats.

stylewarning requested changes Jan 9, 2022

View reviewed changes


		(in-package #:magicl-tests)

		(defmacro swapping-arguments-is ((predicate arg1 arg2))

Conversation

svspire commented Dec 12, 2021

Uh oh!

svspire commented Dec 12, 2021

Uh oh!

ghost commented Jan 7, 2022

Uh oh!

ghost commented Jan 8, 2022 • edited by ghost Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ghost commented Jan 8, 2022

Uh oh!

svspire commented Jan 8, 2022

Uh oh!

ghost commented Jan 8, 2022

Uh oh!

stylewarning left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ghost commented Jan 27, 2022

Uh oh!

stylewarning commented Jan 31, 2022

Uh oh!

ghost commented Feb 4, 2022 • edited by ghost Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

ghost commented Jan 8, 2022 •

edited by ghost

Loading

ghost commented Feb 4, 2022 •

edited by ghost

Loading