ch6 tests behave differently depending on seemingly equivalent series g/-

[sritchie]

Here's a reproduction. You can make this test fail by redefining the existing seq:- in series.cljc:

(defn- seq:- [f g]
  (seq:+ f (seq:negate g)))

to:

(defn- seq:- [f g]
  (map g/- f g))

Which should be equivalent. Here's the check:

(deftest ch6-repro
  "Evaluate this in sicmutils.sicm.ch6_test.cljc. If you defined g/- for series
  as `(map g/- f g)`, this test will return a different result than if you
  define it as `(g/+ f (g/negate g))`."
  (let [L (e/Lie-derivative
           (fn [[_ theta ptheta]]
             (/ (* -1 'a 'b (sin theta))
                ptheta)))
        H (e/series
           (fn [[_ _ ptheta]]
             (/ (square ptheta)
                (* 2 'a)))
           (fn [[_ theta _]]
             (* -1 'e 'b (cos theta))))
        E (((exp (* 'e L)) H) (up 't 'theta 'p_theta))]
    (is (= '((/ (expt p_theta 2) (* 2 a))
             0
             (/ (* a (expt b 2) (expt e 2) (expt (sin theta) 2))
                (* 2 (expt p_theta 2))))
           (take 3 (simplify E))))))

With the existing definition, the 3rd sequence entry is

(/ (* a (expt b 2) (expt e 2) (expt (sin theta) 2)) (* 2 (expt p_theta 2)))

If you use the new g/- based definition, you pick up a factor of 3 in the third term:

(/ (* 3N a (expt b 2) (expt e 2) (expt (sin theta) 2)) (* 2 (expt p_theta 2)))

Seems like a bug with the simplifier.

[sritchie]

@littleredcomputer , does this ring any bells?

[littleredcomputer]

It does not. This needs investigation, which I can do starting 11/9, given other commitments. The details of the simplifier is another field I want to look at. Doing some quickie, `printf` style instrumentation, I found that _a lot_ of time is spent, for example, canonicalizing the order of mixed partial derivatives, essentially, bubble-sorting them by slot number. You can't really partially differentiate something that many times in the Universe in which we're living, so the O(n^2) isn't really a problem, but rather the fact that we keep checking for it. Made me think a good optimization would be just to handle that at construction time so that there was no way to create a form in which the slot variables were unsorted, so I want to look at that too. I will check the patterns for the trig substitutions, since it's clear there may be a sign error in there

…

On Fri, Nov 6, 2020 at 7:38 AM Sam Ritchie ***@***.***> wrote: @littleredcomputer <https://github.com/littleredcomputer> , does this ring any bells? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#151 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAB3Y7PQ3OX7VI2JQ3DZ5T3SOQKALANCNFSM4TM2I4YQ> .

[sritchie]

Yes, @littleredcomputer another thing I noticed, and wasn't sure where to investigate, is that the simplifier sorts * and + calls in different order. Obviously we need to be careful here for non-commutative operations, but I know that there are many tests where we get different results (just a different order) if we don't use the hermetic-simplify-fixture.

[sritchie]

@littleredcomputer now that I've stared at some rules code, I see the same comment from Alexey (bubble sort!) in his "commutativity" code. He handles this by matching the entire list of arguments and sorting the whole thing at once.
https://github.com/axch/rules/blob/master/simplifiers.scm#L52

Can we do something similar, and match a full stack of partials? I think the original scmutils code has rules to:

• turn nests into products, handling partials that are exponentiated
• THEN sort the terms inside the product

But the latter is bubble-sorted, so we can still make it better.

Here's the scheme:

(def canonicalize-partials*
  (ruleset
   ;; First turn nests into products.
   ((partial :i*) ((partial :j*) :f))
   =>
   ((* (partial :i*) (partial :j*)) :f)

   ((partial :i*) ((* (partial :j*) :more*) :f))
   =>
   ((* (partial :i*) (partial :j*) :more*) :f)

   ;; Exponentiation of operators makes things hairy
   ((expt (partial :i*) :n) ((partial :j*) :f))
   =>
   ((* (expt (partial :i*) :n) (partial :j*)) :f)

   ((partial :i*) ((expt (partial :j*) :n) :f))
   =>
   ((* (partial :i*) (expt (partial :j*) :n)) :f)

   ((expt (partial :i*) :n) ((expt (partial :j*) :m) :f))
   =>
   ((* (expt (partial :i*) :n) (expt (partial :j*) :m)) :f)


   ;; Already some accumulation
   ((partial :i*) ((* (partial :j*) :more*) :f))
   =>
   ((* (partial :i*) (partial :j*) :more*) :f)

   ((expt (partial :i*) :n) ((* (partial :j*) :more*) :f))
   =>
   ((* (expt (partial :i*) :n) (partial :j*) :more*) :f)

   ((partial :i*) ((* (expt (partial :j*) :m) :more*) :f))
   =>
   ((* (partial :i*) (expt (partial :j*) :m) :more*) :f)

   ((expt (partial :i*) :n) ((* (expt (partial :j*) :m) :more*) :f))
   =>
   ((* (expt (partial :i*) :n) (expt (partial :j*) :m) :more*) :f)

   ;; Next sort products, if OK
   #_#_#_(((* :xs* (partial :i*) :ys* (partial :j*) :zs*)
           ;; TODO check that this is waht we mean with symbol
           (:? f symbol?))
          :args*)
   (let ((args (expression args)))
     (and commute-partials?
          (symb:elementary-access? i args)
          (symb:elementary-access? j args)
          (list< j i)))
   (((* :xs* (partial :j*) :ys* (partial :i*) :zs*) :f)
    :args*)
   ))

[sritchie]

This seems to have... just fixed itself, as a result of some other twitches we've had in the codebase! #285 changes the series implementation and marks the close of this mystery bug.