emmy

Author: sritchie

deps.edn

@@ -16,8 +16,8 @@
         ;; some statistical routines
         kixi/stats {:mvn/version "0.5.4"}
 
-        ;; SICMUtils numerical and physics routines
-        sicmutils/sicmutils {:mvn/version "0.22.0"}
+        ;; Emmy numerical and physics routines
+        io.github.mentat-collective/emmy {:git/sha "b98fef51d80d3edcff3100562b929f9980ff34d7"}
 
         ;; semantic web goodies and box/arrow graphs
         io.github.jackrusher/mundaneum {:git/sha "d2c934a12388d88ddb3e53fef92ec2eef97d6140"}
@@ -44,7 +44,7 @@
                                                    "notebooks/images.clj"
                                                    "notebooks/logo.clj"
                                                    "notebooks/semantic.clj"
-                                                   "notebooks/sicmutils.clj"
+                                                   "notebooks/emmy.clj"
                                                    "notebooks/rule_30.clj"
                                                    "notebooks/zipper_with_scars.clj"]}
                                :main-opts ["-m" "babashka.cli.exec"]}

dev/user.clj

@@ -17,7 +17,7 @@
   ;; or call `clerk/show!` explicitly
   (clerk/show! "notebooks/introduction.clj")
   (clerk/show! "notebooks/data_science.clj")
-  (clerk/show! "notebooks/sicmutils.clj")
+  (clerk/show! "notebooks/emmy.clj")
   (clerk/show! "notebooks/rule_30.clj")
   (clerk/show! "notebooks/semantic.clj")
   (clerk/show! "notebooks/images.clj")

index.md

@@ -39,7 +39,7 @@
         :description "Some examples about dealing with images."}
        {:title "The double pendulum"
         :preview "https://cdn.nextjournal.com/data/Qmdhk9WEogAvH9cgehnpq3REkATa91JiWgYNUtoSW74Q9A?filename=CleanShot%202021-11-30%20at%[email protected]&content-type=image/png"
-        :path "notebooks/sicmutils.clj"
+        :path "notebooks/emmy.clj"
         :description "Simulate and visualize physical systems from the REPL."}
        {:title "Rule 30"
         :preview "https://cdn.nextjournal.com/data/QmQCrqkdYtKfNm9CGbXhzY4cy6qG8xhpWaCRPF5m6biLgV?filename=CleanShot%202021-11-30%20at%[email protected]&content-type=image/png"

notebooks/sicmutils.clj renamed to notebooks/emmy.clj

@@ -4,16 +4,18 @@
 ;; 1.44](https://tgvaughan.github.io/sicm/chapter001.html#Exe_1-44) from Sussman
 ;; and Wisdom's [Structure and Interpretation of Classical
 ;; Mechanics](https://tgvaughan.github.io/sicm/), using
-;; the [SICMUtils](https://github.com/sicmutils/sicmutils) Clojure library and
+;; the [Emmy](https://github.com/mentat-collective/emmy) Clojure library and
 ;; the Clerk rendering environment.
 
 #_{:clj-kondo/ignore [:refer-all]}
-(ns sicmutils
+(ns emmy
   (:refer-clojure
-   :exclude [+ - * / partial ref zero? numerator denominator compare = run!])
+   :exclude [+ - * / partial ref zero? numerator denominator compare = run!
+             abs infinite?])
   (:require [nextjournal.clerk :as clerk]
-            [sicmutils.env :as e :refer :all]
-            [sicmutils.expression.render :as xr]))
+            [emmy.env :as e :refer :all]
+            [emmy.expression.compile :as xc]
+            [emmy.expression.render :as xr]))
 
 ;; ## Lagrangian
 ;;
@@ -39,7 +41,6 @@
        (* 1/2 m2 (+ (square xdot2)
                     (square ydot2))))))
 
-
 ;; `V` describes a uniform gravitational potential with coefficient `g`, acting
 ;; on two particles with masses of, respectively, `m1` and `m2`. Again, this is
 ;; written in rectangular coordinates.
@@ -86,6 +87,8 @@
 
 ;; And here are the equations of motion for the system:
 
