From cac1791c03115b042cca9c709d913f677e4175d2 Mon Sep 17 00:00:00 2001 From: Sam Ritchie Date: Tue, 13 Aug 2024 11:18:18 -0400 Subject: [PATCH] Change default derivative impl to reverse mode (#170) - #170: - changes `D` and `partial` to use reverse mode automatic differentiation by default, and fixes all associated tests - adds `emmy.generic/{zero?,one?,identity?}` implementations (all false) to `emmy.tape/Completed`, in case some collection type tries to simplify these during reverse-mode AD --- CHANGELOG.md | 9 ++++ src/emmy/calculus/derivative.cljc | 16 +++--- src/emmy/calculus/manifold.cljc | 1 - src/emmy/calculus/vector_calculus.cljc | 5 +- src/emmy/dual.cljc | 7 +++ src/emmy/mechanics/hamilton.cljc | 18 +++---- src/emmy/numbers.cljc | 4 +- test/emmy/abstract/number_test.cljc | 11 ++-- test/emmy/calculus/covariant_test.cljc | 72 ++++++++++++------------- test/emmy/calculus/form_field_test.cljc | 40 +++++++------- test/emmy/calculus/map_test.cljc | 20 +++---- test/emmy/fdg/ch4_test.cljc | 29 +++++----- test/emmy/fdg/ch7_test.cljc | 4 +- test/emmy/fdg/ch9_test.cljc | 8 +-- test/emmy/fdg/einstein_test.cljc | 2 +- test/emmy/mechanics/hamilton_test.cljc | 60 ++++++++++----------- test/emmy/mechanics/rigid_test.cljc | 25 +++++---- test/emmy/series_test.cljc | 8 +-- test/emmy/sicm/ch2_test.cljc | 10 ++-- 19 files changed, 182 insertions(+), 167 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ea065097..147d61cb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,15 @@ ## [unreleased] +- #170: + + - changes `D` and `partial` to use reverse mode automatic differentiation by + default, and fixes all associated tests + + - adds `emmy.generic/{zero?,one?,identity?}` implementations (all false) to + `emmy.tape/Completed`, in case some collection type tries to simplify these + during reverse-mode AD + - #185: - adds a dynamic variable `emmy.calculus.derivative/*mode*` that allows the diff --git a/src/emmy/calculus/derivative.cljc b/src/emmy/calculus/derivative.cljc index 57868007..23036e13 100644 --- a/src/emmy/calculus/derivative.cljc +++ b/src/emmy/calculus/derivative.cljc @@ -237,7 +237,7 @@ ;; implementation for the components. I vote to back out this `::s/structure` ;; installation. -(def ^:dynamic *mode* d/FORWARD-MODE) +(def ^:dynamic *mode* d/REVERSE-MODE) (doseq [t [::v/function ::s/structure]] (defmethod g/partial-derivative [t v/seqtype] [f selectors] @@ -266,9 +266,9 @@ the [Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) of `f`." (o/make-operator - (fn [x] + (fn [f] (binding [*mode* d/FORWARD-MODE] - (g/partial-derivative x []))) + (g/partial-derivative f []))) g/derivative-symbol)) (def ^{:arglists '([f])} @@ -282,9 +282,9 @@ the [Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) of `f`." (o/make-operator - (fn [x] + (fn [f] (binding [*mode* d/REVERSE-MODE] - (g/partial-derivative x []))) + (g/partial-derivative f []))) g/derivative-symbol)) (def ^{:arglists '([f])} @@ -299,9 +299,9 @@ of `f`. The related [[emmy.env/Grad]] returns a function that produces a structure of - the opposite orientation as [[D]]. Both of these functions use forward-mode + the opposite orientation as [[D]]. Both of these functions use reverse-mode automatic differentiation." - D-forward) + D-reverse) (defn D-as-matrix [F] (fn [s] @@ -337,7 +337,7 @@ "Returns an operator that, when applied to a function `f`, produces a function that uses forward-mode automatic differentiation to compute the partial derivative of `f` at the (zero-based) slot index provided via `selectors`." - partial-forward) + partial-reverse) ;; ## Derivative Utilities ;; diff --git a/src/emmy/calculus/manifold.cljc b/src/emmy/calculus/manifold.cljc index 4a90ce64..43e5de11 100644 --- a/src/emmy/calculus/manifold.cljc +++ b/src/emmy/calculus/manifold.cljc @@ -656,7 +656,6 @@ (g/+ (g/square x) (g/square y) (g/square z)))] - (println "r is" r) (when (g/zero? r) (u/illegal-state "SphericalCylindrical singular")) (-> rep diff --git a/src/emmy/calculus/vector_calculus.cljc b/src/emmy/calculus/vector_calculus.cljc index af1947da..6ddd566d 100644 --- a/src/emmy/calculus/vector_calculus.cljc +++ b/src/emmy/calculus/vector_calculus.cljc @@ -30,10 +30,9 @@ The related [[emmy.env/D]] operator returns a function that produces a structure of the opposite orientation as [[Grad]]. Both of these functions use - forward-mode automatic differentiation." + reverse-mode automatic differentiation." (-> (fn [f] - (f/compose s/opposite - (g/partial-derivative f []))) + (f/compose s/opposite (d/D f))) (o/make-operator 'Grad))) (defn gradient diff --git a/src/emmy/dual.cljc b/src/emmy/dual.cljc index eae236a3..16d799dc 100644 --- a/src/emmy/dual.cljc +++ b/src/emmy/dual.cljc @@ -371,6 +371,13 @@ (def ^:const REVERSE-MODE ::reverse) (def ^:const REVERSE-EMPTY (->Completed {})) +;; These are here to handle cases where some collection type might see instances +;; of [[Completed]] and try and handle them during simplification. + +(defmethod g/zero? [Completed] [_] false) +(defmethod g/one? [Completed] [_] false) +(defmethod g/identity? [Completed] [_] false) + ;; `replace-tag` exists to handle subtle bugs that can arise in the case of ;; functional return values. See the "Amazing Bug" sections ;; in [[emmy.calculus.derivative-test]] for detailed examples on how this might diff --git a/src/emmy/mechanics/hamilton.cljc b/src/emmy/mechanics/hamilton.cljc index 00ecee07..4f797a8f 100644 --- a/src/emmy/mechanics/hamilton.cljc +++ b/src/emmy/mechanics/hamilton.cljc @@ -224,10 +224,9 @@ (g/zero? (g/simplify (g/determinant M)))) - (do (println "determinant" (g/determinant M)) - (throw - (ex-info "Legendre Transform Failure: determinant = 0" - {:F F :w w}))) + (throw + (ex-info "Legendre Transform Failure: determinant = 0" + {:F F :w w})) (let [v (g/solve-linear-left M (- w b))] (- (* w v) (F v))))))] (let [Dpg (D putative-G)] @@ -442,11 +441,12 @@ "p.326" [C] (fn [s] - ((- J-func - (f/compose (Phi ((D C) s)) - J-func - (Phi* ((D C) s)))) - (s/compatible-shape s)))) + (let [s-syms (s/compatible-shape s)] + ((- J-func + (f/compose (Phi ((D C) s)) + J-func + (Phi* ((D C) s)))) + s-syms)))) ;; Time-Varying code ;; diff --git a/src/emmy/numbers.cljc b/src/emmy/numbers.cljc index 5d3fab2c..6e292a4c 100644 --- a/src/emmy/numbers.cljc +++ b/src/emmy/numbers.cljc @@ -95,8 +95,8 @@ (g/infinite? a) 0 :else (g// (g/sin a) a))) -(defmethod g/sin [::v/real] [a] (Math/sin a)) -(defmethod g/cos [::v/real] [a] (Math/cos a)) +(defmethod g/sin [::v/real] [a] (if (g/zero? a) 0 (Math/sin a))) +(defmethod g/cos [::v/real] [a] (if (g/zero? a) 1 (Math/cos a))) (defmethod g/tan [::v/real] [a] (Math/tan a)) (defmethod g/cosh [::v/real] [a] (Math/cosh a)) diff --git a/test/emmy/abstract/number_test.cljc b/test/emmy/abstract/number_test.cljc index 35b6863f..28edbf62 100644 --- a/test/emmy/abstract/number_test.cljc +++ b/test/emmy/abstract/number_test.cljc @@ -153,10 +153,13 @@ values like (cos 0) can be exactly evaluated."))) (checking "inexact numbers" 100 [x sg/native-integral] - (is (= (+ 1 (Math/cos x)) - (g/+ 1 (g/cos x))) - "You get a floating-point inexact result by calling generic fns - on a number directly, by comparison.")) + (if (g/zero? x) + (is (= 2 (g/+ 1 (g/cos x))) + "special-cased exact output at 0") + (is (= (+ 1 (Math/cos x)) + (g/+ 1 (g/cos x))) + "You get a floating-point inexact result by calling generic fns + on a number directly, by comparison."))) (testing "literal-number properly prints wrapped sequences, even if they're lazy." (is (= "(* 2 x)" diff --git a/test/emmy/calculus/covariant_test.cljc b/test/emmy/calculus/covariant_test.cljc index dcc36fa9..9d0db6e9 100644 --- a/test/emmy/calculus/covariant_test.cljc +++ b/test/emmy/calculus/covariant_test.cljc @@ -48,24 +48,24 @@ (-> (simplify expr) (x/substitute '(up x0 y0 z0) 'p)))] - (is (= '(+ (* (Y↑0 p) (((partial 0) w_0) p) (X↑0 p)) - (* (Y↑0 p) (w_0 p) (((partial 0) X↑0) p)) - (* (Y↑0 p) (w_1 p) (((partial 0) X↑1) p)) - (* (Y↑0 p) (w_2 p) (((partial 0) X↑2) p)) - (* (Y↑0 p) (((partial 1) w_0) p) (X↑1 p)) - (* (Y↑0 p) (((partial 2) w_0) p) (X↑2 p)) - (* (X↑0 p) (Y↑1 p) (((partial 0) w_1) p)) - (* (X↑0 p) (Y↑2 p) (((partial 0) w_2) p)) - (* (w_0 p) (Y↑1 p) (((partial 1) X↑0) p)) - (* (w_0 p) (Y↑2 p) (((partial 2) X↑0) p)) - (* (w_1 p) (Y↑1 p) (((partial 1) X↑1) p)) - (* (w_1 p) (Y↑2 p) (((partial 2) X↑1) p)) + (is (= '(+ (* (w_2 p) (Y↑2 p) (((partial 2) X↑2) p)) (* (w_2 p) (Y↑1 p) (((partial 1) X↑2) p)) - (* (w_2 p) (Y↑2 p) (((partial 2) X↑2) p)) - (* (X↑1 p) (Y↑1 p) (((partial 1) w_1) p)) - (* (X↑1 p) (Y↑2 p) (((partial 1) w_2) p)) - (* (X↑2 p) (Y↑1 p) (((partial 2) w_1) p)) - (* (X↑2 p) (Y↑2 p) (((partial 2) w_2) p))) + (* (w_2 p) (Y↑0 p) (((partial 0) X↑2) p)) + (* (Y↑2 p) (((partial 0) w_2) p) (X↑0 p)) + (* (Y↑2 p) (w_1 p) (((partial 2) X↑1) p)) + (* (Y↑2 p) (w_0 p) (((partial 2) X↑0) p)) + (* (Y↑2 p) (((partial 1) w_2) p) (X↑1 p)) + (* (Y↑2 p) (((partial 2) w_2) p) (X↑2 p)) + (* (Y↑1 p) (X↑0 p) (((partial 0) w_1) p)) + (* (Y↑1 p) (w_1 p) (((partial 1) X↑1) p)) + (* (Y↑1 p) (w_0 p) (((partial 1) X↑0) p)) + (* (Y↑1 p) (X↑1 p) (((partial 1) w_1) p)) + (* (Y↑1 p) (X↑2 p) (((partial 2) w_1) p)) + (* (Y↑0 p) (X↑0 p) (((partial 0) w_0) p)) + (* (Y↑0 p) (w_1 p) (((partial 0) X↑1) p)) + (* (Y↑0 p) (w_0 p) (((partial 0) X↑0) p)) + (* (Y↑0 p) (X↑1 p) (((partial 1) w_0) p)) + (* (Y↑0 p) (X↑2 p) (((partial 2) w_0) p))) (present ((((g/Lie-derivative X) w) Y) R3-rect-point)))) @@ -92,14 +92,14 @@ (-> (simplify expr) (x/substitute '(up x0 y0) 'p)))] - (is (= '(+ (* -1 (Y↑0 p) (((partial 0) f) p) (((partial 0) X↑0) p)) - (* -1 (Y↑0 p) (((partial 1) f) p) (((partial 0) X↑1) p)) - (* (((partial 0) f) p) (X↑1 p) (((partial 1) Y↑0) p)) - (* (((partial 0) f) p) (((partial 0) Y↑0) p) (X↑0 p)) - (* -1 (((partial 0) f) p) (Y↑1 p) (((partial 1) X↑0) p)) - (* (X↑1 p) (((partial 1) f) p) (((partial 1) Y↑1) p)) - (* (((partial 1) f) p) (X↑0 p) (((partial 0) Y↑1) p)) - (* -1 (((partial 1) f) p) (Y↑1 p) (((partial 1) X↑1) p))) + (is (= '(+ (* (((partial 1) f) p) (((partial 0) Y↑1) p) (X↑0 p)) + (* -1 (((partial 1) f) p) (Y↑1 p) (((partial 1) X↑1) p)) + (* -1 (((partial 1) f) p) (Y↑0 p) (((partial 0) X↑1) p)) + (* (((partial 1) f) p) (((partial 1) Y↑1) p) (X↑1 p)) + (* (X↑0 p) (((partial 0) f) p) (((partial 0) Y↑0) p)) + (* -1 (Y↑1 p) (((partial 0) f) p) (((partial 1) X↑0) p)) + (* -1 (Y↑0 p) (((partial 0) f) p) (((partial 0) X↑0) p)) + (* (X↑1 p) (((partial 0) f) p) (((partial 1) Y↑0) p))) (present ((((g/Lie-derivative X) Y) f) R2-rect-point)))) @@ -114,7 +114,7 @@ ;; we only need linear terms in phi_t(x) -;; Perhaps + ;; Perhaps ;; phi_t(x) = (I + t v(I))(x) @@ -186,14 +186,14 @@ present (fn [expr] (-> (simplify expr) (x/substitute '(up q_x q_y) 'p)))] - (is (= '(+ (* -1 (Y↑0 p) (((partial 0) f) p) (((partial 0) X↑0) p)) - (* -1 (Y↑0 p) (((partial 1) f) p) (((partial 0) X↑1) p)) - (* (((partial 0) f) p) (((partial 0) Y↑0) p) (X↑0 p)) - (* -1 (((partial 0) f) p) (Y↑1 p) (((partial 1) X↑0) p)) - (* (((partial 0) f) p) (((partial 1) Y↑0) p) (X↑1 p)) - (* (((partial 1) f) p) (X↑0 p) (((partial 0) Y↑1) p)) + (is (= '(+ (* (((partial 1) f) p) (((partial 0) Y↑1) p) (X↑0 p)) (* -1 (((partial 1) f) p) (Y↑1 p) (((partial 1) X↑1) p)) - (* (((partial 1) f) p) (X↑1 p) (((partial 1) Y↑1) p))) + (* -1 (((partial 1) f) p) (Y↑0 p) (((partial 0) X↑1) p)) + (* (((partial 1) f) p) (((partial 1) Y↑1) p) (X↑1 p)) + (* (X↑0 p) (((partial 0) f) p) (((partial 0) Y↑0) p)) + (* -1 (Y↑1 p) (((partial 0) f) p) (((partial 1) X↑0) p)) + (* -1 (Y↑0 p) (((partial 0) f) p) (((partial 0) X↑0) p)) + (* (X↑1 p) (((partial 0) f) p) (((partial 1) Y↑0) p))) (present ((D (fn [t] (- ((Y f) ((phiX t) m_0)) @@ -207,7 +207,7 @@ ((Y (compose f (phiX t))) m_0)))) 0))))) -(is (= 0 (simplify + (is (= 0 (simplify (- result-via-Lie ((D (fn [t] (- ((Y f) ((phiX t) m_0)) @@ -750,7 +750,7 @@ (is (= '(down (+ (* -1 (cos (theta t)) (sin (theta t)) (expt ((D phi) t) 2)) (((expt D 2) theta) t)) - (+ (* 2 (cos (theta t)) ((D theta) t) (sin (theta t)) ((D phi) t)) + (+ (* 2 (cos (theta t)) (sin (theta t)) ((D phi) t) ((D theta) t)) (* (expt (sin (theta t)) 2) (((expt D 2) phi) t)))) (g/freeze (simplify @@ -813,7 +813,7 @@ ;; So \Gammar_{\theta \theta} = -r, \Gamma\theta_{\theta \theta} = 0 ;; ;; These are correct Christoffel symbols... -))))) + ))))) (defn CD "Computation of Covariant derivatives by difference quotient. [[CD]] is parallel diff --git a/test/emmy/calculus/form_field_test.cljc b/test/emmy/calculus/form_field_test.cljc index 2385de00..1f26fcbb 100644 --- a/test/emmy/calculus/form_field_test.cljc +++ b/test/emmy/calculus/form_field_test.cljc @@ -331,18 +331,18 @@ ((((g/square ff/d) (m/literal-scalar-field 'f R3-rect)) X Y) R3-cyl-point)))) - (is (= '(+ (* (X↑0 p) (Y↑1 p) (((partial 0) w_1) p)) - (* -1 (X↑0 p) (Y↑1 p) (((partial 1) w_0) p)) - (* (X↑0 p) (Y↑2 p) (((partial 0) w_2) p)) + (is (= '(+ (* (X↑0 p) (Y↑2 p) (((partial 0) w_2) p)) (* -1 (X↑0 p) (Y↑2 p) (((partial 2) w_0) p)) - (* -1 (X↑1 p) (((partial 0) w_1) p) (Y↑0 p)) - (* (X↑1 p) (((partial 1) w_0) p) (Y↑0 p)) + (* (X↑0 p) (Y↑1 p) (((partial 0) w_1) p)) + (* -1 (X↑0 p) (Y↑1 p) (((partial 1) w_0) p)) (* (X↑1 p) (Y↑2 p) (((partial 1) w_2) p)) (* -1 (X↑1 p) (Y↑2 p) (((partial 2) w_1) p)) - (* -1 (X↑2 p) (Y↑1 p) (((partial 1) w_2) p)) - (* (X↑2 p) (Y↑1 p) (((partial 2) w_1) p)) + (* -1 (X↑1 p) (((partial 0) w_1) p) (Y↑0 p)) + (* (X↑1 p) (((partial 1) w_0) p) (Y↑0 p)) (* -1 (X↑2 p) (((partial 0) w_2) p) (Y↑0 p)) - (* (X↑2 p) (((partial 2) w_0) p) (Y↑0 p))) + (* (X↑2 p) (((partial 2) w_0) p) (Y↑0 p)) + (* -1 (X↑2 p) (Y↑1 p) (((partial 1) w_2) p)) + (* (X↑2 p) (Y↑1 p) (((partial 2) w_1) p))) (-> (((ff/d w) X Y) R3-rect-point) (simplify) (x/substitute '(up x0 y0 z0) 'p)))) @@ -353,24 +353,24 @@ (ff/wedge dy dz)) (* (m/literal-scalar-field 'omega_2 R3-rect) (ff/wedge dz dx)))] - (is (= '(+ (* (X↑0 p) (Y↑1 p) (Z↑2 p) (((partial 0) omega_1) p)) - (* (X↑0 p) (Y↑1 p) (Z↑2 p) (((partial 1) omega_2) p)) - (* (X↑0 p) (Y↑1 p) (Z↑2 p) (((partial 2) omega_0) p)) - (* -1 (X↑0 p) (Y↑2 p) (((partial 0) omega_1) p) (Z↑1 p)) - (* -1 (X↑0 p) (Y↑2 p) (((partial 1) omega_2) p) (Z↑1 p)) - (* -1 (X↑0 p) (Y↑2 p) (((partial 2) omega_0) p) (Z↑1 p)) + (is (= '(+ (* -1 (X↑0 p) (Y↑2 p) (Z↑1 p) (((partial 0) omega_1) p)) + (* -1 (X↑0 p) (Y↑2 p) (Z↑1 p) (((partial 1) omega_2) p)) + (* -1 (X↑0 p) (Y↑2 p) (Z↑1 p) (((partial 2) omega_0) p)) + (* (X↑0 p) (Y↑1 p) (((partial 0) omega_1) p) (Z↑2 p)) + (* (X↑0 p) (Y↑1 p) (((partial 1) omega_2) p) (Z↑2 p)) + (* (X↑0 p) (Y↑1 p) (((partial 2) omega_0) p) (Z↑2 p)) (* (X↑1 p) (Y↑2 p) (((partial 0) omega_1) p) (Z↑0 p)) (* (X↑1 p) (Y↑2 p) (((partial 1) omega_2) p) (Z↑0 p)) (* (X↑1 p) (Y↑2 p) (((partial 2) omega_0) p) (Z↑0 p)) - (* -1 (X↑1 p) (Y↑0 p) (Z↑2 p) (((partial 0) omega_1) p)) - (* -1 (X↑1 p) (Y↑0 p) (Z↑2 p) (((partial 1) omega_2) p)) - (* -1 (X↑1 p) (Y↑0 p) (Z↑2 p) (((partial 2) omega_0) p)) + (* -1 (X↑1 p) (Y↑0 p) (((partial 0) omega_1) p) (Z↑2 p)) + (* -1 (X↑1 p) (Y↑0 p) (((partial 1) omega_2) p) (Z↑2 p)) + (* -1 (X↑1 p) (Y↑0 p) (((partial 2) omega_0) p) (Z↑2 p)) (* -1 (X↑2 p) (Y↑1 p) (((partial 0) omega_1) p) (Z↑0 p)) (* -1 (X↑2 p) (Y↑1 p) (((partial 1) omega_2) p) (Z↑0 p)) (* -1 (X↑2 p) (Y↑1 p) (((partial 2) omega_0) p) (Z↑0 p)) - (* (X↑2 p) (Y↑0 p) (((partial 0) omega_1) p) (Z↑1 p)) - (* (X↑2 p) (Y↑0 p) (((partial 1) omega_2) p) (Z↑1 p)) - (* (X↑2 p) (Y↑0 p) (((partial 2) omega_0) p) (Z↑1 p))) + (* (X↑2 p) (Y↑0 p) (Z↑1 p) (((partial 0) omega_1) p)) + (* (X↑2 p) (Y↑0 p) (Z↑1 p) (((partial 1) omega_2) p)) + (* (X↑2 p) (Y↑0 p) (Z↑1 p) (((partial 2) omega_0) p))) (-> (((ff/d omega) X Y Z) R3-rect-point) (simplify) (x/substitute '(up x0 y0 z0) 'p)))) diff --git a/test/emmy/calculus/map_test.cljc b/test/emmy/calculus/map_test.cljc index c12113f7..a47fb673 100644 --- a/test/emmy/calculus/map_test.cljc +++ b/test/emmy/calculus/map_test.cljc @@ -90,10 +90,10 @@ and the differentials of coordinate functions." μ (m/literal-manifold-map 'μ R1-rect R2-rect) f (man/literal-manifold-function 'f R2-rect)] - (is (= '(+ (* (((partial 0) f) (up (μ↑0 τ) (μ↑1 τ))) - ((D μ↑0) τ)) - (* (((partial 1) f) (up (μ↑0 τ) (μ↑1 τ))) - ((D μ↑1) τ))) + (is (= '(+ (* (((partial 1) f) (up (μ↑0 τ) (μ↑1 τ))) + ((D μ↑1) τ)) + (* (((partial 0) f) (up (μ↑0 τ) (μ↑1 τ))) + ((D μ↑0) τ))) (simplify ((((m/differential μ) d:dt) f) ((man/point R1-rect) 'τ))))) @@ -120,10 +120,10 @@ and the differentials of coordinate functions." ;; "However, if we kludge the correct argument it gives the expected ;; answer." - (is (= '(/ (+ (* ((D μ↑0) τ) (e1↑1 (up x0 y0))) - (* -1 ((D μ↑1) τ) (e1↑0 (up x0 y0)))) - (+ (* (e1↑1 (up x0 y0)) (e0↑0 (up x0 y0))) - (* -1 (e1↑0 (up x0 y0)) (e0↑1 (up x0 y0))))) + (is (= '(/ (+ (* -1 ((D μ↑1) τ) (e1↑0 (up x0 y0))) + (* ((D μ↑0) τ) (e1↑1 (up x0 y0)))) + (+ (* -1 (e1↑0 (up x0 y0)) (e0↑1 (up x0 y0))) + (* (e1↑1 (up x0 y0)) (e0↑0 (up x0 y0))))) (simplify (((nth edual 0) (vf/procedure->vector-field @@ -142,8 +142,8 @@ and the differentials of coordinate functions." (man/chart R1-rect)) f (f/compose (af/literal-function 'f '(-> (UP Real Real) Real)) (man/chart S2-spherical))] - (is (= '(+ (* (((partial 0) f) (up (θ τ) (φ τ))) ((D θ) τ)) - (* (((partial 1) f) (up (θ τ) (φ τ))) ((D φ) τ))) + (is (= '(+ (* (((partial 1) f) (up (θ τ) (φ τ))) ((D φ) τ)) + (* (((partial 0) f) (up (θ τ) (φ τ))) ((D θ) τ))) (simplify ((((m/differential μ) d:dt) f) ((man/point R1-rect) 'τ))))) diff --git a/test/emmy/fdg/ch4_test.cljc b/test/emmy/fdg/ch4_test.cljc index 46be8a60..3414004d 100644 --- a/test/emmy/fdg/ch4_test.cljc +++ b/test/emmy/fdg/ch4_test.cljc @@ -178,29 +178,28 @@ ;; equations" alluded to on p.48 of FDG. By equating corresponding ;; entries, we may verify the solution of a, b, c given there. (is (= '(matrix-by-rows - [(+ (* a (sin phi) (sin psi) (sin theta)) - (* -1 b (sin psi) (cos phi) (cos theta)) + [(+ (* a (sin psi) (sin phi) (sin theta)) + (* -1 b (sin psi) (cos theta) (cos phi)) (* -1 c (sin phi) (cos theta) (cos psi)) (* -1 b (sin phi) (cos psi)) (* -1 c (sin psi) (cos phi))) - (+ - (* a (sin phi) (sin theta) (cos psi)) - (* -1 b (cos phi) (cos theta) (cos psi)) - (* c (sin phi) (sin psi) (cos theta)) - (* b (sin phi) (sin psi)) - (* -1 c (cos phi) (cos psi))) + (+ (* a (sin phi) (sin theta) (cos psi)) + (* -1 b (cos theta) (cos phi) (cos psi)) + (* c (sin psi) (sin phi) (cos theta)) + (* b (sin psi) (sin phi)) + (* -1 c (cos phi) (cos psi))) (+ (* a (sin phi) (cos theta)) (* b (sin theta) (cos phi)))] [(+ (* -1 a (sin psi) (sin theta) (cos phi)) - (* -1 b (sin phi) (sin psi) (cos theta)) - (* c (cos phi) (cos theta) (cos psi)) + (* -1 b (sin psi) (sin phi) (cos theta)) + (* c (cos theta) (cos phi) (cos psi)) (* b (cos phi) (cos psi)) - (* -1 c (sin phi) (sin psi))) + (* -1 c (sin psi) (sin phi))) (+ (* -1 a (sin theta) (cos phi) (cos psi)) (* -1 b (sin phi) (cos theta) (cos psi)) - (* -1 c (sin psi) (cos phi) (cos theta)) + (* -1 c (sin psi) (cos theta) (cos phi)) (* -1 b (sin psi) (cos phi)) (* -1 c (sin phi) (cos psi))) - (+ (* -1 a (cos phi) (cos theta)) (* b (sin phi) (sin theta)))] + (+ (* -1 a (cos theta) (cos phi)) (* b (sin phi) (sin theta)))] [(+ (* a (sin psi) (cos theta)) (* c (sin theta) (cos psi))) (+ (* a (cos theta) (cos psi)) (* -1 c (sin psi) (sin theta))) (* -1 a (sin theta))]) @@ -211,8 +210,8 @@ [(* -1 (sin psi) (sin theta)) (* -1 (sin theta) (cos psi)) (* -1 (cos theta))] - [(+ (* (sin psi) (cos phi) (cos theta)) (* (sin phi) (cos psi))) - (+ (* (cos phi) (cos theta) (cos psi)) (* -1 (sin phi) (sin psi))) + [(+ (* (sin psi) (cos theta) (cos phi)) (* (sin phi) (cos psi))) + (+ (* (cos theta) (cos phi) (cos psi)) (* -1 (sin psi) (sin phi))) (* -1 (sin theta) (cos phi))]) (simplify ((D h) 0)))))))) diff --git a/test/emmy/fdg/ch7_test.cljc b/test/emmy/fdg/ch7_test.cljc index 4ddd0d0a..7c8e7d1b 100644 --- a/test/emmy/fdg/ch7_test.cljc +++ b/test/emmy/fdg/ch7_test.cljc @@ -350,9 +350,9 @@ (((expt D 2) alpha) t)) (+ (* 2 (cos (alpha t)) - ((D alpha) t) (sin (alpha t)) - ((D beta) t)) + ((D beta) t) + ((D alpha) t)) (* (expt (sin (alpha t)) 2) (((expt D 2) beta) t)))) (simplify diff --git a/test/emmy/fdg/ch9_test.cljc b/test/emmy/fdg/ch9_test.cljc index 9bb073db..4738b963 100644 --- a/test/emmy/fdg/ch9_test.cljc +++ b/test/emmy/fdg/ch9_test.cljc @@ -147,9 +147,9 @@ (let [q (up x y)] (is (= '(down - (+ (* ((D x) (f t)) (m_00 (up (x (f t)) (y (f t)))) (((expt D 2) f) t)) - (* (m_01 (up (x (f t)) (y (f t)))) ((D y) (f t)) (((expt D 2) f) t))) - (+ (* ((D x) (f t)) (m_01 (up (x (f t)) (y (f t)))) (((expt D 2) f) t)) + (+ (* (m_01 (up (x (f t)) (y (f t)))) ((D y) (f t)) (((expt D 2) f) t)) + (* ((D x) (f t)) (m_00 (up (x (f t)) (y (f t)))) (((expt D 2) f) t))) + (+ (* (m_01 (up (x (f t)) (y (f t)))) ((D x) (f t)) (((expt D 2) f) t)) (* ((D y) (f t)) (m_11 (up (x (f t)) (y (f t)))) (((expt D 2) f) t)))) (simplify ((- (compose (e/Euler-Lagrange-operator L2) @@ -188,10 +188,10 @@ (is (= '(/ (+ (* (expt c 2) (((expt (partial 0) 2) V) (up y z t))) (* (expt c 2) (((expt (partial 1) 2) V) (up y z t))) (* (expt c 2) (((expt (partial 2) 2) V) (up y z t))) - (* -1 (expt (((partial 0) V) (up y z t)) 2)) (* 2 (V (up y z t)) (((expt (partial 0) 2) V) (up y z t))) (* 2 (V (up y z t)) (((expt (partial 1) 2) V) (up y z t))) (* 2 (V (up y z t)) (((expt (partial 2) 2) V) (up y z t))) + (* -1 (expt (((partial 0) V) (up y z t)) 2)) (* -1 (expt (((partial 1) V) (up y z t)) 2)) (* -1 (expt (((partial 2) V) (up y z t)) 2))) (+ (expt c 2) (* 2N (V (up y z t))))) diff --git a/test/emmy/fdg/einstein_test.cljc b/test/emmy/fdg/einstein_test.cljc index b2c1a26f..a2c87c7f 100644 --- a/test/emmy/fdg/einstein_test.cljc +++ b/test/emmy/fdg/einstein_test.cljc @@ -193,7 +193,7 @@ (e/Christoffel->Cartan (e/metric->Christoffel-2 metric basis))) ws (e/basis->oneform-basis basis)] - (is (= '[(/ (+ (* (expt c 2) ((D rho) t) (R t)) + (is (= '[(/ (+ (* (expt c 2) (R t) ((D rho) t)) (* 3 (expt c 2) ((D R) t) (rho t)) (* 3 ((D R) t) (p t))) (* (expt c 2) (R t))) diff --git a/test/emmy/mechanics/hamilton_test.cljc b/test/emmy/mechanics/hamilton_test.cljc index d004c355..be01c8a1 100644 --- a/test/emmy/mechanics/hamilton_test.cljc +++ b/test/emmy/mechanics/hamilton_test.cljc @@ -697,38 +697,38 @@ ((H/time-independent-canonical? (H/F->CT L/p->r)) (up 't (L/coordinate-tuple 'r 'phi) - (L/momentum-tuple 'p_r 'p_phi)))))) - - ;; but not all transforms are canonical: - (testing "non-canonical-transform" - (with-redefs [gensym (a/monotonic-symbol-generator 1 "x")] - (letfn [(a-non-canonical-transform [[t theta p]] - (let [x (* p (g/sin theta)) - p_x (* p (g/cos theta))] - (up t x p_x)))] - (is (= '(up 0 - (+ (* -1 p x3) x3) - (+ (* p x2) (* -1 x2))) - (simplify - ((H/time-independent-canonical? a-non-canonical-transform) - (up 't 'theta 'p)))) - "The genysym redef is ugly but makes this test repeatable, at + (L/momentum-tuple 'p_r 'p_phi))))))) + + ;; but not all transforms are canonical: + (testing "non-canonical-transform" + (with-redefs [gensym (a/monotonic-symbol-generator 2 "x")] + (letfn [(a-non-canonical-transform [[t theta p]] + (let [x (* p (g/sin theta)) + p_x (* p (g/cos theta))] + (up t x p_x)))] + (is (= '(up 0 + (+ (* -1 p x03) x03) + (+ (* p x02) (* -1 x02))) + (simplify + ((H/time-independent-canonical? a-non-canonical-transform) + (up 't 'theta 'p)))) + "The genysym redef is ugly but makes this test repeatable, at least.") - (is (= '(matrix-by-rows - [0 0 0] - [0 0 (+ (* -1 p) 1)] - [0 (+ p -1) 0]) - (simplify - ((H/symplectic? a-non-canonical-transform) - ['t 'theta 'p])))) - - (is (= '(matrix-by-rows - [0 (+ (* -1 p) 1)] - [(+ p -1) 0]) - (simplify - ((H/symplectic-transform? a-non-canonical-transform) - ['t 'theta 'p]))))))))) + (is (= '(matrix-by-rows + [0 0 0] + [0 0 (+ (* -1 p) 1)] + [0 (+ p -1) 0]) + (simplify + ((H/symplectic? a-non-canonical-transform) + ['t 'theta 'p])))) + + (is (= '(matrix-by-rows + [0 (+ (* -1 p) 1)] + [(+ p -1) 0]) + (simplify + ((H/symplectic-transform? a-non-canonical-transform) + ['t 'theta 'p])))))))) (deftest symplectic-tests (testing "unit" diff --git a/test/emmy/mechanics/rigid_test.cljc b/test/emmy/mechanics/rigid_test.cljc index 6c1a7785..ebd3573f 100644 --- a/test/emmy/mechanics/rigid_test.cljc +++ b/test/emmy/mechanics/rigid_test.cljc @@ -24,10 +24,10 @@ (deftest rigid-tests (f/with-literal-functions [theta phi psi] (is (= '(column-matrix - (+ (* (sin (theta t)) ((D phi) t) (sin (psi t))) - (* (cos (psi t)) ((D theta) t))) - (+ (* (cos (psi t)) (sin (theta t)) ((D phi) t)) - (* -1 ((D theta) t) (sin (psi t)))) + (+ (* (sin (psi t)) (sin (theta t)) ((D phi) t)) + (* ((D theta) t) (cos (psi t)))) + (+ (* (sin (theta t)) ((D phi) t) (cos (psi t))) + (* -1 (sin (psi t)) ((D theta) t))) (+ (* (cos (theta t)) ((D phi) t)) ((D psi) t))) (simplify ((rig/Euler->omega-body @@ -35,12 +35,11 @@ 't)))) (is (= '(column-matrix - (+ (* (sin (theta t)) ((D phi) t) (sin (psi t))) - (* (cos (psi t)) ((D theta) t))) - (+ (* (cos (psi t)) (sin (theta t)) ((D phi) t)) - (* -1 ((D theta) t) (sin (psi t)))) - (+ (* (cos (theta t)) ((D phi) t)) - ((D psi) t))) + (+ (* (sin (psi t)) (sin (theta t)) ((D phi) t)) + (* ((D theta) t) (cos (psi t)))) + (+ (* (sin (theta t)) ((D phi) t) (cos (psi t))) + (* -1 (sin (psi t)) ((D theta) t))) + (+ (* (cos (theta t)) ((D phi) t)) ((D psi) t))) (simplify (((rig/M-of-q->omega-body-of-t rotation/Euler->M) (up theta phi psi)) @@ -60,10 +59,10 @@ (let [an-Euler-state (up 't (up 'theta 'phi 'psi) (up 'thetadot 'phidot 'psidot))] - (is (= '(+ (* A phidot (expt (sin psi) 2) (expt (sin theta) 2)) + (is (= '(+ (* A phidot (expt (sin theta) 2) (expt (sin psi) 2)) (* B phidot (expt (sin theta) 2) (expt (cos psi) 2)) - (* A thetadot (sin psi) (sin theta) (cos psi)) - (* -1 B thetadot (sin psi) (sin theta) (cos psi)) + (* A thetadot (sin theta) (sin psi) (cos psi)) + (* -1N B thetadot (sin theta) (sin psi) (cos psi)) (* C phidot (expt (cos theta) 2)) (* C psidot (cos theta))) (simplify diff --git a/test/emmy/series_test.cljc b/test/emmy/series_test.cljc index ba40b18d..ae68af0a 100644 --- a/test/emmy/series_test.cljc +++ b/test/emmy/series_test.cljc @@ -666,10 +666,10 @@ (deftest cos-1-square-terms (testing "Generate terms for the power series expansion of $cos(z)-1$." - ;; (take 10 (g/- s/cos-series 1)) - ;; => (0 0 -1/2 0 1/24 0 -1/720 0 1/40320 0) - ;; From this we can see that we may regard the expansion as a series in x^2, with a - ;; zero constant term. + ;; (take 10 (g/- s/cos-series 1)) + ;; => (0 0 -1/2 0 1/24 0 -1/720 0 1/40320 0) + ;; From this we can see that we may regard the expansion as a series in x^2, with a + ;; zero constant term. (let [terms (->> (g/- s/cos-series 1) (remove g/zero?) ;; eliminate the useless zero coefficients of the odd powers of x (map double) ;; we don't want the rational arithmetic to survive diff --git a/test/emmy/sicm/ch2_test.cljc b/test/emmy/sicm/ch2_test.cljc index f4a5b9e4..532266be 100644 --- a/test/emmy/sicm/ch2_test.cljc +++ b/test/emmy/sicm/ch2_test.cljc @@ -50,18 +50,18 @@ (((r/M-of-q->omega-body-of-t Euler->M) q) 't))))) (is (= '(column-matrix - (+ (* φdot (sin θ) (sin ψ)) (* θdot (cos ψ))) - (+ (* φdot (cos ψ) (sin θ)) (* -1 θdot (sin ψ))) + (+ (* φdot (sin ψ) (sin θ)) (* θdot (cos ψ))) + (+ (* φdot (sin θ) (cos ψ)) (* -1 θdot (sin ψ))) (+ (* φdot (cos θ)) ψdot)) (e/freeze (simplify ((r/M->omega-body Euler->M) Euler-state)))))))) (deftest section-2-9 - (is (v/= '(+ (* A φdot (expt (sin ψ) 2) (expt (sin θ) 2)) + (is (v/= '(+ (* A φdot (expt (sin θ) 2) (expt (sin ψ) 2)) (* B φdot (expt (sin θ) 2) (expt (cos ψ) 2)) - (* A θdot (sin ψ) (sin θ) (cos ψ)) - (* -1 B θdot (sin ψ) (sin θ) (cos ψ)) + (* A θdot (sin θ) (sin ψ) (cos ψ)) + (* -1N B θdot (sin θ) (sin ψ) (cos ψ)) (* C φdot (expt (cos θ) 2)) (* C ψdot (cos θ))) (simplify