mentat-collective / emmy Goto Github PK
View Code? Open in Web Editor NEWThe Emmy Computer Algebra System.
Home Page: https://emmy.mentat.org
License: GNU General Public License v3.0
The Emmy Computer Algebra System.
Home Page: https://emmy.mentat.org
License: GNU General Public License v3.0
If we do this, and take a tip from clojure.algo.generic
, it will be easier to write generic functions that do stuff like replace-tag
from sicmutils.differential
.
Many types implement something like this, thanks to @littleredcomputer . Getting them together under one roof is a nice move.
The Quadrature interface contributed in sicmutils/sicmutils#129 lets you specify in the options map whether your interval is open, closed or half-open.
Certain integrators can only handle certain cases - the Midpoint and Milne rules, or the open version of Romberg, handle open cases, while all of the trapezoid-based methods only handle closed intervals.
"Adaptive" can take an open and closed integrator and make a new integrator that can handle any of these combinations.
We might modify the defintegrator
interface in sicmutils.numerical.quadrature.common
to take a Clojure Spec declaring what kinds of intervals the defined integrator can handle.
We could also make a simpler version of adaptive that only splits a single time, to allow integrators like Romberg to handle all of the kinds of intervals (open-closed would split at an internal midpoint and attack the left and right with fully open and fully closed, respectively.)
We have univariate power series. Here's some prior art on multivariate power series: https://github.com/michibo/pwrsrs
For example, derivative becomes:
(r/choice
(r/rule (D ?f) r/=> ((expt D 2) ?f))
(r/rule ((expt D ?n) ?f) r/=> ((expt D (? #(inc (% '?n)))) ?f))
(r/rule ?f r/=> (D ?f)))
This is faster and more readable!
Here's a nice library that does this now: https://github.com/Kerbaya/ieee754lib
Related Spire issue: typelevel/spire#236
In scmutils, this lives in numerics/quad/quad-flo.scm
:
;;; Fundamental quad manipulation
;;; given double arguments.
;;; Knuth's theorem B provides a method of adding 2 floats
;;; to produce a sum and an error estimate.
;;; It takes 6 floating operations and no comparisions.
(define (+B u v k) ;k=(lambda (sum error) ...)
(let ((s (flo:+ u v)))
(let ((up (flo:- s v)))
(let ((vpp (flo:- s up)))
(k s (flo:+ (flo:- u up) (flo:- v vpp)))))))
(define (-B u v k) ;k=(lambda (sum error) ...)
(let ((s (flo:- u v)))
(let ((up (flo:+ s v)))
(let ((vpp (flo:- s up)))
(k s
(flo:- (flo:- u up)
(flo:+ v vpp)))))))
;;; Multiplication of two floating operations, transcribed from
;;; Quinn & Tremaine 1990.
(define IEEE-P 53) ;see standard.
(define tremaine-j
(ceiling (/ IEEE-p 2)))
(define tremaine-c
(exact->inexact (+ (expt 2 tremaine-j) 1)))
(define (*T u v k) ;k=(lambda (prod error) ...)
(let* ((up (flo:* u tremaine-c))
(u1 (flo:- u up))
(u1 (flo:+ u1 up))
(u2 (flo:- u u1))
(vp (flo:* v tremaine-c))
(v1 (flo:- v vp))
(v1 (flo:+ v1 vp))
(v2 (flo:- v v1))
(prod (flo:* u v))
(e1 (flo:* u1 v1))
(e2 (flo:* u2 v2))
(e12 (flo:* u1 v2))
(e21 (flo:* u2 v1))
(err (flo:- e1 prod))
(err (flo:+ err e12))
(err (flo:+ err e21))
(perr (flo:+ err e2)))
(k prod perr)))
(define (/T u v k) ;k=(lambda (prod error) ...)
(let ((quot (flo:/ u v)))
(*T quot v
(lambda (rat err)
(let* ((qerr (flo:- u rat))
(qerr (flo:/ (flo:- qerr err) v)))
(k quot qerr))))))
;;; Quad arithmetic -- a quad is a cons of the high and low part.
(define (quad:+ q1 q2)
(let ((q1h (car q1))
(q1l (cdr q1))
(q2h (car q2))
(q2l (cdr q2)))
(+B q1h q2h
(lambda (s e)
(+B s (flo:+ (flo:+ q1l q2l) e)
cons)))))
(define (quad:- q1 q2)
(let ((q1h (car q1))
(q1l (cdr q1))
(q2h (car q2))
(q2l (cdr q2)))
(-B q1h q2h
(lambda (s e)
(+B s (flo:+ (flo:- q1l q2l) e)
cons)))))
(define (quad:* q1 q2)
(let ((q1h (car q1))
(q1l (cdr q1))
(q2h (car q2))
(q2l (cdr q2)))
(*T q1h q2h
(lambda (p e)
(+B p
(flo:+ (flo:+ (flo:* q1l q2h) (flo:* q1h q2l)) e)
cons)))))
(define (quad:/ q1 q2)
(let ((q1h (car q1))
(q1l (cdr q1))
(q2h (car q2))
(q2l (cdr q2)))
(/T q1h q2h
(lambda (q e)
(+B q
(flo:+ (flo:+ (flo:/ q1l q2h)
(flo:* -1.0 (flo:* q (/ q2l q2h))))
e)
cons)))))
(define (->quad r)
(cons (exact->inexact r) 0.0))
(define quad:0 (->quad 0.0))
(define quad:1 (->quad 1.0))
(define quad:2 (->quad 2.0))
(define quad:3 (->quad 3.0))
(define quad:4 (->quad 4.0))
(define *quad-epsilon*
(flo:* *machine-epsilon*
*machine-epsilon*))
(define (quad:square x)
(quad:* x x))
(define (quad:sqrt x)
(let ((xh (car x))
(xl (cdr x)))
(let lp ((guess (->quad (sqrt xh))))
(let ((new (quad:/ (quad:+ guess (quad:/ x guess)) quad:2)))
(if (= (car guess) (car new))
new
(lp new))))))
(define quad:sin-bugger-factor
(expt 2 (* -1 IEEE-p)))
(define (quad:sin x)
(let ((xh (quad:high x)))
(if (< (flo:abs xh) quad:sin-bugger-factor)
x
(let ((y (quad:sin (quad:/ x quad:3))))
(quad:- (quad:* quad:3 y)
(quad:* quad:4
(quad:* y
(quad:* y y))))))))
(define (alternate-sine x)
(let ((xsq (quad:* x x)))
(let lp ((n 1) (xp x) (f quad:1) (sum quad:0))
(let ((term (quad:/ xp f)))
(if (< (abs (quad:high term)) *quad-epsilon*)
sum
(lp (+ n 2)
(quad:* (->quad -1.0) (quad:* xsq xp))
(quad:* f (->quad (* (+ n 1) (+ n 2))))
(quad:+ sum term)))))))
(define quad:high car)
(define quad:low cdr)
#| probably wrong
(define rat:pi 314159265358979323846264338327950288419716939937510/100000000000000000000000000000000000000000000000000)
(define (quad:floor x)
(floor (quad:high x)))
(define (quad:truncate x)
(->quad (truncate (quad:high x)))) ;;; what if the error has an integer part?
(define continued-wallp true)
(define (best-fraction-from-quad x maxint)
(let* ((a0 (inexact->exact (floor (quad:high x))))
(r0 (quad:- x (quad:truncate x))))
(if (<= (quad:high r0) (* (abs (quad:high x)) *quad-epsilon*))
a0
(let* ((d (quad:/ quad:1 r0))
(a1 (inexact->exact (truncate (quad:high d)))))
(let lp ((s d) (last-p a0) (last-q 1) (p (+ (* a0 a1) 1)) (q a1))
(if continued-wallp
(write-line (list s last-p last-q p q)))
(if (> q maxint)
(/ last-p last-q)
(let ((r (quad:- s (quad:truncate s))))
(if (<= (quad:high r) (* (quad:high s) *quad-epsilon*))
(/ p q)
(let* ((new-d (quad:/ quad:1 r))
(a (inexact->exact (truncate (quad:high new-d)))))
(lp new-d
p
q
(+ last-p (* a p))
(+ last-q (* a q))))))))))))
|#
(define (cquad:+ z1 z2)
(map quad:+ z1 z2))
(define (cquad:- z1 z2)
(map quad:- z1 z2))
(define (cquad:* z1 z2)
(list (quad:- (quad:* (car z1) (car z2))
(quad:* (cadr z1) (cadr z2)))
(quad:+ (quad:* (car z1) (cadr z2))
(quad:* (cadr z1) (car z2)))))
(define (cquad:/ z1 z2)
(let ((m (quad:+ (quad:* (car z2) (car z2))
(quad:* (cadr z2) (cadr z2)))))
(list (quad:/ (quad:+ (quad:* (car z1) (car z2))
(quad:* (cadr z1) (cadr z2)))
m)
(quad:/ (quad:- (quad:* (cadr z1) (car z2))
(quad:* (car z1) (cadr z2)))
m))))
(define (q->cq q)
(list q quad:0))
(define (c->cq c)
(list (->quad (real-part c))
(->quad (imag-part c))))
(define (quad:magnitude z)
(quad:sqrt (quad:+ (quad:* (car z) (car z)) (quad:* (cadr z) (cadr z)))))
;;; some other trig functions
(define (compute-quad:pi)
(let lp ((guess (->quad pi)))
(write-line guess)
(let ((new-guess (quad:+ guess (quad:sin guess))))
(if (< (abs (quad:high (quad:- new-guess guess))) *quad-epsilon*)
new-guess
(lp new-guess)))))
(define quad:pi '(3.141592653589793 . 1.2246467991473535e-16))
(define quad:pi/2 '(1.5707963267948966 . 6.123233995736767e-17))
(define (quad:cos x)
(quad:sin (quad:+ x quad:pi/2)))
(define (quad:tan x) (quad:/ (quad:sin x) (quad:cos x)))
(define (quad:atan x)
(let ((fp (quad:+ quad:1 (quad:* x x)))
(guess (->quad (atan (quad:high x)))))
(let lp ((guess guess))
(let ((new-guess (quad:- guess (quad:/ (quad:- (quad:tan guess) x) fp))))
(if (< (abs (quad:high (quad:- new-guess guess))) *quad-epsilon*)
new-guess
(lp new-guess))))))
The current implementation of g/partial-derivative
for ::s/structure
instances is the same implementation used by ::function
.
The reason for this, I think (and the reason it works!) Is that all of the machinery that derivative.cljc
develops to allow for structured inputs (Jacobian support) gets to work a SINGLE time and call the structure as a function for each perturbed input in the incoming structure.
The alternative would be to implement partial-derivative
on structures like this:
(defmethod g/partial-derivative [::s/structure v/seqtype]
[structure selectors]
(s/mapr #(g/partial-derivative % selectors) structure))
This would fix two problems I've noticed:
PowerSeries
, or really any structure whose derivative we can take cheaply WITHOUT calling derivative
and using forward-mode AD, we currently always bypass that implementation. a PowerSeries
in particular can take its own derivative easily, and would stay a PowerSeries
. The current implementation actually calls it, dropping it back down to Series
(and not appropriately shifting all of the entries!!)IPersistentVector
that you haven't turned into a structure, the whole current approach will fail because a vector implements IFn
by looking up an entry, no applying its args to contained functions.Compare these two cases for an example:
sicmutils.calculus.derivative> ((D [g/sin g/cos]) 5)
Execution error (IllegalArgumentException) at sicmutils.calculus.derivative/derivative$fn (derivative.cljc:191).
Key must be integer
sicmutils.calculus.derivative> ((D (s/up g/sin g/cos)) 5)
(up 0.28366218546322625 0.9589242746631385)
It would be really nice if it were possible to build our matrix operations on top of Clojure's core.matrix protocols.
This may be too much! But a big advantage here is that we could plug in to very fast matrix math for our state evolution.
The other possible matrix target to support is https://github.com/uncomplicate/neanderthal. https://github.com/whilo/denisovan exists, if someone wants to use neanderthal
as a core.matrix
implementation... and of course they can use neanderthal directly here. What's important is that the matrix representation be extensible.
This symbolic manipulation library claims core.matrix support: https://github.com/clojure-numerics/expresso and might be worth looking at.
I mentioned goog.math.Matrix
in sicmutils/placeholder#82 , but I'm not sure this is terribly fast. gl-matrix
is very quick, on the other hand. It would be nice to let users do the FAST stuff they want to, and then interact with other sicmutils data types via the generics. If only as a template for how folks can fuse into the ecosystem.
Thanks @qobi for the reference.
This is the solution to the problem of an exponential number of derivatives being calculated for (((expt D n) f) 'x)
== the size of the power set of a set of size n
(see the comment below).
The task is to implement the method described in this paper: https://engineering.purdue.edu/~qobi/papers/popl2007a.pdf
A method is presented for computing all higher-order partial derivatives of a multivariate function R^n → R. This method works by evaluating the function under a nonstandard interpretation, lifting reals to multivariate power series. Multivariate power series, with potentially an infinite number of terms with nonzero coefficients, are represented using a lazy data structure constructed out of linear terms. A complete implementation of this method in SCHEME is presented, along with a straightforward exposition, based on Taylor expansions, of the method’s correctness.
And THEN to figure out how to expose this behavior inside the D
operator as it exists currently.
These three are obviously missing generics; the implementations are not necessarily obvious!
The goal here is to get the bessel functions J0, Y0, J1 and Y1 into Clojure. Once we do this, we can use the recurrence relations in bessel.scm
to implement the binary functions.
I found a nice paper by John Harrison called "Fast and Accurate Bessel Function Computation", with an associated presentation here.
He used a power series approach for part of it, so I spent some time this morning getting that working using our power series library. Too good to pass up, right?
This task is not done, at all, but this is a fast, promising approach that someone should feel free to take on as a challenge. Happy to help stare at the code!
Here's how far I got (also stashed at https://github.com/sicmutils/sicmutils/tree/sritchie/bessel_exploration):
(ns sicmutils.numerical.bessel
(:require [sicmutils.generic :as g]
[sicmutils.util :as u]
[sicmutils.util.stream :as us]
[sicmutils.series :as s]
[sicmutils.value :as v]))
(defn coef
"This inefficiently implements `(n, m)` from section 4 (Asymptotic Expansions)
of the paper, just to check that we have the right idea."
[n m]
(let [num (reduce (fn [acc i]
(* acc (- (* 4 (* n n))
(g/square
(dec (* 2 (inc i)))))))
1
(range 0 m))
den (* (g/expt 2 (* 2 m))
(g/factorial m))]
(/ num den)))
(comment
;; next, these are the example expansions for $P_0(x)$ and $Q_0(x)$, just to
;; verify:
(let [p0-series (-> (s/power-series*
(map (fn [m]
(* (g/expt (- 1) m)
(/ (coef 0 (* 2 m))
(g/expt 2 (* 2 m)))))
(range)))
(s/inflate 2))
q0-series (-> (s/power-series*
(map (fn [m]
(* (g/expt (- 1) m)
(/ (coef 0 (inc (* 2 m)))
(g/expt 2 (inc (* 2 m))))))
(range)))
(s/inflate 2))]
(= '(1
0
(/ -9 128)
0
(/ 3675 32768)
0
(/ -2401245 4194304)
0
(/ 13043905875 2147483648))
(v/freeze (take 9 p0-series)))
(= '((/ -1 8)
0
(/ 75 1024)
0
(/ -59535 262144)
0
(/ 57972915 33554432)
0)
(v/freeze (take 8 q0-series)))))
(defn updater
"Returns a function that generates $(n, m+2)$ given a triple of the numerator,
denominator and `m` from the previous evaluation.
(Returns a triple of the same form.)"
[n]
(let [four*n**2 (* 4 n n)]
(fn [[num denom m]]
(let [two-m (* 2 m)
num' (* num
(- four*n**2 (g/square (+ two-m 1)))
(- four*n**2 (g/square (+ two-m 3))))
den' (* denom
16
(+ m 1)
(+ m 2))]
[num' den' (+ 2 m)]))))
(defn nm-seq
"Given a value of `n` and an optional initial `m` (defaults to 0, only 0 or 1
accepted), returns an infinite sequence of:
[(n, m), (n, m+2), (n, m+4), ...]
Using the definition of (n, m) from the paper."
([n]
(raw-coef** n 0))
([n m]
{:pre [(#{0 1} m)]}
(let [num (if (zero? m)
1
(- (* 4 n n) (g/square (- (* 2 m) 1))))
den (g/expt 2 (* 2 m))
m (u/bigint m)]
(iterate (updater n) [num den m]))))
(defn p-series
"Returns the power series $P_n(x)$ from section 4 of the paper."
[n]
(-> (s/power-series*
(map-indexed
(fn [i [num den two-m]]
(* (g/expt (- 1) i)
(/ num (* den (g/expt 2 two-m)))))
(nm-seq n 0)))
(s/inflate 2)))
(defn q-series
"Returns the power series $Q_n(x)$ from section 4 of the paper.
TODO this uses a goofy `cons 0` to do an efficient multiplication by `x`...
bake this into the power series library."
[n]
(let [unshifted (-> (s/power-series*
(map-indexed
(fn [i [num den two-m+1]]
(* (g/expt (- 1) i)
(/ num (* den (g/expt 2 two-m+1)))))
(raw-coef** n 1)))
(s/inflate 2))]
(s/power-series*
(cons 0 unshifted))))
(defn alpha-series
"Returns the power series $\\alpha(n)$ defined in the 'modified expansions'
section of the paper."
[n]
(s/compose s/atan-series
(g// (g/- (q-series n))
(p-series n))))
(defn beta-series
"Returns the power series $\\beta(n)$ defined in the 'modified expansions'
section of the paper."
[n]
(g/sqrt (g/+ (g/square (p-series n))
(g/square (q-series n)))))
(comment
;; Once again, validate that our efficient version calculates the correct
;; series from the paper.
(= '(1
0
(/ -9 128)
0
(/ 3675 32768)
0
(/ -2401245 4194304)
0
(/ 13043905875 2147483648))
(v/freeze (take 9 (p-series 0))))
(= '((/ -1 8)
0
(/ 75 1024)
0
(/ -59535 262144)
0
(/ 57972915 33554432)
0)
(v/freeze
(take 8 (q-series 0)))))
;; ## Locating the Zeros
;;
;; This doesn't QUITE read like the paper because those power series are in
;; $\frac{1}{x}$, so we need to be a bit careful when trying to multiply both
;; sides by $x$. `identity` here stands for $\frac{1}{x}$.
(def x-factor
(g/+ 1 (g/* s/identity (g/- (alpha-series 0)))))
(def p-factor
(g/invert x-factor))
(def one-over-p
(g/* s/identity p-factor))
(def one-over-x
(s/revert one-over-p))
(def p-over-x
(g// one-over-x s/identity))
(def x-over-p
(g/invert p-over-x))
(defn x-as-fn-of-inv-p
"This is a little wonky; because the series is in 1/p, we have to evaluate it
first then add the final p term back on."
[p]
(g/+ (s/series 0 p)
((g// (g/- x-over-p 1) s/identity)
(g/invert p))))
(comment
(= (0 (/ (+ (* 8N (expt p 2)) 1) (* 8N p))
0
(/ -31N (* 384N (expt p 3)))
0
(/ 3779N (* 15360N (expt p 5)))
0
(/ -6277237N (* 3440640N (expt p 7)))
0
(/ 2092163573N (* 82575360N (expt p 9))))
((g/simplify (take 10 (x-as-fn-of-inv-p 'p))))))
This function checks if the access-chain terminates in another structure or a final symbolic thing. Stare at this more to understand what it is doing... the args
list is a bit confusing. But check the call site to see.
This particular version assumes that we are dealing with a frozen structure.
This is used ONLY in rules.cljc
. Look there in the partial derivative simplifications for where check occurs.
(define (symb:elementary-access? access-chain args)
(define (sea chain thing)
(cond ((and (not (null? chain))
(pair? thing)
(or (eq? (car thing) 'up)
(eq? (car thing) 'down)))
(sea (cdr chain)
(list-ref (cdr thing) (car chain))))
((null? chain)
(or (not (pair? thing))
(not (or (eq? (car thing) 'up)
(eq? (car thing) 'down)))))
(else #f)))
(if (= (length args) 1)
(sea access-chain (car args))
(sea (cdr access-chain)
(list-ref args (car access-chain)))))
We have enough implementations now that this makes sense!
all have sensible arg-shift
and arg-scale
, that keep their types around. Rational functions, series and polynomials can compose without needing to convert to functions too.
NOTE: Please also
These live in the kernel of scmutils. Find these by looking for enable-constructor-simplifications?
and seeing which branch that gives. This is ON by default in scmutils, so we might find some behavior change. Here are some, but please check https://github.com/Tipoca/scmutils/blob/master/src/kernel/numsymb.scm#L34 for more.
;;; The following are more hairy construction simplifications for the
;;; most important functions.
(define (symb1:+ a1 a2)
(cond ((sum? a1)
(cond ((sum? a2)
(addup-args (append (operands a1) (operands a2)) '()))
((difference? a2)
(if (null? (cdr (operands a2)))
(addup-args (operands a1) (operands a2))
(addup-args (append (operands a1)
(list (car (operands a2))))
(cdr (operands a2)))))
(else (addup-args (append (operands a1)
(list a2)) '()))))
((difference? a1)
(if (null? (cdr (operands a1)))
(cond ((sum? a2) (addup-args (operands a2) (operands a1)))
((difference? a2)
(if (null? (cdr (operands a2)))
(addup-args '() (append (operands a1) (operands a2)))
(addup-args (list (car (operands a2)))
(append (operands a1) (cdr (operands a2))))))
(else (addup-args (list a2) (operands a1))))
(cond ((sum? a2)
(addup-args (append (list (car (operands a1))) (operands a2))
(cdr (operands a1))))
((difference? a2)
(if (null? (cdr (operands a2)))
(addup-args (list (car (operands a1)))
(append (cdr (operands a1))
(operands a2)))
(addup-args (list (car (operands a1))
(car (operands a2)))
(append (cdr (operands a1))
(cdr (operands a2))))))
(else (addup-args (list (car (operands a1)) a2)
(cdr (operands a1)))))))
(else
(cond ((sum? a2)
(addup-args (append (list a1) (operands a2)) '()))
((difference? a2)
(if (null? (cdr (operands a2)))
(addup-args (list a1) (operands a2))
(addup-args (append (list a1) (list (car (operands a2))))
(cdr (operands a2)))))
(else (addup-args (list a1 a2) '()))))))
(define (addup-args pos neg)
(define (make-answer sum pos neg)
(if (zero? sum)
(if (null? pos)
(if (null? neg)
:zero
(if (null? (cdr neg))
`(- ,(car neg))
`(- (+ ,@neg))))
(if (null? neg)
(if (null? (cdr pos))
(car pos)
`(+ ,@pos))
(if (null? (cdr pos))
(if (null? (cdr neg))
(if (equal? (car pos) (car neg))
:zero
`(- ,(car pos) ,(car neg)))
`(- ,(car pos) (+ ,@neg)))
(if (null? (cdr neg))
`(- (+ ,@pos) ,(car neg))
(if (equal? pos neg)
:zero
`(- (+ ,@pos) (+ ,@neg)))))))
(if (null? pos)
(if (null? neg)
sum
(if (null? (cdr neg))
`(- ,sum ,(car neg))
`(- ,sum (+ ,@neg))))
(if (null? neg)
`(+ ,sum ,@pos)
(if (null? (cdr neg))
`(- (+ ,sum ,@pos) ,(car neg))
`(- (+ ,sum ,@pos) (+ ,@neg)))))))
(let plp ((p pos) (sum :zero) (respos '()))
(cond ((null? p)
(let nlp ((n neg) (sum sum) (resneg '()))
(cond ((null? n)
(make-answer sum
(reverse respos)
(reverse resneg)))
((number? (car n))
(nlp (cdr n) (- sum (car n)) resneg))
(else
(nlp (cdr n) sum (cons (car n) resneg))))))
((number? (car p))
(plp (cdr p) (+ sum (car p)) respos))
(else
(plp (cdr p) sum (cons (car p) respos))))))
(define (symb1:* a1 a2)
(cond ((product? a1)
(cond ((product? a2)
(mulup-args (append (operands a1) (operands a2)) '()))
((quotient? a2)
(if (null? (cdr (operands a2)))
(mulup-args (operands a1) (operands a2))
(mulup-args (append (operands a1)
(list (car (operands a2))))
(cdr (operands a2)))))
(else (mulup-args (append (operands a1) (list a2)) '()))))
((quotient? a1)
(if (null? (cdr (operands a1)))
(cond ((product? a2) (mulup-args (operands a2) (operands a1)))
((quotient? a2)
(if (null? (cdr (operands a2)))
(mulup-args '() (append (operands a1) (operands a2)))
(mulup-args (list (car (operands a2)))
(append (operands a1) (cdr (operands a2))))))
(else (mulup-args (list a2) (operands a1))))
(cond ((product? a2)
(mulup-args (append (list (car (operands a1))) (operands a2))
(cdr (operands a1))))
((quotient? a2)
(if (null? (cdr (operands a2)))
(mulup-args (list (car (operands a1)))
(append (cdr (operands a1))
(operands a2)))
(mulup-args (list (car (operands a1))
(car (operands a2)))
(append (cdr (operands a1))
(cdr (operands a2))))))
(else (mulup-args (list (car (operands a1)) a2)
(cdr (operands a1)))))))
(else
(cond ((product? a2)
(mulup-args (append (list a1) (operands a2)) '()))
((quotient? a2)
(if (null? (cdr (operands a2)))
(mulup-args (list a1) (operands a2))
(mulup-args (append (list a1) (list (car (operands a2))))
(cdr (operands a2)))))
(else (mulup-args (list a1 a2) '()))))))
(define (mulup-args pos neg)
(define (make-product numfact factors)
(if (null? factors)
numfact
(if (one? numfact)
(if (null? (cdr factors))
(car factors)
`(* ,@factors))
(if (null? (cdr factors))
`(* ,numfact ,(car factors))
`(* ,numfact ,@factors)))))
(define (make-answer pfactor pos nfactor neg)
(let ((num (make-product pfactor pos))
(den (make-product nfactor neg)))
(cond ((and (number? den) (zero? den))
(error "zero divide in mulup-args"))
((and (number? num) (zero? num))
:zero)
((and (number? den) (one? den))
num)
((and (number? num) (number? den))
(/ num den))
((equal? num den)
:one)
(else
`(/ ,num ,den)))))
(let plp ((p pos) (pfactor :one) (respos '()))
(cond ((null? p)
(let nlp ((n neg) (nfactor :one) (resneg '()))
(cond ((null? n)
(make-answer pfactor
(reverse respos)
nfactor
(reverse resneg)))
((number? (car n))
(nlp (cdr n) (* nfactor (car n)) resneg))
(else
(nlp (cdr n) nfactor (cons (car n) resneg))))))
((number? (car p))
(plp (cdr p) (* pfactor (car p)) respos))
(else
(plp (cdr p) pfactor (cons (car p) respos))))))
(define (symb1:- a1 a2)
(cond ((sum? a1)
(cond ((sum? a2)
(addup-args (operands a1) (operands a2)))
((difference? a2)
(if (null? (cdr (operands a2)))
(addup-args (append (operands a1) (operands a2)) '())
(addup-args (append (operands a1) (cdr (operands a2)))
(list (car (operands a2))))))
(else (addup-args (operands a1) (list a2)))))
((difference? a1)
(if (null? (cdr (operands a1)))
(cond ((sum? a2)
(addup-args '() (append (operands a1) (operands a2))))
((difference? a2)
(if (null? (cdr (operands a2)))
(addup-args (operands a2) (operands a1))
(addup-args (cdr (operands a2))
(append (operands a1)
(list (car (operands a2)))))))
(else (addup-args '() (append (operands a1) (list a2)))))
(cond ((sum? a2)
(addup-args (list (car (operands a1)))
(append (cdr (operands a1)) (operands a2))))
((difference? a2)
(if (null? (cdr (operands a2)))
(addup-args (append (list (car (operands a1)))
(operands a2))
(cdr (operands a1)))
(addup-args (cons (car (operands a1))
(cdr (operands a2)))
(append (cdr (operands a1))
(list (car (operands a2)))))))
(else (addup-args (list (car (operands a1)))
(append (cdr (operands a1)) (list a2)))))))
(else
(cond ((sum? a2)
(addup-args (list a1) (operands a2)))
((difference? a2)
(if (null? (cdr (operands a2)))
(addup-args (append (list a1) (operands a2)) '())
(addup-args (append (list a1) (cdr (operands a2)))
(list (car (operands a2))))))
(else (addup-args (list a1) (list a2)))))))
(define (symb1:/ a1 a2)
(cond ((product? a1)
(cond ((product? a2)
(mulup-args (operands a1) (operands a2)))
((quotient? a2)
(if (null? (cdr (operands a2)))
(mulup-args (append (operands a1) (operands a2)) '())
(mulup-args (append (operands a1) (cdr (operands a2)))
(list (car (operands a2))))))
(else (mulup-args (operands a1) (list a2)))))
((quotient? a1)
(if (null? (cdr (operands a1)))
(cond ((product? a2)
(mulup-args '() (append (operands a1) (operands a2))))
((quotient? a2)
(if (null? (cdr (operands a2)))
(mulup-args (operands a2) (operands a1))
(mulup-args (cdr (operands a2))
(append (operands a1)
(list (car (operands a2)))))))
(else (mulup-args '() (append (operands a1) (list a2)))))
(cond ((product? a2)
(mulup-args (list (car (operands a1)))
(append (cdr (operands a1)) (operands a2))))
((quotient? a2)
(if (null? (cdr (operands a2)))
(mulup-args (append (list (car (operands a1)))
(operands a2))
(cdr (operands a1)))
(mulup-args (cons (car (operands a1))
(cdr (operands a2)))
(append (cdr (operands a1))
(list (car (operands a2)))))))
(else (mulup-args (list (car (operands a1)))
(append (cdr (operands a1)) (list a2)))))))
(else
(cond ((product? a2)
(mulup-args (list a1) (operands a2)))
((quotient? a2)
(if (null? (cdr (operands a2)))
(mulup-args (append (list a1) (operands a2)) '())
(mulup-args (append (list a1) (cdr (operands a2)))
(list (car (operands a2))))))
(else (mulup-args (list a1) (list a2)))))))
We need to:
(sicmutils.structure/orientation->separator :sicmutils.structure/up)
to some easier variable and expose it(raise-sym 'x 3)
that will generate x↑3
So users don't have to type that in!
@alexgian sent me a lovely paper by the legendary Bill Gosper paper on the continued fraction representation of real numbers.
This would be an AWESOME data type to add to SICMUtils, and a great target for a literate-style essay, like the numerical methods work or Power Series.
From the site:
"Multiprecision arithmetic algorithms usually represent real numbers as decimals, or perhaps as their base-2n analogues. But this representation has some puzzling properties. For example, there is no exact representation of even as simple a number as one-third. Continued fractions are a practical but little-known alternative.
"Continued fractions are a representation of the real numbers that are in many ways more mathematically natural than the usual decimal or binary representations. All rational numbers have simple representations, and so do many irrational numbers, such as sqrt(2) and e1. One reason that continued fractions are not often used, however, is that it's not clear how to involve them in basic operations like addition and multiplication. This was an unsolved problem until 1972, when Bill Gosper found practical algorithms for continued fraction arithmetic."
These come packaged with Google Closure; it would be great to have these available as another matrix library in SICMUtils. If we don't depend on the namespace by default, then anyone compiling Clojurescript won't see it show up in their output.
@littleredcomputer , wdyt... is this the sort of thing that would be better served in a new, small project? sicmutils-closure
, to pick a very confusing name choice?
EDIT I think emmy.goog.math
is a good namespace here.
I saw a reference to "Arb: Efficient Arbitrary-Precision Midpoint-Radius Interval Arithmetic" on hackernews; it would be awesome to get proper interval arithmetic in to sicmutils.
More references:
Imagine an instantly-firing-up terminal calculator. We can do this GraalVM. Some resources for thinking this through:
Here's a nice paper we can add: https://blog.acolyer.org/2020/10/02/toward-an-api-for-the-real-numbers/
I THINK this is implemented in Spire: https://github.com/typelevel/spire/blob/master/core/src/main/scala/spire/math/Real.scala
I have my eyes on https://github.com/JuliaApproximation/FastGaussQuadrature.jl as a nice target to clone here. This would be an awesome addition to our evenly-spaced-points methods we have now!
It would be lovely if Nextjournal and friends had some sicmutils-version
variable inside env
that they could query to figure out what is actually running! I'm not sure of the best pattern here. Does anyone have any pointers?
Here is how to do it on the JVM, modeled after *clojure.core/clojure-version*
:
(let [version-string (slurp (.getResourceAsStream
(clojure.lang.RT/baseLoader)
"SICMUTILS_VERSION"))
[_ major minor incremental qualifier snapshot]
(re-matches
#"(\d+)\.(\d+)\.(\d+)(?:-([a-zA-Z0-9_]+))?(?:-(SNAPSHOT))?"
version-string)
version {:major (Integer/valueOf ^String major)
:minor (Integer/valueOf ^String minor)
:incremental (Integer/valueOf ^String incremental)
:qualifier (if (= qualifier "SNAPSHOT") nil qualifier)}]
(def ^:dynamic *sicmutils-version*
(if (.contains version-string "SNAPSHOT")
(assoc version :interim true)
version)))
Right now the compile
functions use a dictionary called compiled-fn-whitelist to associate symbols with functions. After simplification and common subexpression elimination, the compiler swaps these functions back in for their associated symbols and calls eval
on the result.
There are a few alternatives a user might want to the current compiled-fn-whitelist:
Geex looks like it might be the easiest to integrate.
The code for this lives in the simplify package, in
Candidates:
UPDATE: I started working on this a bit as part of the monster sicmutils/sicmutils#341 effort, then backed off. Here's the code I had started sketching to make this work:
:else
clause in GCD dispatch...
(def ^:dynamic *euclid-breakpoint-arity* 3)
(def ^:dynamic *gcd-cut-losses* nil)
(or (with-limited-time [1.0 :seconds]
(fn [] (sparse-gcd u v)))
(with-limited-time [1.0 :seconds]
(fn [] (classical-gcd u v)))
(with-limited-time [100.0 :seconds]
(fn [] (sparse-gcd u v)))
(and (< arity *euclid-breakpoint-arity*)
(with-limited-time [100.0 :seconds]
(fn [] (classical-gcd u v))))
(sparse-gcd u v)
(if *gcd-cut-losses*
(or (with-limited-time *gcd-cut-losses*
(fn [] (classical-gcd u v)))
1)
(classical-gcd u v)))
The task here is to get a writeup of what we want interactive math and physics explorations, on their own and in the context of a textbook, to feel like with a Clojure library like sicmutils. We can start this in a discussion here, but the goal is to populate:
add trig functions for js/BigInt
, goog.math.{Long, Integer}
.
These work as of sicmutils/sicmutils#168 , but they work by converting these types to js/Number
, which won't work if they're too large.
This is fine if we can first confine the numbers back down to a reasonable range.
See Roylance's "Expressing Mathematical Routines Constructively" for a scheme implementation that does this.
to match the CLJ behavior. Here's an example of how clj does this:
TODO gotta get promoting arithmetic in cljs! here is how https://github.com/clojure/clojure/blob/master/src/jvm/clojure/lang/Numbers.java#L488 file a ticket.
This should be user-controllable, at the expression level and the system level.
You really just want simplification the FIRST time. Then if you call "factor" or something, you want that to not run through the simplifier again and UN-factor. So explicit "simplify" style calls should not trigger further simplification. How to make that happen?
Pros:
Cons:
I'm making a ticket about this because during my deep dive on Brent's Method for minimization, I discovered that the algorithm described on the wiki page is in fact wrong! https://en.wikipedia.org/wiki/Brent%27s_method
So the task here is to:
Right now, the Complex
datatype is based on
js/Number
in Clojurescript due to the Complex.js dependencyjosdejong/mathjs#694 has some discussion about the problems this has created in math.js. The big one is that you can't extend complex arithmetic to arbitrary precision.
We COULD probably do this if we rebuilt Complex in Clojurescript, on top of the generic arithmetic. (I think you might be able to get complex matrices for free??)
https://github.com/infusion/Complex.js/blob/master/complex.js is not that big of a library. porting it over to cljc
, and backing it with generic operations, would be a really nice project.
This MathOverflow comment mentions a nice way to do this:
Also, Sussman in "Abstraction in Numerical Methods" (bottom of page 6) notes:
By the way, if the sequence
$p_1, p_2, p_3, ...$ is not known in advance... Alternatively, the appropriate exponent$p$ for a given column can be inferred numerically from the early terms of that column; we have done that in our library, but do not pursue it here.
I can't find the actual library he's mentioning, so we'll have to figure this out ourselves!
@kovasb has some excellent experiments in https://github.com/kovasb/combinator on how to do very fast term rewriting in clj/cljs, 2-3x faster than Mathematica is able to support. @littleredcomputer , you may have studied this already (please close the issue if so!) but I'd love your thoughts on the zipper-based implementation here, and whether it might give us some advantages over the current implementation.
Some other Clojure term-rewriting systems:
Nice notes describing this and related papers:
/**
* Algebraic provides an exact number type for algebraic numbers. Algebraic
* numbers are roots of polynomials with rational coefficients. With it, we can
* represent expressions involving addition, multiplication, division, n-roots
* (eg. `sqrt` or `cbrt`), and roots of rational polynomials. So, it is similar
* [[Rational]], but adds roots as a valid, exact operation. The cost is that
* this will not be as fast as [[Rational]] for many operations.
*
* In general, you can assume all operations on this number type are exact,
* except for those that explicitly construct approximations to an Algebraic
* number, such as `toBigDecimal`.
*
* For an overview of the ideas, algorithms, and proofs of this number type,
* you can read the following papers:
*
* - "On Guaranteed Accuracy Computation." C. K. Yap.
* - "Recent Progress in Exact Geometric Computation." C. Li, S. Pion, and C. K. Yap.
* - "A New Constructive Root Bound for Algebraic Expressions" by C. Li & C. K. Yap.
* - "A Separation Bound for Real Algebraic Expressions." C. Burnikel, et al.
*/
Currently this code just works on functions. If you memoize something else, it will get wrapped in a function.
Make memoize generic, and install it for:
fmap
MetaFn
and Fn
etcThen do some perf testing to see if this slows things down! Maybe we have a function-specific memoize
and then a dispatching one too, so if you KNOW you are using a function you can specialize.
rcf.scm has a compose
function that converts to, roughly, this:
;; The following function still reads slightly mystical; this comes from rcf.scm
;; in the scmutils library. I'm sure with some staring and understanding we
;; could build a version that allows composition with any variable in the
;; left [[RationalFunction]] instance.
(defn compose
"only plugs r2 in for the principal indeterminate.
TODO we COULD do a version that composes with a different one??"
[r1 r2]
{:pre [(rational-function? r1)]}
(if-not (rational-function? r2)
(g/div (p/evaluate (bare-u r1) [r2])
(p/evaluate (bare-v r1) [r2]))
(let [nr1 (bare-u r1)
nr2 (bare-u r2)
dr1 (bare-v r1)
dr2 (bare-v r2)
dn (p/degree nr1)
dd (p/degree dr1)
narity (+ (p/arity dr1) 1)
nnr1 (p/extend 1 (p/reciprocal nr1))
ndr1 (p/extend 1 (p/reciprocal dr1))
scales [(second (p/new-variables narity)) 1]
pn (p/evaluate (p/reciprocal
(p/arg-scale nnr1 scales))
[nr2 dr2])
pd (p/evaluate (p/reciprocal
(p/arg-scale ndr1 scales))
[nr2 dr2])]
(cond (> dn dd) (g/div pn (p/mul (p/expt dr2 (- dn dd)) pd))
(< dn dd) (g/div (p/mul (p/expt dr2 (- dd dn)) pn) pd)
:else (g/div pn pd)))))
But it's a bit mystical and I don't see how it works yet. The task is to get this working, and get a similar version working for polynomials.
Hey!
disclaimer: There are many moving pieces of sicmutils I don't understand; I find this to be quite difficult.
NOTE: this paper, which discusses "dimensional algebra": https://fermatslibrary.com/s/can-one-take-the-logarithm-or-the-sine-of-a-dimensioned-quantity-or-a-unit-dimensional-analysis-involving-transcendental-functions#email-newslettter
Also note the scmutils files that implement this live in scmutils/units:
I'd love to see first class unit support. I asked in the Clojurians Slack #sicmutils channel, and got some pointers from @sritchie. Messages pasted for context:
Have you given any thought to putting SI units on numbers; allowing 1 km / 30 minutes = 0.5 km/h and similar calculations? You talked a bit about the extensibility of the number systems (multimethods for +, *, ...) on the scicloj meetup.
Or, one could define plain SI units symbolically, and derived SI units in terms of base SI units. Perhaps adding defmethods and a new type is the wrong way to go.
it’s a brilliant idea! In fact the original scmutils has a “units” package that I would love to mentor / help someone port over into sicmutils
it takes a while to learn how scmutils works so it would be a tough go on your own, but if you’re interested this would be an AWESOME project
https://github.com/Tipoca/scmutils/blob/master/src/units/with-units.scm#L243
that would let us put all these constants in! https://github.com/Tipoca/scmutils/blob/master/src/units/constants.scm#L108
Another note: the https://github.com/Tipoca/scmutils/tree/master/src/units[`units`] folder hosts a bunch of physical constants, plus arithmetic on units and the ability to work with numbers and expressions that have units attached.
This is how to do this:
sicmutils.env> (freeze (with-meta (D (literal-function 'f)) {:name "cake"}))
"cake"
sicmutils.env> (with-meta (D (literal-function 'f)) {:name "cake"})
#function[clojure.lang.AFunction/1]
@littleredcomputer , curious for your thoughts here.
As I work through the cljs port sicmutils/sicmutils#40 , clear layers are emerging:
Complex
and friends, that rely on dependencies.Matrix
, Series
, etcFor some reason it's much more natural in other build systems like SBT to offer modules: https://github.com/twitter/algebird so maybe this is not a good idea here... but it does feel like it might be a nice move to explicitly break the library into these modules, so that someone could depend on JUST the numerical stuff, or JUST the automatic differentiation, without pulling in the rest of the dependencies for the book.
If we start adding better tools for visualization into the library, like a Quil dependency, I can see this becoming even more true.
EDIT: Here's an example, from the Clojurians slack, of how this might work: https://github.com/juxt/crux/blob/master/project.clj
This is how the operation is defined in Scheme:
(define (g:gcd:n args)
(cond ((null? args) :zero)
((null? (cdr args)) (car args))
(else
(let lp
((as (cddr args))
(ans (g:gcd:bin (car args) (cadr args))))
(cond ((null? as) ans)
((g:one? ans) ans)
(else
(lp (cdr as) (g:gcd:bin ans (car as)))))))))
The task here would be to
g/-gcd
and g/-lcm
g/gcd
, g/lcm
that wraps them and does this trick with both functions.I've been thinking about what other renderers might make sense, in addition to "infix" and "TeX"... a couple of thoughts (cc @mk).
Here is some prior art for the latter. Here's an example of sympy doing this trick.
We need these to get proper comparisons working between numbers and literal-number
and Differential
instances. We have sicmutils.value/compare
, so these should be simple to add and document.
I think they should replace the clojure.core versions in sicmutils.env.
Also: This takes us into territory well covered by clojure.algo.generic
https://github.com/clojure/algo.generic/blob/master/src/main/clojure/clojure/algo/generic/comparison.clj. We should in fact add the full suite of operations, following their lead.
Right now these print as the thing that's wrapped, leading to confusing error messages like
0 does not equal 0
in tests. If we make these print as #sicm/expression
, it will all be much clearer. Here's what it will take:
;;print-method
(let [rep (str "#" (.-sym (.-type s))
" " (.toString s))]
(.write w ^String rep))
;; clojurescript version
(let [rep (str "#" (.-sym (type s))
" " (pr-str expression))]
(-write writer rep))
Plus an entry in data_readers.clj
.
NOTE that we also want print-dup
methods for some of the types!
This would let us do ((partial :key) f) {:key 10 :val 100}
. I am thinking more about supporting maps on the input...
Of course this ticket also implies that D
should return a Jacobian shaped like a map. That breaks the contract that the output is compatible for contraction with an increment in the input.
Do we restrict this behavior to Grad
, which keeps orientations preserved?
This is an idea from scmutils that is not actually used there. The idea is to add metadata to Literal
instances that states more specific claims we know about those expressions.
Actually, while typing this I see that we can potentially add assumptions like positive?
etc so forms as they're simplified. So this is a bigger idea.
But for now, consider metadata that would make a literal respond to real?
, integer?
etc, perhaps with a different v/kind
.
Here are notes from a previous think on this:
known-real?
is a predicate in numbers.scm
- this is a system for registering symbols as "real". The symbolic environment (from numsymb.scm
) uses this information in the call to g/conjugate
for symbols and just mirrors the term back out. We could of course do this in the other complex calls.
We could ALSO take this unused idea from express.scm
:
(define (make-real-literal expression)
(let ((e (make-numerical-literal expression)))
(add-property! e 'real #t)
e))
and use the metadata field that's now in Literal
. BUT thinking now, we can actually just attach metadata to symbols themselves and have them respond directly also to known-real?
.
UPDATE: symbolic expressions and symbols now both respond to with-meta
. All we need here is:
known-real?
function:real? true
and metadata like this.(define *known-reals* '())
(define (known-real? z)
(cond ((structure? z)
(s:forall known-real? z))
((matrix? z)
(let ((m (m:num-rows matrix))
(n (m:num-cols matrix))
(mat (matrix->array z)))
(let rowlp ((i 0))
(if (fix:= i m)
#t
(let collp ((j 0))
(if (fix:= j n)
(rowlp (fix:+ i 1))
(if (known-real? (array-ref mat i j))
(collp (fix:+ j 1))
#f)))))))
((differential? z)
(for-all? (differential->terms z)
(lambda (term)
(known-real? (differential-coefficient term)))))
(else
(there-exists? *known-reals*
(lambda (w)
(or (equal? w z)
(let ((diff
(ignore-errors
(lambda ()
(simplify (g:- w z))))))
(and (not (condition? diff))
(exact-zero? diff)))))))))
;;; Permanent declaration
(define (declare-known-reals . stuff)
(set! *known-reals* (list-union stuff *known-reals*)))
(define (declare-unknown-reals . stuff)
(set! *known-reals* (list-difference *known-reals* stuff)))
;;; Temporary declaration
(define (with-known-reals stuff thunk)
(fluid-let ((*known-reals*
(list-union stuff *known-reals*)))
(thunk)))
Right now, compile-fn
(link) can handle functions of different arities, and the GUTS of compile-state-fn
can handle arguments with any structure.
It would be great to unify these, so compile-state-fn
called the guts of compile-fn
, and the latter could take an arbitrary argument template instead of just n
for its arity.
Some nice design work waiting here!
@hcarvalhoaves and CTW in Clojurians Slack noted that this would be great.
Some tasks:
=
and inequality expressions::v/numeric
->infix
etcWe might want inequalities represented with a new record type, rather than an expression with an embedded <=
. But they should certainly freeze to that.
This is hard because we'll need to detect differentials in the input and respond there, if we want to do this. It all works fine IF you restrict everything to the numerator...
sd/IPerturbed
(perturbed? [_]
(sd/perturbed? u))
(replace-tag [this old new]
(RationalFunction. arity
(sd/replace-tag u old new)
v
m))
(extract-tangent [this tag]
(RationalFunction. arity
(sd/extract-tangent u tag)
v
m))
But this is not very genuine.
Here's a test that shows how this can get wacky:
(testing "IPerturbed"
(letfn [(f [x]
(rf/make
(p/make [1 2 (g/square x) 3])
(p/make [5 (g/square x)])))
(g [x y]
(let [n (g/+ 1 (g/* 2 y)
(g/* (g/square x) (g/square y))
(g/* 3 (g/cube y)))
d (g/+ 5 (g/* (g/square x) y))]
(g// n d)))]
(((D f) 'x) 'y)
(((cd/partial 0) g) 'x 'y)))
This is a nice system to take some derivation and return an unapplied functional version of the result. Here is the code, and some examples from Scheme.
Once we have this we'll need to
present
function in the testssuppress-arguments
with-suppressed-args
version that does the same but only inside of a let binding. Or maybe just call it suppress
...(define *suppressed-argument-list* '())
(define *suppressed-argument-list-counter* 0)
(define (suppress-arguments arguments)
(let ((n (+ (length *suppressed-argument-list*) 1)))
(set! *suppressed-argument-list-counter*
(+ *suppressed-argument-list-counter* 1))
(set! *suppressed-argument-list*
(cons (cons arguments
(symbol "args."
*suppressed-argument-list-counter*))
*suppressed-argument-list*))
n))
(define (show-suppressed-arguments)
(pp (map (lambda (al)
`(,(cdr al) = ,@(car al)))
*suppressed-argument-list*)))
(define (clear-arguments)
(set! *suppressed-argument-list* '())
(set! *suppressed-argument-list-counter* 0)
0)
(define (arg-suppressor expression)
(if (pair? expression)
(let ((v (assoc (cdr expression) *suppressed-argument-list*)))
(if v
(list (arg-suppressor (car expression)) (cdr v))
(cons (arg-suppressor (car expression))
(arg-suppressor (cdr expression)))))
expression))
;;; For example
(let ((t 't) (xy (up 'x 'y)) (uv (up 'r 's)))
(* (((partial 2) Hp) (up t uv (- (((partial 2) F1) t xy uv))))
(((partial 2) ((partial 1) F1)) 't xy uv)))
#|
(down
(+
(*
(((partial 1 0) ((partial 2 1) F1)) t (up x y) (up r s))
(((partial 2 1) Hp)
(up t
(up r s)
(down (* -1 (((partial 2 0) F1) t (up x y) (up r s)))
(* -1 (((partial 2 1) F1) t (up x y) (up r s)))))))
...mess...)
...mess...)
|#
;;; We choose arguments to suppress:
(suppress-arguments '((up t
(up r s)
(down (* -1 (((partial 2 0) F1) t (up x y) (up r s)))
(* -1 (((partial 2 1) F1) t (up x y) (up r s)))))))
#| 1 |#
(suppress-arguments '(t (up x y) (up r s)))
#| 2 |#
;;; Now look at the pretty result:
(let ((t 't) (xy (up 'x 'y)) (uv (up 'r 's)))
(* (((partial 2) Hp) (up t uv (- (((partial 2) F1) t xy uv))))
(((partial 2) ((partial 1) F1)) 't xy uv)))
#|
(down
(+ (* (((partial 2 0) Hp) args.1) (((partial 1 0) ((partial 2 0) F1)) args.2))
(* (((partial 2 1) Hp) args.1) (((partial 1 0) ((partial 2 1) F1)) args.2)))
(+ (* (((partial 2 0) Hp) args.1) (((partial 1 1) ((partial 2 0) F1)) args.2))
(* (((partial 2 1) Hp) args.1) (((partial 1 1) ((partial 2 1) F1)) args.2))))
|#
(show-suppressed-arguments)
((args.2 = t (up x y) (up r s))
(args.1 = (up t
(up r s)
(down (* -1 (((partial 2 0) F1) t (up x y) (up r s)))
(* -1 (((partial 2 1) F1) t (up x y) (up r s)))))))
scmutils supports the following ODE solvers in its numerics package:
qcceuler
?)We currently only have support for Bulirsch-Stoer (the best, default!), via Colin Smith's "ODEX" library in Javascript and the Apache-Commons implementation on the JVM.
Fleshing this out with the care we used for the quadrature package would be brilliant.
here are the files:
We do support such literal functions, but with an undocumented syntax for domain and range. This should be harmonized with the book.
@barak and @qobi presented on a series of AD benchmarks that they performed as part of their work on Stalingrad:
scmutils didn't do terribly well for a few reasons, one of which is no compilation of the derivatives... but no matter!
The task here is to get the SICMUtils code for these benchmarks tidied up, measured and going faster.
I took a pass at getting the code mechanically ported from the two straightforward benchmarks that worked in scmutils, particle-FF
and saddle-FF
. These are described in the first paper linked above.
NOTE: probabilistic-lambda-calculus
, probabilistic-prolog
and backprop
all have forward-mode examples too that would work with some porting. Also once I get sicmutils/sicmutils#226 up and running we can get the rest of the benchmarks running here too. I know they'd be faster with a more intelligent compile
implementation that could handle multiple arguments more gracefully. This would need to generate multiple nested functions in forward mode.
;; Run these from sicmutils.env... I have some defensive g/ prefixes that we probably don't need.
(defn multivariate-argmin [f x]
(let [g (D f)]
(loop [x x
fx (f x)
gx (g/transpose (g x))
eta 1e-5
i 0]
(cond (<= (compare (g/abs gx) 1e-5) 0)
x
(= i 10)
(recur x fx gx (g/* 2.0 eta) 0)
:else
(let [x-prime (g/- x (g/* eta gx))]
(if (<= (compare (g/abs (- x x-prime)) 1e-5) 0)
x
(let [fx-prime (f x-prime)]
(if (< (compare fx-prime fx) 0)
(recur x-prime fx-prime (g/transpose (g x-prime)) eta (+ i 1))
(recur x fx gx (/ eta 2.0) 0)))))))))
(defn naive-euler [w]
(let [charges [[10.0 (- 10.0 w)] [10.0 0.0]]
x-initial [0.0 8.0]
xdot-initial [0.75 0.0]
delta-t 1e-1
p (fn [x]
(transduce (map (fn [c] (/ 1.0 (g/abs (- x c)))))
g/+
0.0
charges))]
(loop [x x-initial
xdot xdot-initial]
(let [[x' y' :as x-new] (g/+ x (g/* delta-t xdot))]
(if (g/negative? y')
(let [delta-t-f (g// (g/- 0.0 (second x))
(second xdot))
x-t-f (g/+ x (g/* delta-t-f xdot))]
(g/square (first x-t-f)))
(let [xddot (g/* -1.0 (g/transpose ((D p) x)))]
(recur x-new (g/+ xdot (g/* delta-t xddot)))))))))
(defn multivariate-argmax [f x]
(multivariate-argmin (fn [x] (g/- (f x))) x))
(defn multivariate-max [f x]
(f (multivariate-argmax f x)))
(defn particle-FF []
(multivariate-argmin naive-euler 0.0))
(defn saddle-FF []
(let [start [1.0 1.0]
f (fn [[x1 y1] [x2 y2]]
(- (+ (g/square x1) (g/square y1))
(+ (g/square x2) (g/square y2))))
x1* (multivariate-argmin
(fn [x1]
(multivariate-max
(fn [x2] (f x1 x2))
start))
start)
x2* (multivariate-argmax
(fn [x2] (f x1* x2)) start)]
[x1* x2*]))
And the results (with nothing to compare it to, since this is a new M1 mac and I haven't run any other benchmarks here):
sicmutils.env> (time (saddle-FF))
"Elapsed time: 9141.064291 msecs"
[(up 8.246324826140356E-6 8.246324826140356E-6) (up 8.246324826140356E-6 8.246324826140356E-6)]
sicmutils.env> (time (particle-FF))
"Elapsed time: 28148.853958 msecs"
0.20719187464861197
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.