+;; TODO this currently causes a notebook failure.
+#_
 (let [L (L-double-pendulum 'm_1 'm_2 'l_1 'l_2 'g)]
   (binding [xr/*TeX-vertical-down-tuples* true]
     (render-eq
@@ -172,7 +175,8 @@
 ;; state. Chaotic first:
 
 (def raw-chaotic-data
-  (run! step horizon chaotic-initial-q))
+  (time
+   (run! step horizon chaotic-initial-q)))
 
 ;; Looks good:
 
@@ -181,12 +185,14 @@
 ;; Next, the regular initial condition:
 
 (def raw-regular-data
-  (run! step horizon regular-initial-q))
+  (time
+   (run! step horizon regular-initial-q)))
 
 ;; Peek at the final state:
 
 (peek raw-regular-data)
 
+
 ;; ## Measurements, Data Transformation
 
 ;; Next we'll chart the measurements trapped in those sequences of state tuples.
@@ -225,8 +231,16 @@
 #_{:clj-kondo/ignore [:unresolved-symbol]}
 (defn transform-data [xs]
   (let [energy-fn (L-energy m1 m2 l1 l2 g)
-        monitor   (energy-monitor energy-fn (first xs))
-        xform     (angles->rect l1 l2)
+        monitor   (xc/compile-state-fn
+                   (energy-monitor energy-fn (first xs))
+                   false
+                   (first xs)
+                   {:calling-convention :structure})
+        xform     (xc/compile-state-fn
+                   (angles->rect l1 l2)
+                   false
+                   (first xs)
+                   {:calling-convention :structure})
         pv        (principal-value Math/PI)]
     (map (fn [[t [theta1 theta2] [thetadot1 thetadot2] :as state]]
            (let [[x1 y1 x2 y2] (xform state)]
@@ -245,11 +259,11 @@
 ;; The following forms transform the raw data for each initial condition and
 ;; bind the results to `chaotic-data` and `regular-data` for exploration.
 
-(def chaotic-data
+(defonce chaotic-data
   (doall
    (transform-data raw-chaotic-data)))
 
-(def regular-data
+(defonce regular-data
   (doall
    (transform-data raw-regular-data)))
 
@@ -307,9 +321,9 @@
 ;; a helper function that should be in clojure.core
 (defn deep-merge [v & vs]
   (letfn [(rec-merge [v1 v2]
-                     (if (and (map? v1) (map? v2))
-                       (merge-with deep-merge v1 v2)
-                       v2))]
+            (if (and (map? v1) (map? v2))
+              (merge-with deep-merge v1 v2)
+              v2))]
     (when (some identity vs)
       (reduce #(rec-merge %1 %2) v vs))))
 
@@ -439,9 +453,9 @@
 (defn L-double-double-pendulum [m1 m2 l1 l2 g]
   (fn [[t [thetas1 thetas2] [qdots1 qdots2]]]
     (let [s1 (up t thetas1 qdots1)
-            s2 (up t thetas2 qdots2)]
-        (+ ((L-double-pendulum m1 m2 l1 l2 g) s1)
-           ((L-double-pendulum m1 m2 l1 l2 g) s2)))))
+          s2 (up t thetas2 qdots2)]
+      (+ ((L-double-pendulum m1 m2 l1 l2 g) s1)
+         ((L-double-pendulum m1 m2 l1 l2 g) s2)))))
 
 (def dd-state-derivative
   (compose

notebooks/emmy_repro.clj

@@ -0,0 +1,37 @@
+(ns emmy-repro
+  (:require [emmy.env :as e]
+            [emmy.mechanics.lagrange]))
+
+;; ## BUG 1:
+
+;; This notebook takes close to 2 seconds to evaluate:
+
+;; Clerk evaluated '/Users/sritchie/code/clj/clerk-demo/notebooks/emmy_repro.clj' in 1853.674042ms.
+
+(defn angles->rect [theta1]
+  (e/sin theta1))
+
+
+;; Form the final Langrangian in generalized coordinates (the angles of each
+;; segment) by composing `L-rect` with a properly transformed `angles->rect`
+;; coordinate transform!
+
+;; ## BUG 2:
+;;
+;; The following form:
+
+#_
+(let [L (emmy.mechanics.lagrange/L-pendulum 'g 'm 'l)]
+  (((e/Lagrange-equations L)
+    (e/literal-function 'theta_1))
+   't))
+
+;; Evaluates to this:
+(e/literal-number
+ '(- (* 1/2 m 2 l (((expt D 2) theta_1) t) l) (* g m l (- (sin (theta_1 t))))))
+
+
+;; But if I include it in a notebook, I get this:
+
+;; Execution error (NullPointerException) at clojure.tools.analyzer.jvm.utils/members* (utils.clj:272).
+;; Cannot invoke "java.lang.Class.getName()" because the return value of "clojure.lang.IFn.invoke(Object)" is null

notebooks/index.md

@@ -3,6 +3,6 @@
 ## 🥁 [Intro](#/notebooks/introduction.clj)
 ## 📈 [Data Science](#/notebooks/data_science.clj)
 ## 🕸 [Semantic Queries](#/notebooks/semantic.clj)
-## 🔬 [SICMUtils Double Pendulum](#/notebooks/sicmutils.clj)
+## 🔬 [Emmy Double Pendulum](#/notebooks/emmy.clj)
 ## 🕹 [Rule 30](#/notebooks/rule_30.clj)
 ## ✏️ [Markdown Files](#/notebooks/markdown.md)