Mercurial > dylan

     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>

    80 

    81 </head>

    82 <body>

    83 

    84 <div id="content">

    85 

    86 

    87 

    88 <div class="header">

    89   <div class="float-right">	

    90     <!-- 

    91     <form>

    92       <input type="text"/><input type="submit" value="search the blog &raquo;"/> 

    93     </form>

    94     -->

    95   </div>

    96 

    97   <h1>aurellem <em>&#x2609;</em></h1>

    98   <ul class="nav">

    99     <li><a href="/">read the blog &raquo;</a></li>

   100     <!-- li><a href="#">learn about us &raquo;</a></li-->

   101   </ul>

   102 </div>

   103 

   104 <h1 class="title">Building a Classical Mechanics Library in Clojure</h1>

   105 

   106 

   107 

   108 

   109 

   110 

   111 

   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 &ldquo;sequences with spin&rdquo; </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>

   138 

   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">

   142 

   143 

   144 

   145 </div>

   146 

   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">

   150 

   151 

   152 </div>

   153 

   154 </div>

   155 

   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">

   159 

   160 

   161 </div>

   162 

   163 </div>

   164 

   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">

   168 

   169 

   170 

   171 </div>

   172 

   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 &ldquo;sequences with spin&rdquo; </h4>

   175 <div class="outline-text-4" id="text-1-3-1">

   176 

   177 

   178 

   179 

   180 

   181 <pre class="src src-clojure">(<span style="color: #f0dfaf; font-weight: bold;">in-ns</span> 'sicm.utils)

   182 

   183 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">Let some objects have spin</span>

   184 

   185 (<span style="color: #f0dfaf; font-weight: bold;">defprotocol</span> <span style="color: #f0dfaf;">Spinning</span>

   186   (up? [this])

   187   (down? [this]))

   188 

   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>))

   193 

   194 

   195 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">DEFINITION: A tuple is a sequence with spin</span>

   196 

   197 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">Tuple</span>

   198   [spin coll]

   199 

   200   clojure.lang.Seqable

   201   (<span style="color: #8cd0d3;">seq</span> [this] (<span style="color: #8cd0d3;">seq</span> (.coll this)))

   202 

   203   clojure.lang.Counted

   204   (<span style="color: #8cd0d3;">count</span> [this] (<span style="color: #8cd0d3;">count</span> (.coll this)))

   205 

   206   Spinning

   207   (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))))

   209 

   210 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">print-method</span> Tuple

   211   [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))

   217 

   218 

   219 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">CONSTRUCTORS</span>

   220 

   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))       

   225 

   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>

   231 

   232 

   233 

   234 

   235 

   236 

   237 </div>

   238 

   239 </div>

   240 

   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">

   244 

   245 

   246 

   247 

   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>

   250 

   251 

   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)) true

   254       (<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)))))

   256 

   257 

   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 )

   267 

   268 (<span style="color: #8cd0d3;">extend-protocol</span> Matrix incanter.Matrix

   269   (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)))

   275 

   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))

   278 

   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))

   281 

   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)))

   284 

   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))

   290 

   291 

   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   [&amp; 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))))

   299 

   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   [&amp; 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)))))

   306 

   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))

   312 

   313 </pre>

   314 

   315 

   316 

   317 

   318 </div>

   319 </div>

   320 

   321 </div>

   322 

   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">

   326 

   327 

   328 

   329 

   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.arithmetic

   332      'clojure.contrib.generic.collection

   333      'clojure.contrib.generic.functor

   334      'clojure.contrib.generic.math-functions)

   335 

   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   [&amp; xs]

   339   (<span style="color: #8cd0d3;">every?</span> number? xs))

   340 

   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 &amp; 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)))

   347 

   348 

   349 

   350 <span style="color: #708070;">;;; </span><span style="color: #7f9f7f;">CONTRACTION collapses two compatible tuples into a number.</span>

   351 

   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))

   363    

   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))))

   369 

   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))))

   385 

   386 

   387 

   388 

   389 

   390 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">conj</span> Tuple

   391   [tuple &amp; xs]

   392   (tuple-surgery tuple #(<span style="color: #8cd0d3;">apply</span> conj % xs)))

   393 

   394 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">fmap</span> Tuple

   395   [f tuple]

   396   (tuple-surgery tuple (<span style="color: #8cd0d3;">partial</span> map f)))

   397 

   398 

   399 

   400 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">TODO: define Scalar, and add it to the hierarchy above Number and Complex</span>

   401 

   402                                         

   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)))

   408 

   409 

   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))

   412 

   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))

   415 

   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))

   418 

   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))

   421 

   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))

   425 

   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     ))

   432 

   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     ))

   439 

   440 

   441 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">exp</span> incanter.Matrix

   442   [M]

   443   (incanter.core/exp M))

   444 

   445 </pre>

   446 

   447 

   448 

   449 

   450 </div>

   451 </div>

   452 

   453 </div>

   454 

   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">

   458 

   459 

   460 

   461 

   462 </div>

   463 

   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">

   467 

   468 

   469 </div>

   470 

   471 </div>

   472 

   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">

   476 

   477 

   478 

   479 

   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.seq

   482      'clojure.contrib.generic.arithmetic

   483      'clojure.contrib.generic.collection

   484      'clojure.contrib.generic.math-functions)

   485 

   486 <span style="color: #708070;">;;</span><span style="color: #7f9f7f;">&#8706;</span>

   487 

   488 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">DEFINITION : Differential Term</span>

   489 

   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>

   494 

   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>

   497 

   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>

   501 

   502 (<span style="color: #f0dfaf; font-weight: bold;">deftype</span> <span style="color: #f0dfaf;">DifferentialTerm</span> [coefficient partials])

   503 

   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>))))

   510 

   511 

   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>

   515 

   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.Associative

   520   (<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))))

   533 

   534 

   535 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">fmap</span> DifferentialSeq

   536   [f dseq]

   537   (DifferentialSeq. (fmap f (.terms dseq))))

   538 

   539 

   540 <span style="color: #708070;">;; </span><span style="color: #7f9f7f;">BUILDING DIFFERENTIAL OBJECTS</span>

   541 

   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 &amp; 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))))))

   556   

   557 

   558 

   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) dseq

   565       (<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))))))

   573 

   574 

   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) 0

   581       (<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))))))

   588 

   589 

   590 

   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> concat

   593           (<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]))))

   596 

   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         m

   607         (<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)))))

   610 

   611 

   612 

   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))

   616 

   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))))

   621 

   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))

   625 

   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>

   630 

   631 

   632 

   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   )

   651 

   652 

   653 

   654 

   655 

   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]))))

   659 

   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      0

   666     )))

   667 

   668 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">exp</span> DifferentialSeq [x]

   669            ((linear-approximator exp exp) x))

   670 

   671 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">sin</span> DifferentialSeq

   672   [x]

   673   ((linear-approximator sin cos) x))

   674 

   675 (<span style="color: #f0dfaf; font-weight: bold;">defmethod</span> <span style="color: #f0dfaf;">cos</span> DifferentialSeq

   676   [x]

   677   ((linear-approximator cos #(<span style="color: #8cd0d3;">-</span> (sin %))) x))

   678      

   679 

   680 

   681 </pre>

   682 

   683 

   684 

   685 

   686 

   687 

   688 

   689 

   690 

   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>
author	Dylan Holmes <ocsenave@gmail.com>
date	Thu, 27 Oct 2011 23:02:59 -0500
parents	b4de894a1e2e
children