Mercurial > dylan
view sicm/utils.html @ 4:10c30f787f4b
minor change to quandry
| author | Robert McIntyre <rlm@mit.edu> |
|---|---|
| date | Fri, 28 Oct 2011 00:07:59 -0700 |
| parents | b4de894a1e2e |
| children |
line wrap: on
line source
1 <?xml version="1.0" encoding="utf-8"?>2 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"3 "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">4 <html xmlns="http://www.w3.org/1999/xhtml"5 lang="en" xml:lang="en">6 <head>7 <title>Building a Classical Mechanics Library in Clojure</title>8 <meta http-equiv="Content-Type" content="text/html;charset=utf-8"/>9 <meta name="generator" content="Org-mode"/>10 <meta name="generated" content="2011-08-11 04:10:36 EDT"/>11 <meta name="author" content="Dylan Holmes"/>12 <meta name="description" content=""/>13 <meta name="keywords" content=""/>14 <style type="text/css">15 <!--/*--><![CDATA[/*><!--*/16 html { font-family: Times, serif; font-size: 12pt; }17 .title { text-align: center; }18 .todo { color: red; }19 .done { color: green; }20 .tag { background-color: #add8e6; font-weight:normal }21 .target { }22 .timestamp { color: #bebebe; }23 .timestamp-kwd { color: #5f9ea0; }24 .right {margin-left:auto; margin-right:0px; text-align:right;}25 .left {margin-left:0px; margin-right:auto; text-align:left;}26 .center {margin-left:auto; margin-right:auto; text-align:center;}27 p.verse { margin-left: 3% }28 pre {29 border: 1pt solid #AEBDCC;30 background-color: #F3F5F7;31 padding: 5pt;32 font-family: courier, monospace;33 font-size: 90%;34 overflow:auto;35 }36 table { border-collapse: collapse; }37 td, th { vertical-align: top; }38 th.right { text-align:center; }39 th.left { text-align:center; }40 th.center { text-align:center; }41 td.right { text-align:right; }42 td.left { text-align:left; }43 td.center { text-align:center; }44 dt { font-weight: bold; }45 div.figure { padding: 0.5em; }46 div.figure p { text-align: center; }47 textarea { overflow-x: auto; }48 .linenr { font-size:smaller }49 .code-highlighted {background-color:#ffff00;}50 .org-info-js_info-navigation { border-style:none; }51 #org-info-js_console-label { font-size:10px; font-weight:bold;52 white-space:nowrap; }53 .org-info-js_search-highlight {background-color:#ffff00; color:#000000;54 font-weight:bold; }55 /*]]>*/-->56 </style>57 <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />58 <script type="text/javascript">59 <!--/*--><![CDATA[/*><!--*/60 function CodeHighlightOn(elem, id)61 {62 var target = document.getElementById(id);63 if(null != target) {64 elem.cacheClassElem = elem.className;65 elem.cacheClassTarget = target.className;66 target.className = "code-highlighted";67 elem.className = "code-highlighted";68 }69 }70 function CodeHighlightOff(elem, id)71 {72 var target = document.getElementById(id);73 if(elem.cacheClassElem)74 elem.className = elem.cacheClassElem;75 if(elem.cacheClassTarget)76 target.className = elem.cacheClassTarget;77 }78 /*]]>*///-->79 </script>81 </head>82 <body>84 <div id="content">88 <div class="header">89 <div class="float-right">90 <!--91 <form>92 <input type="text"/><input type="submit" value="search the blog »"/>93 </form>94 -->95 </div>97 <h1>aurellem <em>☉</em></h1>98 <ul class="nav">99 <li><a href="/">read the blog »</a></li>100 <!-- li><a href="#">learn about us »</a></li-->101 </ul>102 </div>104 <h1 class="title">Building a Classical Mechanics Library in Clojure</h1>112 <div id="table-of-contents">113 <h2>Table of Contents</h2>114 <div id="text-table-of-contents">115 <ul>116 <li><a href="#sec-1">1 Useful data types </a>117 <ul>118 <li><a href="#sec-1-1">1.1 Complex numbers </a></li>119 <li><a href="#sec-1-2">1.2 Power series </a></li>120 <li><a href="#sec-1-3">1.3 Tuples and tensors </a>121 <ul>122 <li><a href="#sec-1-3-1">1.3.1 Tuples are “sequences with spin” </a></li>123 <li><a href="#sec-1-3-2">1.3.2 Matrices </a></li>124 </ul>125 </li>126 <li><a href="#sec-1-4">1.4 Generic Operations </a></li>127 </ul>128 </li>129 <li><a href="#sec-2">2 Operators and Differentiation </a>130 <ul>131 <li><a href="#sec-2-1">2.1 Operators </a></li>132 <li><a href="#sec-2-2">2.2 Differential Terms and Sequences </a></li>133 </ul>134 </li>135 </ul>136 </div>137 </div>139 <div id="outline-container-1" class="outline-2">140 <h2 id="sec-1"><span class="section-number-2">1</span> Useful data types </h2>141 <div class="outline-text-2" id="text-1">145 </div>147 <div id="outline-container-1-1" class="outline-3">148 <h3 id="sec-1-1"><span class="section-number-3">1.1</span> Complex numbers </h3>149 <div class="outline-text-3" id="text-1-1">152 </div>154 </div>156 <div id="outline-container-1-2" class="outline-3">157 <h3 id="sec-1-2"><span class="section-number-3">1.2</span> Power series </h3>158 <div class="outline-text-3" id="text-1-2">161 </div>163 </div>165 <div id="outline-container-1-3" class="outline-3">166 <h3 id="sec-1-3"><span class="section-number-3">1.3</span> Tuples and tensors </h3>167 <div class="outline-text-3" id="text-1-3">171 </div>173 <div id="outline-container-1-3-1" class="outline-4">174 <h4 id="sec-1-3-1"><span class="section-number-4">1.3.1</span> Tuples are “sequences with spin” </h4>175 <div class="outline-text-4" id="text-1-3-1">181 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)183 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">Let some objects have spin</span>185 (<span style="color: #f0dfaf; font-weight: bold;">defprotocol</span> <span style="color: #f0dfaf;">Spinning</span>186 (up? [this])187 (down? [this]))189 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">spin</span>190 <span style="color: #8fb28f;">"Returns the spin of the Spinning s, either :up or :down"</span>191 [<span style="color: #dfdfbf; font-weight: bold;">#^Spinning</span> s]192 (<span style="color: #f0dfaf; font-weight: bold;">cond</span> (up? s) <span style="color: #8cd0d3;">:up</span> (down? s) <span style="color: #8cd0d3;">:down</span>))195 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">DEFINITION: A tuple is a sequence with spin</span>197 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">Tuple</span>198 [spin coll]200 clojure.lang.Seqable201 (<span style="color: #8cd0d3;">seq</span> [this] (<span style="color: #8cd0d3;">seq</span> (.coll this)))203 clojure.lang.Counted204 (<span style="color: #8cd0d3;">count</span> [this] (<span style="color: #8cd0d3;">count</span> (.coll this)))206 Spinning207 (up? [this] (<span style="color: #8cd0d3;">=</span> <span style="color: #8cd0d3;">::up</span> (.spin this)))208 (down? [this] (<span style="color: #8cd0d3;">=</span> <span style="color: #8cd0d3;">::down</span> (.spin this))))210 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">print-method</span> Tuple211 [o w]212 (<span style="color: #8cd0d3;">print-simple</span>213 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (up? o)214 (<span style="color: #8cd0d3;">str</span> <span style="color: #cc9393;">"u"</span> (.coll o))215 (<span style="color: #8cd0d3;">str</span> <span style="color: #cc9393;">"d"</span> (<span style="color: #8cd0d3;">vec</span>(.coll o))))216 w))219 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">CONSTRUCTORS</span>221 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">up</span>222 <span style="color: #8fb28f;">"Create a new up-tuple containing the contents of coll."</span>223 [coll]224 (Tuple. <span style="color: #8cd0d3;">::up</span> coll))226 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">down</span>227 <span style="color: #8fb28f;">"Create a new down-tuple containing the contents of coll."</span>228 [coll]229 (Tuple. <span style="color: #8cd0d3;">::down</span> coll))230 </pre>237 </div>239 </div>241 <div id="outline-container-1-3-2" class="outline-4">242 <h4 id="sec-1-3-2"><span class="section-number-4">1.3.2</span> Matrices </h4>243 <div class="outline-text-4" id="text-1-3-2">248 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)249 (<span style="color: #8cd0d3;">require</span> 'incanter.core) <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">use incanter's fast matrices</span>252 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">all-equal?</span> [coll]253 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">empty?</span> (<span style="color: #8cd0d3;">rest</span> coll)) true254 (<span style="color: #f0dfaf; font-weight: bold;">and</span> (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">first</span> coll) (<span style="color: #8cd0d3;">second</span> coll))255 (<span style="color: #f0dfaf; font-weight: bold;">recur</span> (<span style="color: #8cd0d3;">rest</span> coll)))))258 (<span style="color: #f0dfaf; font-weight: bold;">defprotocol</span> <span style="color: #f0dfaf;">Matrix</span>259 (rows [matrix])260 (cols [matrix])261 (diagonal [matrix])262 (trace [matrix])263 (determinant [matrix])264 (transpose [matrix])265 (conjugate [matrix])266 )268 (<span style="color: #8cd0d3;">extend-protocol</span> Matrix incanter.Matrix269 (rows [rs] (<span style="color: #8cd0d3;">map</span> down (<span style="color: #8cd0d3;">apply</span> map vector (<span style="color: #8cd0d3;">apply</span> map vector rs))))270 (cols [rs] (<span style="color: #8cd0d3;">map</span> up (<span style="color: #8cd0d3;">apply</span> map vector rs)))271 (diagonal [matrix] (incanter.core/diag matrix) )272 (determinant [matrix] (incanter.core/det matrix))273 (trace [matrix] (incanter.core/trace matrix))274 (transpose [matrix] (incanter.core/trans matrix)))276 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">count-rows</span> [matrix]277 ((<span style="color: #8cd0d3;">comp</span> count rows) matrix))279 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">count-cols</span> [matrix]280 ((<span style="color: #8cd0d3;">comp</span> count cols) matrix))282 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">square?</span> [matrix]283 (<span style="color: #8cd0d3;">=</span> (count-rows matrix) (count-cols matrix)))285 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">identity-matrix</span>286 <span style="color: #8fb28f;">"Define a square matrix of size n-by-n with 1s along the diagonal and</span>287 <span style="color: #8fb28f;"> 0s everywhere else."</span>288 [n]289 (incanter.core/identity-matrix n))292 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">matrix-by-rows</span>293 <span style="color: #8fb28f;">"Define a matrix by giving its rows."</span>294 [& rows]295 (<span style="color: #f0dfaf; font-weight: bold;">if</span>296 (<span style="color: #8cd0d3;">not</span> (all-equal? (<span style="color: #8cd0d3;">map</span> count rows)))297 (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"All rows in a matrix must have the same number of elements."</span>))298 (incanter.core/matrix (<span style="color: #8cd0d3;">vec</span> rows))))300 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">matrix-by-cols</span>301 <span style="color: #8fb28f;">"Define a matrix by giving its columns"</span>302 [& cols]303 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not</span> (all-equal? (<span style="color: #8cd0d3;">map</span> count cols)))304 (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"All columns in a matrix must have the same number of elements."</span>))305 (incanter.core/matrix (<span style="color: #8cd0d3;">vec</span> (<span style="color: #8cd0d3;">apply</span> map vector cols)))))307 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">identity-matrix</span>308 <span style="color: #8fb28f;">"Define a square matrix of size n-by-n with 1s along the diagonal and</span>309 <span style="color: #8fb28f;"> 0s everywhere else."</span>310 [n]311 (incanter.core/identity-matrix n))313 </pre>318 </div>319 </div>321 </div>323 <div id="outline-container-1-4" class="outline-3">324 <h3 id="sec-1-4"><span class="section-number-3">1.4</span> Generic Operations </h3>325 <div class="outline-text-3" id="text-1-4">330 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)331 (<span style="color: #8cd0d3;">use</span> 'clojure.contrib.generic.arithmetic332 'clojure.contrib.generic.collection333 'clojure.contrib.generic.functor334 'clojure.contrib.generic.math-functions)336 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">numbers?</span>337 <span style="color: #8fb28f;">"Returns true if all arguments are numbers, else false."</span>338 [& xs]339 (<span style="color: #8cd0d3;">every?</span> number? xs))341 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">tuple-surgery</span>342 <span style="color: #8fb28f;">"Applies the function f to the items of tuple and the additional</span>343 <span style="color: #8fb28f;"> arguments, if any. Returns a Tuple of the same type as tuple."</span>344 [tuple f & xs]345 ((<span style="color: #f0dfaf; font-weight: bold;">if</span> (up? tuple) up down)346 (<span style="color: #8cd0d3;">apply</span> f (<span style="color: #8cd0d3;">seq</span> tuple) xs)))350 <span style="color: #708070;">;;; </span><span style="color: #7f9f7f;">CONTRACTION collapses two compatible tuples into a number.</span>352 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">contractible?</span>353 <span style="color: #8fb28f;">"Returns true if the tuples a and b are compatible for contraction,</span>354 <span style="color: #8fb28f;"> else false. Tuples are compatible if they have the same number of</span>355 <span style="color: #8fb28f;"> components, they have opposite spins, and their elements are</span>356 <span style="color: #8fb28f;"> pairwise-compatible."</span>357 [a b]358 (<span style="color: #f0dfaf; font-weight: bold;">and</span>359 (<span style="color: #8cd0d3;">isa?</span> (<span style="color: #8cd0d3;">type</span> a) Tuple)360 (<span style="color: #8cd0d3;">isa?</span> (<span style="color: #8cd0d3;">type</span> b) Tuple)361 (<span style="color: #8cd0d3;">=</span> (<span style="color: #8cd0d3;">count</span> a) (<span style="color: #8cd0d3;">count</span> b))362 (<span style="color: #8cd0d3;">not=</span> (spin a) (spin b))364 (<span style="color: #8cd0d3;">not-any?</span> false?365 (<span style="color: #8cd0d3;">map</span> #(<span style="color: #f0dfaf; font-weight: bold;">or</span>366 (numbers? %1 %2)367 (contractible? %1 %2))368 a b))))370 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">contract</span>371 <span style="color: #8fb28f;">"Contracts two tuples, returning the sum of the</span>372 <span style="color: #8fb28f;"> products of the corresponding items. Contraction is recursive on</span>373 <span style="color: #8fb28f;"> nested tuples."</span>374 [a b]375 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not</span> (contractible? a b))376 (<span style="color: #f0dfaf; font-weight: bold;">throw</span>377 (Exception. <span style="color: #cc9393;">"Not compatible for contraction."</span>))378 (<span style="color: #8cd0d3;">reduce</span> +379 (<span style="color: #8cd0d3;">map</span>380 (<span style="color: #8cd0d3;">fn</span> [x y]381 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (numbers? x y)382 (<span style="color: #8cd0d3;">*</span> x y)383 (contract x y)))384 a b))))390 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">conj</span> Tuple391 [tuple & xs]392 (tuple-surgery tuple #(<span style="color: #8cd0d3;">apply</span> conj % xs)))394 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">fmap</span> Tuple395 [f tuple]396 (tuple-surgery tuple (<span style="color: #8cd0d3;">partial</span> map f)))400 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">TODO: define Scalar, and add it to the hierarchy above Number and Complex</span>403 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [Tuple Tuple] <span style="color: #708070;">; </span><span style="color: #7f9f7f;">tuple*tuple</span>404 [a b]405 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (contractible? a b)406 (contract a b)407 (<span style="color: #8cd0d3;">map</span> (<span style="color: #8cd0d3;">partial</span> * a) b)))410 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [java.lang.Number Tuple] <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">scalar * tuple</span>411 [a x] (fmap (<span style="color: #8cd0d3;">partial</span> * a) x))413 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [Tuple java.lang.Number]414 [x a] (<span style="color: #8cd0d3;">*</span> a x))416 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [java.lang.Number incanter.Matrix] <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">scalar * matrix</span>417 [x M] (incanter.core/mult x M))419 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [incanter.Matrix java.lang.Number]420 [M x] (<span style="color: #8cd0d3;">*</span> x M))422 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [incanter.Matrix incanter.Matrix] <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">matrix * matrix</span>423 [M1 M2]424 (incanter.core/mmult M1 M2))426 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [incanter.Matrix Tuple] <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">matrix * tuple</span>427 [M v]428 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #f0dfaf; font-weight: bold;">and</span> (<span style="color: #8cd0d3;">apply</span> numbers? v) (up? v))429 (<span style="color: #8cd0d3;">*</span> M (matrix-by-cols v))430 (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"Currently, you can only multiply a matrix by a tuple of *numbers*"</span>))431 ))433 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [Tuple incanter.Matrix] <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">tuple * Matrix</span>434 [v M]435 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #f0dfaf; font-weight: bold;">and</span> (<span style="color: #8cd0d3;">apply</span> numbers? v) (down? v))436 (<span style="color: #8cd0d3;">*</span> (matrix-by-rows v) M)437 (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"Currently, you can only multiply a matrix by a tuple of *numbers*"</span>))438 ))441 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">exp</span> incanter.Matrix442 [M]443 (incanter.core/exp M))445 </pre>450 </div>451 </div>453 </div>455 <div id="outline-container-2" class="outline-2">456 <h2 id="sec-2"><span class="section-number-2">2</span> Operators and Differentiation </h2>457 <div class="outline-text-2" id="text-2">462 </div>464 <div id="outline-container-2-1" class="outline-3">465 <h3 id="sec-2-1"><span class="section-number-3">2.1</span> Operators </h3>466 <div class="outline-text-3" id="text-2-1">469 </div>471 </div>473 <div id="outline-container-2-2" class="outline-3">474 <h3 id="sec-2-2"><span class="section-number-3">2.2</span> Differential Terms and Sequences </h3>475 <div class="outline-text-3" id="text-2-2">480 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)481 (<span style="color: #8cd0d3;">use</span> 'clojure.contrib.seq482 'clojure.contrib.generic.arithmetic483 'clojure.contrib.generic.collection484 'clojure.contrib.generic.math-functions)486 <span style="color: #708070;">;;</span><span style="color: #7f9f7f;">∂</span>488 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">DEFINITION : Differential Term</span>490 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The</span>491 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">coefficient of the quantity is returned by the 'coefficient' method,</span>492 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">while the sequence of differential parameters is returned by the</span>493 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">method 'partials'.</span>495 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">Instead of using (potentially ambiguous) letters to denote</span>496 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].</span>498 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">The coefficient can be any arithmetic object; the</span>499 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">partials must be a nonrepeating sorted sequence of nonnegative</span>500 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">integers.</span>502 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">DifferentialTerm</span> [coefficient partials])504 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">differential-term</span>505 <span style="color: #8fb28f;">"Make a differential term given pairs of coefficients and partials."</span>506 [coefficient partials]507 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #f0dfaf; font-weight: bold;">and</span> (<span style="color: #8cd0d3;">coll?</span> partials) (<span style="color: #8cd0d3;">every?</span> #(<span style="color: #f0dfaf; font-weight: bold;">and</span> (<span style="color: #8cd0d3;">integer?</span> %) (<span style="color: #8cd0d3;">not</span>(<span style="color: #8cd0d3;">neg?</span> %))) partials))508 (DifferentialTerm. coefficient (<span style="color: #8cd0d3;">set</span> partials))509 (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (java.lang.IllegalArgumentException. <span style="color: #cc9393;">"Partials must be a collection of integers."</span>))))512 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">DEFINITION : Differential Sequence</span>513 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">A differential sequence is a sequence of differential terms, all with different partials.</span>514 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">Internally, it is a map from the partials of each term to their coefficients.</span>516 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">DifferentialSeq</span>517 [terms]518 <span style="color: #708070;">;;</span><span style="color: #7f9f7f;">clojure.lang.IPersistentMap</span>519 clojure.lang.Associative520 (<span style="color: #8cd0d3;">assoc</span> [this key val]521 (DifferentialSeq.522 (<span style="color: #8cd0d3;">cons</span> (differential-term val key) terms)))523 (<span style="color: #8cd0d3;">cons</span> [this x]524 (DifferentialSeq. (<span style="color: #8cd0d3;">cons</span> x terms)))525 (containsKey [this key]526 (<span style="color: #8cd0d3;">not</span>(<span style="color: #8cd0d3;">nil?</span> (find-first #(<span style="color: #8cd0d3;">=</span> (.partials %) key) terms))))527 (<span style="color: #8cd0d3;">count</span> [this] (<span style="color: #8cd0d3;">count</span> (.terms this)))528 (<span style="color: #8cd0d3;">empty</span> [this] (DifferentialSeq. []))529 (entryAt [this key]530 ((<span style="color: #8cd0d3;">juxt</span> #(.partials %) #(.coefficient %))531 (find-first #(<span style="color: #8cd0d3;">=</span> (.partials %) key) terms)))532 (<span style="color: #8cd0d3;">seq</span> [this] (<span style="color: #8cd0d3;">seq</span> (.terms this))))535 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">fmap</span> DifferentialSeq536 [f dseq]537 (DifferentialSeq. (fmap f (.terms dseq))))540 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">BUILDING DIFFERENTIAL OBJECTS</span>542 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">differential-seq</span>543 <span style="color: #8fb28f;">"Define a differential sequence by specifying coefficient/partials</span>544 <span style="color: #8fb28f;">pairs, which are used to create differential terms to populate the</span>545 <span style="color: #8fb28f;">sequence."</span>546 ([coefficient partials]547 (DifferentialSeq. {(<span style="color: #8cd0d3;">set</span> partials) coefficient}))548 ([coefficient partials & cps]549 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">odd?</span> (<span style="color: #8cd0d3;">count</span> cps))550 (<span style="color: #f0dfaf; font-weight: bold;">throw</span> (Exception. <span style="color: #cc9393;">"differential-seq requires an even number of terms."</span>))551 (DifferentialSeq.552 (<span style="color: #8cd0d3;">reduce</span>553 #(<span style="color: #8cd0d3;">assoc</span> %1 (<span style="color: #8cd0d3;">set</span> (<span style="color: #8cd0d3;">second</span> %2)) (<span style="color: #8cd0d3;">first</span> %2))554 {(<span style="color: #8cd0d3;">set</span> partials) coefficient}555 (<span style="color: #8cd0d3;">partition</span> 2 cps))))))559 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">big-part</span>560 <span style="color: #8fb28f;">"Returns the part of the differential sequence that is finite,</span>561 <span style="color: #8fb28f;"> i.e. not infinitely small, as a differential sequence. If given a</span>562 <span style="color: #8fb28f;"> non-differential object, returns that object."</span>563 [dseq]564 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not=</span> (<span style="color: #8cd0d3;">type</span> dseq) DifferentialSeq) dseq565 (<span style="color: #f0dfaf; font-weight: bold;">let</span> [m (.terms dseq)566 keys (<span style="color: #8cd0d3;">sort-by</span> count (<span style="color: #8cd0d3;">keys</span> m))567 smallest-var (<span style="color: #8cd0d3;">last</span> (<span style="color: #8cd0d3;">last</span> keys))]568 (DifferentialSeq.569 (<span style="color: #8cd0d3;">reduce</span>570 #(<span style="color: #8cd0d3;">assoc</span> %1 (<span style="color: #8cd0d3;">first</span> %2) (<span style="color: #8cd0d3;">second</span> %2))571 {}572 (<span style="color: #8cd0d3;">remove</span> #((<span style="color: #8cd0d3;">first</span> %) smallest-var) m))))))575 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">small-part</span>576 <span style="color: #8fb28f;">"Returns the part of the differential sequence that infinitely</span>577 <span style="color: #8fb28f;"> small, as a differential sequence. If given a non-differential</span>578 <span style="color: #8fb28f;"> object, returns that object."</span>579 [dseq]580 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not=</span> (<span style="color: #8cd0d3;">type</span> dseq) DifferentialSeq) 0581 (<span style="color: #f0dfaf; font-weight: bold;">let</span> [m (.terms dseq)582 keys (<span style="color: #8cd0d3;">sort-by</span> count (<span style="color: #8cd0d3;">keys</span> m))583 smallest-var (<span style="color: #8cd0d3;">last</span> (<span style="color: #8cd0d3;">last</span> keys))]584 (DifferentialSeq.585 (<span style="color: #8cd0d3;">reduce</span>586 #(<span style="color: #8cd0d3;">assoc</span> %1 (<span style="color: #8cd0d3;">first</span> %2) (<span style="color: #8cd0d3;">second</span> %2)) {}587 (<span style="color: #8cd0d3;">filter</span> #((<span style="color: #8cd0d3;">first</span> %) smallest-var) m))))))591 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">cartesian-product</span> [set1 set2]592 (<span style="color: #8cd0d3;">reduce</span> concat593 (<span style="color: #f0dfaf; font-weight: bold;">for</span> [x set1]594 (<span style="color: #f0dfaf; font-weight: bold;">for</span> [y set2]595 [x y]))))597 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">differential-multiply</span>598 <span style="color: #8fb28f;">"Multiply two differential sequences. The square of any differential</span>599 <span style="color: #8fb28f;"> variable is zero since differential variables are infinitesimally</span>600 <span style="color: #8fb28f;"> small."</span>601 [dseq1 dseq2]602 (DifferentialSeq.603 (<span style="color: #8cd0d3;">reduce</span>604 (<span style="color: #8cd0d3;">fn</span> [m [[vars1 coeff1] [vars2 coeff2]]]605 (<span style="color: #f0dfaf; font-weight: bold;">if</span> (<span style="color: #8cd0d3;">not</span> (<span style="color: #8cd0d3;">empty?</span> (clojure.set/<span style="color: #dfdfbf; font-weight: bold;">intersection</span> vars1 vars2)))606 m607 (<span style="color: #8cd0d3;">assoc</span> m (clojure.set/<span style="color: #dfdfbf; font-weight: bold;">union</span> vars1 vars2) (<span style="color: #8cd0d3;">*</span> coeff1 coeff2))))608 {}609 (cartesian-product (.terms dseq1) (.terms dseq2)))))613 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [DifferentialSeq DifferentialSeq]614 [dseq1 dseq2]615 (differential-multiply dseq1 dseq2))617 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">+</span> [DifferentialSeq DifferentialSeq]618 [dseq1 dseq2]619 (DifferentialSeq.620 (<span style="color: #8cd0d3;">merge-with</span> + (.terms dseq1) (.terms dseq2))))622 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [java.lang.Number DifferentialSeq]623 [x dseq]624 (fmap (<span style="color: #8cd0d3;">partial</span> * x) dseq))626 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">*</span> [DifferentialSeq java.lang.Number]627 [dseq x]628 (fmap (<span style="color: #8cd0d3;">partial</span> * x) dseq))629 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">DERIVATIVES</span>633 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">linear-approximator</span>634 <span style="color: #8fb28f;">"Returns an operator that linearly approximates the given function."</span>635 ([f df|dx]636 (<span style="color: #8cd0d3;">fn</span> [x]637 (<span style="color: #f0dfaf; font-weight: bold;">let</span> [big-part (big-part x)638 small-part (small-part x)]639 (<span style="color: #8cd0d3;">+</span> (fmap f big-part)640 (<span style="color: #8cd0d3;">*</span> (fmap df|dx big-part) small-part)))))641 ([f df|dx df|dy]642 (<span style="color: #8cd0d3;">fn</span> [x y]643 (<span style="color: #f0dfaf; font-weight: bold;">let</span> [X (big-part x)644 Y (big-part y)645 DX (small-part x)646 DY (small-part y)]647 (<span style="color: #8cd0d3;">+</span> (f X Y)648 (<span style="color: #8cd0d3;">+</span> (<span style="color: #8cd0d3;">*</span> DX (df|dx X Y))649 (<span style="color: #8cd0d3;">*</span> DY (df|dy X Y)))))))650 )656 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">D</span>[f]657 (<span style="color: #8cd0d3;">fn</span>[x]658 (f (differential-seq x [] 1 [0]))))660 (<span style="color: #f0dfaf; font-weight: bold;">defn</span> <span style="color: #f0dfaf;">d</span>[f partials]661 (<span style="color: #8cd0d3;">fn</span> [x]662 (<span style="color: #8cd0d3;">get</span>663 (.terms ((D f)x))664 (<span style="color: #8cd0d3;">set</span> partials)665 0666 )))668 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">exp</span> DifferentialSeq [x]669 ((linear-approximator exp exp) x))671 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">sin</span> DifferentialSeq672 [x]673 ((linear-approximator sin cos) x))675 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">cos</span> DifferentialSeq676 [x]677 ((linear-approximator cos #(<span style="color: #8cd0d3;">-</span> (sin %))) x))681 </pre>691 </div>692 </div>693 </div>694 <div id="postamble">695 <p class="date">Date: 2011-08-11 04:10:36 EDT</p>696 <p class="author">Author: Dylan Holmes</p>697 <p class="creator">Org version 7.6 with Emacs version 23</p>698 <a href="http://validator.w3.org/check?uri=referer">Validate XHTML 1.0</a>699 </div>700 </div>701 </body>702 </html>