# HG changeset patch # User Robert McIntyre # Date 1319785385 25200 # Node ID b4de894a1e2e6d6b64abf5f2f1b05a877c46fe97 # Parent 8d8278e09888fe25a7cf33164a849be8e44b71f6 initial import diff -r 8d8278e09888 -r b4de894a1e2e categorical/imperative.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/imperative.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,1207 @@ + + + + +LAMBDA: The Ultimate Imperative + + + + + + + + + + + + +

+ + + +

1 Abstract
2 Introduction +
- 2.1 Simple Loops +
  - 2.1.1 Simple Recursion
  - 2.1.2 Iteration
+
3 Imperative Programming +
+
4 Continuations +
- 4.1 Continuation-Passing Recursion
- 4.2 Escape Expressions
+

+ +

1 Abstract

+ + +

+We demonstrate how to model the following common programming constructs in +terms of an applicative order language similar to LISP: +

Simple Recursion +
Iteration +
Compound Statements and Expressions +
GO TO and Assignment +
Continuation—Passing +
Escape Expressions +
Fluid Variables +
Call by Name, Call by Need, and Call by Reference +

+ +

The models require only (possibly self-referent) lambda application, +conditionals, and (rarely) assignment. No complex data structures such as +stacks are used. The models are transparent, involving only local syntactic +transformations. +Some of these models, such as those for GO TO and assignment, are already well +known, and appear in the work of Landin, Reynolds, and others. The models for +escape expressions, fluid variables, and call by need with side effects are +new. This paper is partly tutorial in intent, gathering all the models +together for purposes of context. +This report describes research done at the Artificial Intelligence Laboratory +of the Massachusetts Institute of Teehnology. Support for the laboratory's +artificial intelligence research is provided in part by the Advanced Research +Projects Agency of the Department of Defense under Office of Naval Research +contract N000l4-75-C-0643. +

+ +

2 Introduction

+ + +

+ We catalogue a number of common programming constructs. For each +construct we examine "typical" usage in well-known programming languages, and +then capture the essence of the semantics of the construct in terms of a +common meta-language. +The lambda calculus {Note Alonzowins} is often used as such a meta- +language. Lambda calculus offers clean semantics, but it is clumsy because it +was designed to be a minimal language rather than a convenient one. All +lambda calculus "functions" must take exactly one "argument"; the only "data +type" is lambda expressions; and the only "primitive operation‘ is variable‘ +substitution. While its utter simplicity makes lambda calculus ideal for +logicians, it is too primitive for use by programmers. The meta-language we +use is a programming language called SCHEME {Note Schemepaper) which is based +on lambda calculus. +SCHEME is a dialect of LISP. [McCarthy 62] It is an expression- +oriented, applicative order, interpreter-based language which allows one to +manipulate programs as data. It differs from most current dialects of LISP in +that it closes all lambda expressions in the environment of their definition +or declaration, rather than in the execution environment. {Note Closures} +This preserves the substitution semantics of lambda calculus, and has the +consequence that all variables are lexically scoped, as in ALGOL. [Naur 63] +Another difference is that SCHEME is implemented in such a way that tail- +recursions execute without net growth of the interpreter stack. {Note +Schemenote} We have chosen to use LISP syntax rather than, say, ALGOL syntax +because we want to treat programs as data for the purpose of describing +transformations on the code. LISP supplies names for the parts of an +executable expression and standard operators for constructing expressions and +extracting their components. The use of LISP syntax makes the structure of +such expressions manifest. We use ALGOL as an expository language, because it +is familiar to many people, but ALGOL is not sufficiently powerful to express +the necessary concepts; in particular, it does not allow functions to return +functions as values. we are thus forced to use a dialect of LISP in many +cases. +We will consider various complex programming language constructs and +show how to model them in terms of only a few simple ones. As far as possible +we will use only three control constructs from SCHEME: LAMBDA expressions, as +in LISP, which are Just functions with lexically scoped free variables; +LABELS, which allows declaration of mutually recursive procedures (Note +Labelsdef}; and IF, a primitive conditional expression. For more complex +modelling we will introduce an assignment primitive (ASET).i We will freely +assume the existence of other comon primitives, such as arithmetic functions. +The constructs we will examine are divided into four broad classes. +The first is sﬁnph?Lo0pl; this contains simple recursions and iterations, and +an introduction to the notion of continuations. The second is hmponﬂivo +Connrucls; this includes compound statements, G0 T0, and simple variable +assignments. The third is continuations, which encompasses the distinction between statements and expressions, escape operators (such as Landin's J- +operator [Landin 65] and Reynold's escape expression [Reynolds 72]). and fluid +(dynamically bound) variables. The fourth is Parameter Passing Mechanisms, such +as ALGOL call—by-name and FORTRAN call-by-location. +Some of the models presented here are already well-known, particularly +those for G0 T0 and assignment. [McCarthy 60] [Landin 65] [Reynolds 72] +Those for escape operators, fluid variables, and call-by-need with side +effects are new. +

+ +

2.1 Simple Loops

+ +

By simple loops we mean constructs which enable programs to execute the +same piece of code repeatedly in a controlled manner. Variables may be made +to take on different values during each repetition, and the number of +repetitions may depend on data given to the program. +

+ +

2.1.1 Simple Recursion

+ +

One of the easiest ways to produce a looping control structure is to +use a recursive function, one which calls itself to perform a subcomputation. +For example, the familiar factorial function may be written recursively in +ALGOL: ' +

+ + +

integer procedure fact(n) value n: integer n:
+fact := if n=0 then 1 else n*fact(n-1);
+
+

+ + + + +

+The invocation fact(n) computes the product of the integers from 1 to n using +the identity n!=n(n-1)! (n>0). If $n$ is zero, 1 is returned; otherwise fact: +calls itself recursively to compute $(n-1)!$ , then multiplies the result by $n$ +and returns it. +

+This same function may be written in Clojure as follows: +

+ + +

(ns categorical.imperative)
+(defn fact[n]
+  (if (= n 0) 1 (* n (fact (dec n))))
+)
+
+

+ + + + + +

+Clojure does not require an assignment to the ``variable'' fan: to return a value +as ALGOL does. The IF primitive is the ALGOL if-then-else rendered in LISP +syntax. Note that the arithmetic primitives are prefix operators in SCHEME. +

+ +

2.1.2 Iteration

+ +

There are many other ways to compute factorial. One important way is +through the use of iteration. +A comon iterative construct is the DO loop. The most general form we +have seen in any programming language is the MacLISP DO [Moon 74]. It permits +the simultaneous initialization of any number of control variables and the +simultaneous stepping of these variables by arbitrary functions at each +iteration step. The loop is terminated by an arbitrary predicate, and an +arbitrary value may be returned. The DO loop may have a body, a series of +expressions executed for effect on each iteration. A version of the MacLISP +DO construct has been adopted in SCHEME. +

+The general form of a SCHEME DO is: +

+ + +

(DO ((<var1> (init1> <stepl>)
+((var2> <init2> <step2>)
+(<varn> <intn> (stepn>))
+(<pred> (value>)
+(optional body>)
+
+

+ + + +

+The semantics of this are that the variables are bound and initialized to the +values of the <initi> expressions, which must all be evaluated in the +environment outside the D0; then the predicate <pred> is evaluated in the new +environment, and if TRUE, the (value) is evaluated and returned. Otherwise +the (optional body) is evaluated, then each of the steppers <stepi> is +evaluated in the current environment, all the variables made to have the +results as their values, the predicate evaluated again, and so on. +Using D0 loops in both ALGOL and SCHEME, we may express FACT by means +of iteration. +

+ + +

integer procedure fact(n); value n; integer n:
+begin
+integer m, ans;
+ans := 1;
+for m := n step -l until 0 do ans := m*ans;
+fact := ans;
+end;
+
+

+ + + + + +

(in-ns 'categorical.imperative)
+(defn fact-do[n]
+)
+
+

+ + + + +

+Note that the SCHEME D0 loop in FACT has no body – the stepping functions do +all the work. The ALGOL DO loop has an assignment in its body; because an +ALGOL DO loop can step only one variable, we need the assignment to step the +the variable "manually'. +In reality the SCHEME DO construct is not a primitive; it is a macro +which expands into a function which performs the iteration by tail—recursion. +Consider the following definition of FACT in SCHEME. Although it appears to +be recursive, since it "calls itself", it is entirely equivalent to the DO +loop given above, for it is the code that the DO macro expands into! It +captures the essence of our intuitive notion of iteration, because execution +of this program will not produce internal structures (e.g. stacks or variable +bindings) which increase in size with the number of iteration steps. +

+ + + +

(in-ns 'categorical.imperative)
+(defn fact-do-expand[n]
+  (let [fact1
+        (fn[m ans]
+          (if (zero? m) ans
+              (recur (dec m) (* m ans))))]
+    (fact1 n 1)))
+
+

+ + + + +

+From this we can infer a general way to express iterations in SCHEME in +a manner isomorphic to the HacLISP D0. The expansion of the general D0 loop +

+ + +

(DO ((<varl> <init1> <step1>)
+(<var2> (init2) <step2>)
+...
+(<varn> <initn> <stepn>))
+(<pred> <va1ue>)
+<body>)
+
+

+ + + +

+is this: +

+ + +

(let [doloop
+(fn [dummy <var1> <var2> ... <varn>)
+(if <pred> <value>
+(recur <body> <step1> <step2> ... <stepn>)))]
+
+(doloop nil <init1> <init2> ... <initn>))
+
+

+ + + +

+The identifiers doloop and dummy are chosen so as not to conflict with any +other identifiers in the program. +Note that, unlike most implementations of D0, there are no side effects +in the steppings of the iteration variables. D0 loops are usually modelled +using assignment statements. For example: +

+ + +

for r :1 0 step b until c do <statement>;
+
+

+ + + + +

+can be modelled as follows: [Naur 63] +

+ + +

begin
+x := a;
+L: if (x-c)*sign(n) > 0 then go to Endloop;
+<statement>;
+x := x+b;
+go to L:
+Endloop:
+end;
+
+

+ + + +

+Later we will see how such assignment statements can in general be +modelled without using side effects. +

+ +

3 Imperative Programming

+ +

Lambda calculus (and related languages, such as ``pure LISP'') is often +used for modelling the applicative constructs of programming languages. +However, they are generally thought of as inappropriate for modelling +imperative constructs. In this section we show how imperative constructs may +be modelled by applicative SCHEME constructs. +

+ +

3.1 Compound Statements

+ +

The simplest kind of imperative construct is the statement sequencer, +for example the compound statement in ALGOL: +

+ + +

begin
+S1;
+S2;
+end
+
+

+ + + + +

+This construct has two interesting properties: +

(1) It performs statement S1 before S2, and so may be used for + sequencing. +
(2) S1 is useful only for its side effects. (In ALGOL, S2 must also + be a statement, and so is also useful only for side effects, but + other languages have compound expressions containing a statement + followed by an expression.) +

+ + +

+The ALGOL compound statement may actually contain any number of statements, +but such statements can be expressed as a series of nested two-statement +compounds. That is: +

+ + +

begin
+S1;
+S2;
+...
+Sn-1;
+Sn;
+end
+
+

+ + + +

+is equivalent to: +

+ + + +

begin
+S1;
+begin
+S2;
+begin
+...
+
+begin
+Sn-1;
+Sn;
+end;
+end;
+end:
+end
+
+

+ + + + +

+It is not immediately apparent that this sequencing can be expressed in a +purely applicative language. We can, however, take advantage of the implicit +sequencing of applicative order evaluation. Thus, for example, we may write a +two-statement sequence as follows: +

+ + + +

((fn[dummy] S2) S1)
+
+

+ + + + +

+where dummy is an identifier not used in S2. From this it is +manifest that the value of S1 is ignored, and so is useful only for +side effects. (Note that we did not claim that S1 is expressed in a +purely applicative language, but only that the sequencing can be so +expressed.) From now on we will use the form (BLOCK S1 S2) as an +abbreviation for this expression, and (BLOCK S1 S2...Sn-1 Sn) as an +abbreviation for +

+ + + +

(BLOCK S1 (BLOCK S2 (BLOCK ... (BLOCK Sn-1 Sn)...))) 
+
+

+ + + + +

+ +

3.2 2.2. The G0 TO Statement

+ + +

+A more general imperative structure is the compound statement with +labels and G0 T0s. Consider the following code fragment due to +Jacopini, taken from Knuth: [Knuth 74] +

+ + +

begin 
+L1: if B1 then go to L2; S1;
+    if B2 then go to L2; S2;
+    go to L1;
+L2: S3;
+
+

+ + + + +

+It is perhaps surprising that this piece of code can be syntactically +transformed into a purely applicative style. For example, in SCHEME we +could write: +

+ + + +

(in-ns 'categorical.imperative)
+(let
+    [L1 (fn[]
+          (if B1 (L2) (do S1
+                          (if B2 (L2) (do S2 (L1))))))
+     L2 (fn[] S3)]
+  (L1))
+
+

+ + + + +

+As with the D0 loop, this transformation depends critically on SCHEME's +treatment of tail-recursion and on lexical scoping of variables. The labels +are names of functions of no arguments. In order to ‘go to" the labeled code, +we merely call the function named by that label. +

+ +

3.3 2.3. Simple Assignment

+ +

Of course, all this sequencing of statements is useless unless the +statements have side effects. An important side effect is assignment. For +example, one often uses assignment to place intermediate results in a named +location (i.e. a variable) so that they may be used more than once later +without recomputing them: +

+ + + +

begin
+real a2, sqrtdisc;
+a2 := 2*a;
+sqrtdisc := sqrt(b^2 - 4*a*c);
+root1 := (- b + sqrtdisc) / a2;
+root2 := (- b - sqrtdisc) / a2;
+print(root1);
+print(root2);
+print(root1 + root2);
+end
+
+

+ + + + +

+It is well known that such naming of intermediate results may be accomplished +by calling a function and binding the formal parameter variables to the +results: +

+ + + +

(fn [a2 sqrtdisc]
+  ((fn[root1 root2]
+    (do (println root1)
+        (println root2)
+        (println (+ root1 root2))))
+  (/ (+ (- b) sqrtdisc) a2)
+  (/ (- (- b) sqrtdisc) a2))
+
+  (* 2 a)
+  (sqrt (- (* b b) (* 4 a c))))
+
+

+ + + +

+This technique can be extended to handle all simple variable assignments which +appear as statements in blocks, even if arbitrary G0 T0 statements also appear +in such blocks. {Note Mccarthywins} +

+For example, here is a program which uses G0 TO statements in the form +presented before; it determines the parity of a non-negative integer by +counting it down until it reaches zero. +

+ + + +

begin
+L1: if a = 0 then begin parity := 0; go to L2; end;
+    a := a - 1;
+    if a = 0 then begin parity := 1; go to L2; end;
+    a := a - 1;
+    go to L1;
+L2: print(parity);
+
+

+ + + + +

+This can be expressed in SCHEME: +

+ + + +

(let
+    [L1 (fn [a parity]
+          (if (zero? a) (L2 a 0)
+              (L3 (dec a) parity)))
+     L3 (fn [a parity]
+          (if (zero? a) (L2 a 1)
+              (L1 (dec a) parity)))
+     L2 (fn [a parity]
+          (println parity))]
+  (L1 a parity))
+
+

+ + + + +

+The trick is to pass the set of variables which may be altered as arguments to +the label functions. {Note Flowgraph} It may be necessary to introduce new +labels (such as L3) so that an assignment may be transformed into the binding +for a function call. At worst, one may need as many labels as there are +statements (not counting the eliminated assignment and GO TO statements). +

+ +

3.4 Compound Expressions '

+ +

At this point we are almost in a position to model the most general +form of compound statement. In LISP, this is called the 'PROG feature". In +addition to statement sequencing and G0 T0 statements, one can return a value +from a PROG by using the RETURN statement. +

+Let us first consider the simplest compound statement, which in SCHEME +we call BLOCK. Recall that +(BLOCK s1 sz) is defined to be ((lambda (dummy) s2) s1) +

+Notice that this not only performs Sl before S2, but also returns the value of +S2. Furthermore, we defined (BLOCK S1 S2 ... Sn) so that it returns the value +of Sn. Thus BLOCK may be used as a compound expression, and models a LISP +PROGN, which is a PROG with no G0 T0 statements and an implicit RETURN of the +last ``statement'' (really an expression). +

+Host LISP compilers compile D0 expressions by macro-expansion. We have +already seen one way to do this in SCHEME using only variable binding. A more +common technique is to expand the D0 into a PROG, using variable assignments +instead of bindings. Thus the iterative factorial program: +

+ + + +

(oarxnz FACT .
+(LAMaoA (n) .
+(D0 ((M N (- H 1))
+(Ans 1 (* M Ans)))
+((- H 0) A"$))))
+
+

+ + + + +

+would expand into: +

+ + + +

(DEFINE FACT
+. (LAMeoA (M)
+(PRO6 (M Ans)
+(sszro M n
+Ans 1) ~
+LP (tr (- M 0) (RETURN Ans))
+(ssero M (- n 1)
+Ans (' M Ans))
+(60 LP))))
+
+

+ + + + +

+where SSETQ is a simultaneous multiple assignment operator. (SSETQ is not a +SCHEME (or LISP) primitive. It can be defined in terms of a single assignment +operator, but we are more interested here in RETURN than in simultaneous +assignment. The SSETQ's will all be removed anyway and modelled by lambda +binding.) We can apply the same technique we used before to eliminate G0 T0 +statements and assignments from compound statements: +

+ + + +

(DEFINE FACT
+(LAHBOA (I)
+(LABELS ((L1 (LAManA (M Ans)
+(LP n 1)))
+(LP (LAMaoA (M Ans)
+(IF (- M o) (nztunn Ans)
+(£2 H An$))))
+(L2 (LAMaoA (M An )
+(LP (- " 1) (' H ﬂN$)))))
+(L1 NIL NlL))))
+
+

+ + +

+ clojure +

+We still haven't done anything about RETURN. Let's see… +

==> the value of (FACT 0) is the value of (Ll NIL NIL) +
==> which is the value of (LP 0 1) +
==> which is the value of (IF (= 0 0) (RETURN 1) (L2 0 1)) +
==> which is the value of (RETURN 1) +

+ + +

+Notice that if RETURN were the identity function (LAMBDA (X) X), we would get the right answer. This is in fact a +general truth: if we Just replace a call to RETURN with its argument, then +our old transformation on compound statements extends to general compound +expressions, i.e. PROG. +

+ +

4 Continuations

+ +

Up to now we have thought of SCHEME's LAMBDA expressions as functions, +and of a combination such as (G (F X Y)) as meaning ``apply the function F to +the values of X and Y, and return a value so that the function G can be +applied and return a value …'' But notice that we have seldom used LAMBDA +expressions as functions. Rather, we have used them as control structures and +environment modifiers. For example, consider the expression: +(BLOCK (PRINT 3) (PRINT 4)) +

+This is defined to be an abbreviation for: +((LAMBDA (DUMMY) (PRINT 4)) (PRINT 3)) +

+We do not care about the value of this BLOCK expression; it follows that we +do not care about the value of the (LAMBDA (DUMMY) ...). We are not using +LAMBDA as a function at all. +

+It is possible to write useful programs in terms of LAHBDA expressions +in which we never care about the value of any lambda expression. We have +already demonstrated how one could represent any "FORTRAN" program in these +terms: all one needs is PROG (with G0 and SETQ), and PRINT to get the answers +out. The ultimate generalization of this imperative programing style is +continuation-passing. (Note Churchwins} . +

+ +

4.1 Continuation-Passing Recursion

+ +

Consider the following alternative definition of FACT. It has an extra +argument, the continuation, which is a function to call with the answer, when +we have it, rather than return a value +

+ + + +

procedure Inc!(n, c); value n, c;
+integer n: procedure c(integer value);
+if n-=0 then c(l) else
+begin
+procedure l(!mp(d) value a: integer a;
+c(mm);
+Iacl(n-1, romp):
+end;
+
+

+ + + + + +

(in-ns 'categorical.imperative)
+(defn fact-cps[n k]
+  (if (zero? n) (k 1)
+      (recur (dec n) (fn[a] (k (* n a))))))
+
+

+ + +

+ clojure +It is fairly clumsy to use this version of‘ FACT in ALGOL; it is necessary to +do something like this: +

+ + + +

begin
+integer ann
+procedure :emp(x); value 2:; integer x;
+ans :1 x;
+]'act(3. temp);
+comment Now the variable "am" has 6;
+end;
+
+

+ + + + +

+Procedure fact does not return a value, nor does temp; we must use a side +effect to get the answer out. +

+FACT is somewhat easier to use in SCHEME. We can call it like an +ordinary function in SCHEME if we supply an identity function as the second +argument. For example, (FACT 3 (LAMBDA (X) X)) returns 6. Alternatively, we +could write (FACT 3 (LAHBDA (X) (PRINT X))); we do not care about the value +of this, but about what gets printed. A third way to use the value is to +write +(FACT 3 (LAMBDA (x) (SQRT x))) +instead of +(SQRT (FACT 3 (LAMBDA (x) x))) +

+In either of these cases we care about the value of the continuation given to +FACT. Thus we care about the value of FACT if and only if we care about the +value of its continuation! +

+We can redefine other functions to take continuations in the same way. +For example, suppose we had arithmetic primitives which took continuations; to +prevent confusion, call the version of "+" which takes a continuation '++", +etc. Instead of writing +(- (+ B Z)(* 4 A C)) +we can write +

+ + +

(in-ns 'categorical.imperative)
+(defn enable-continuation "Makes f take a continuation as an additional argument."[f]
+  (fn[& args] ((fn[k] (k (apply f (reverse (rest (reverse args)))))) (last args)) ))
+(def +& (enable-continuation +))
+(def *& (enable-continuation *))
+(def -& (enable-continuation -))
+
+
+(defn quadratic[a b c k]
+(*& b b
+    (fn[x] (*& 4 a c
+               (fn[y]
+                 (-& x y k))))))
+
+

+ + + + +

+where k is the continuation for the entire expression. +

+This is an obscure way to write an algebraic expression, and we would +not advise writing code this way in general, but continuation-passing brings +out certain important features of the computation: +

¹ The operations to be performed appear in the order in which they are +

+ +

performed. In fact, they must be performed in this +order. Continuation- passing removes the need for the rule about +left-to-right argument evaluation{Note Evalorder) +

² In the usual applicative expression there are two implicit + temporary values: those of (* B B) and (* 4 A C). The first of + these values must be preserved over the computation of the second, + whereas the second is used as soon as it is produced. These facts + are manifest in the appearance and use of the variable X and Y in + the continuation-passing version. +

+ + +

+In short, the continuation-passing version specifies exactly and + explicitly what steps are necessary to compute the value of the + expression. One can think of conventional functional application + for value as being an abbreviation for the more explicit + continuation-passing style. Alternatively, one can think of the + interpreter as supplying to each function an implicit default + continuation of one argument. This continuation will receive the + value of the function as its argument, and then carry on the + computation. In an interpreter this implicit continuation is + represented by the control stack mechanism for function returns. + Now consider what computational steps are implied by: +

+(LAMBDA (A B C ...) (F X Y Z ...)) when we "apply" the LAMBDA + expression we have some values to apply it to; we let the names A, + B, C … refer to these values. We then determine the values of X, + Y, Z … and pass these values (along with "the buck", + i.e. control!) to the lambda expression F (F is either a lambda + expression or a name for one). Passing control to F is an + unconditional transfer. (Note Jrsthack) {Note Hewitthack) Note that + we want values from X, Y, Z, … If these are simple expressions, + such as variables, constants, or LAMBDA expressions, the evaluation + process is trivial, in that no temporary storage is required. In + pure continuation-passing style, all evaluations are trivial: no + combination is nested within another, and therefore no ‘hidden + temporaries" are required. But if X is a combination, say (G P Q), + then we want to think of G as a function, because we want a value + from it, and we will need an implicit temporary place to keep the + result while evaluating Y and Z. (An interpreter usually keeps these + temporary places in the control stack!) On the other hand, we do not + necessarily need a value from F. This is what we mean by tail- + recursion: F is called tail-recursively, whereas G is not. A better + name for tail-recursion would be "tail-transfer", since no real + recursion is implied. This is why we have made such a fuss about + tail-recursion: it can be used for transfer of control without + making any commitment about whether the expression expected to + return a value. Thus it can be used to model statement-like control + structures. Put another way, tail—recursion does not require a + control stack as nested recursion does. In our models of iteration + and imperative style all the LAMBDA expressions used for control (to + simulate GO statements, for example) are called in tail-recursive + fashion. +

+ +

4.2 Escape Expressions

+ +

Reynolds [Reynolds 72] defines the construction += escape x in r = +

+to evaluate the expression $r$ in an environment such that the variable $x$ is +bound to an escape function. If the escape function is never applied, then +the value of the escape expression is the value of $r$. If the escape function +is applied to an argument $a$, however, then evaluation of $r$ is aborted and the +escape expression returns $a$. {Note J-operator} (Reynolds points out that +this definition is not quite accurate, since the escape function may be called +even after the escape expression has returned a value; if this happens, it +“returns again"!) +

+As an example of the use of an escape expression, consider this +procedure to compute the harmonic mean of an array of numbers. If any of the +numbers is zero, we want the answer to be zero. We have a function hannaunl +which will sum the reciprocals of numbers in an array, or call an escape +function with zero if any of the numbers is zero. (The implementation shown +here is awkward because ALGOL requires that a function return its value by +assignment.) +

+ + +

begin
+real procedure Imrmsum(a, n. escfun)§ p
+real array at integer n; real procedure esciun(real);
+begin
+real sum;
+sum := 0;
+for i :1 0 until n-l do
+begin
+if a[i]=0 then cscfun(0);
+sum :1 sum ¢ I/a[i];
+enm
+harmsum SI sum;
+end: .
+real array b[0:99]:
+print(escape x in I00/hm-m.mm(b, 100, x));
+end
+
+

+ + + + +

+If harmsum exits normally, the number of elements is divided by the sum and +printed. Otherwise, zero is returned from the escape expression and printed +without the division ever occurring. +This program can be written in SCHEME using the built-in escape +operator CATCH: +

+ + + +

abc
+
+

+ + + +

Footnotes:

¹ DEFINITION NOT FOUND: 1 +

+ +

² DEFINITION NOT FOUND: 2 +

+ +

Date: 2011-07-15 02:49:04 EDT

Author: Guy Lewis Steele Jr. and Gerald Jay Sussman

Org version 7.6 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e categorical/imperative.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/imperative.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,772 @@ +#+TITLE: LAMBDA: The Ultimate Imperative +#+AUTHOR: Guy Lewis Steele Jr. and Gerald Jay Sussman + +* Abstract + +We demonstrate how to model the following common programming constructs in +terms of an applicative order language similar to LISP: +- Simple Recursion +- Iteration +- Compound Statements and Expressions +- GO TO and Assignment +- Continuation—Passing +- Escape Expressions +- Fluid Variables +- Call by Name, Call by Need, and Call by Reference +The models require only (possibly self-referent) lambda application, +conditionals, and (rarely) assignment. No complex data structures such as +stacks are used. The models are transparent, involving only local syntactic +transformations. +Some of these models, such as those for GO TO and assignment, are already well +known, and appear in the work of Landin, Reynolds, and others. The models for +escape expressions, fluid variables, and call by need with side effects are +new. This paper is partly tutorial in intent, gathering all the models +together for purposes of context. +This report describes research done at the Artificial Intelligence Laboratory +of the Massachusetts Institute of Teehnology. Support for the laboratory's +artificial intelligence research is provided in part by the Advanced Research +Projects Agency of the Department of Defense under Office of Naval Research +contract N000l4-75-C-0643. + +* Introduction + + We catalogue a number of common programming constructs. For each +construct we examine "typical" usage in well-known programming languages, and +then capture the essence of the semantics of the construct in terms of a +common meta-language. +The lambda calculus {Note Alonzowins} is often used as such a meta- +language. Lambda calculus offers clean semantics, but it is clumsy because it +was designed to be a minimal language rather than a convenient one. All +lambda calculus "functions" must take exactly one "argument"; the only "data +type" is lambda expressions; and the only "primitive operation‘ is variable‘ +substitution. While its utter simplicity makes lambda calculus ideal for +logicians, it is too primitive for use by programmers. The meta-language we +use is a programming language called SCHEME {Note Schemepaper) which is based +on lambda calculus. +SCHEME is a dialect of LISP. [McCarthy 62] It is an expression- +oriented, applicative order, interpreter-based language which allows one to +manipulate programs as data. It differs from most current dialects of LISP in +that it closes all lambda expressions in the environment of their definition +or declaration, rather than in the execution environment. {Note Closures} +This preserves the substitution semantics of lambda calculus, and has the +consequence that all variables are lexically scoped, as in ALGOL. [Naur 63] +Another difference is that SCHEME is implemented in such a way that tail- +recursions execute without net growth of the interpreter stack. {Note +Schemenote} We have chosen to use LISP syntax rather than, say, ALGOL syntax +because we want to treat programs as data for the purpose of describing +transformations on the code. LISP supplies names for the parts of an +executable expression and standard operators for constructing expressions and +extracting their components. The use of LISP syntax makes the structure of +such expressions manifest. We use ALGOL as an expository language, because it +is familiar to many people, but ALGOL is not sufficiently powerful to express +the necessary concepts; in particular, it does not allow functions to return +functions as values. we are thus forced to use a dialect of LISP in many +cases. +We will consider various complex programming language constructs and +show how to model them in terms of only a few simple ones. As far as possible +we will use only three control constructs from SCHEME: LAMBDA expressions, as +in LISP, which are Just functions with lexically scoped free variables; +LABELS, which allows declaration of mutually recursive procedures (Note +Labelsdef}; and IF, a primitive conditional expression. For more complex +modelling we will introduce an assignment primitive (ASET).i We will freely +assume the existence of other comon primitives, such as arithmetic functions. +The constructs we will examine are divided into four broad classes. +The first is sﬁnph?Lo0pl; this contains simple recursions and iterations, and +an introduction to the notion of continuations. The second is hmponﬂivo +Connrucls; this includes compound statements, G0 T0, and simple variable +assignments. The third is continuations, which encompasses the distinction between statements and expressions, escape operators (such as Landin's J- +operator [Landin 65] and Reynold's escape expression [Reynolds 72]). and fluid +(dynamically bound) variables. The fourth is Parameter Passing Mechanisms, such +as ALGOL call—by-name and FORTRAN call-by-location. +Some of the models presented here are already well-known, particularly +those for G0 T0 and assignment. [McCarthy 60] [Landin 65] [Reynolds 72] +Those for escape operators, fluid variables, and call-by-need with side +effects are new. +** Simple Loops +By /simple loops/ we mean constructs which enable programs to execute the +same piece of code repeatedly in a controlled manner. Variables may be made +to take on different values during each repetition, and the number of +repetitions may depend on data given to the program. +*** Simple Recursion +One of the easiest ways to produce a looping control structure is to +use a recursive function, one which calls itself to perform a subcomputation. +For example, the familiar factorial function may be written recursively in +ALGOL: ' +#+begin_src algol +integer procedure fact(n) value n: integer n: +fact := if n=0 then 1 else n*fact(n-1); +#+end_src + +The invocation =fact(n)= computes the product of the integers from 1 to n using +the identity n!=n(n-1)! (n>0). If $n$ is zero, 1 is returned; otherwise =fact=: +calls itself recursively to compute $(n-1)!$ , then multiplies the result by $n$ +and returns it. + +This same function may be written in Clojure as follows: +#+begin_src clojure +(ns categorical.imperative) +(defn fact[n] + (if (= n 0) 1 (* n (fact (dec n)))) +) +#+end_src + +#+results: +: #'categorical.imperative/fact + +Clojure does not require an assignment to the ``variable'' fan: to return a value +as ALGOL does. The IF primitive is the ALGOL if-then-else rendered in LISP +syntax. Note that the arithmetic primitives are prefix operators in SCHEME. + +*** Iteration +There are many other ways to compute factorial. One important way is +through the use of /iteration/. +A comon iterative construct is the DO loop. The most general form we +have seen in any programming language is the MacLISP DO [Moon 74]. It permits +the simultaneous initialization of any number of control variables and the +simultaneous stepping of these variables by arbitrary functions at each +iteration step. The loop is terminated by an arbitrary predicate, and an +arbitrary value may be returned. The DO loop may have a body, a series of +expressions executed for effect on each iteration. A version of the MacLISP +DO construct has been adopted in SCHEME. + +The general form of a SCHEME DO is: +#+begin_src nothing +(DO (( (init1> ) +((var2> ) +( (stepn>)) +( (value>) +(optional body>) +#+end_src +The semantics of this are that the variables are bound and initialized to the +values of the expressions, which must all be evaluated in the +environment outside the D0; then the predicate is evaluated in the new +environment, and if TRUE, the (value) is evaluated and returned. Otherwise +the (optional body) is evaluated, then each of the steppers is +evaluated in the current environment, all the variables made to have the +results as their values, the predicate evaluated again, and so on. +Using D0 loops in both ALGOL and SCHEME, we may express FACT by means +of iteration. +#+begin_src algol +integer procedure fact(n); value n; integer n: +begin +integer m, ans; +ans := 1; +for m := n step -l until 0 do ans := m*ans; +fact := ans; +end; +#+end_src + +#+begin_src clojure +(in-ns 'categorical.imperative) +(defn fact-do[n] +) +#+end_src + +Note that the SCHEME D0 loop in FACT has no body -- the stepping functions do +all the work. The ALGOL DO loop has an assignment in its body; because an +ALGOL DO loop can step only one variable, we need the assignment to step the +the variable "manually'. +In reality the SCHEME DO construct is not a primitive; it is a macro +which expands into a function which performs the iteration by tail—recursion. +Consider the following definition of FACT in SCHEME. Although it appears to +be recursive, since it "calls itself", it is entirely equivalent to the DO +loop given above, for it is the code that the DO macro expands into! It +captures the essence of our intuitive notion of iteration, because execution +of this program will not produce internal structures (e.g. stacks or variable +bindings) which increase in size with the number of iteration steps. + +#+begin_src clojure +(in-ns 'categorical.imperative) +(defn fact-do-expand[n] + (let [fact1 + (fn[m ans] + (if (zero? m) ans + (recur (dec m) (* m ans))))] + (fact1 n 1))) +#+end_src + +From this we can infer a general way to express iterations in SCHEME in +a manner isomorphic to the HacLISP D0. The expansion of the general D0 loop +#+begin_src nothing +(DO (( ) +( (init2) ) +... +( )) +( ) +) +#+end_src +is this: +#+begin_src nothing +(let [doloop +(fn [dummy ... ) +(if +(recur ... )))] + +(doloop nil ... )) +#+end_src +The identifiers =doloop= and =dummy= are chosen so as not to conflict with any +other identifiers in the program. +Note that, unlike most implementations of D0, there are no side effects +in the steppings of the iteration variables. D0 loops are usually modelled +using assignment statements. For example: +#+begin_src nothing +for r :1 0 step b until c do ; +#+end_src + +can be modelled as follows: [Naur 63] +#+begin_src nothing +begin +x := a; +L: if (x-c)*sign(n) > 0 then go to Endloop; +; +x := x+b; +go to L: +Endloop: +end; +#+end_src +Later we will see how such assignment statements can in general be +modelled without using side effects. + +* Imperative Programming +Lambda calculus (and related languages, such as ``pure LISP'') is often +used for modelling the applicative constructs of programming languages. +However, they are generally thought of as inappropriate for modelling +imperative constructs. In this section we show how imperative constructs may +be modelled by applicative SCHEME constructs. +** Compound Statements +The simplest kind of imperative construct is the statement sequencer, +for example the compound statement in ALGOL: +#+begin_src algol +begin +S1; +S2; +end +#+end_src + +This construct has two interesting properties: +- (1) It performs statement S1 before S2, and so may be used for + sequencing. +- (2) S1 is useful only for its side effects. (In ALGOL, S2 must also + be a statement, and so is also useful only for side effects, but + other languages have compound expressions containing a statement + followed by an expression.) + +The ALGOL compound statement may actually contain any number of statements, +but such statements can be expressed as a series of nested two-statement +compounds. That is: +#+begin_src algol +begin +S1; +S2; +... +Sn-1; +Sn; +end +#+end_src +is equivalent to: + +#+begin_src algol +begin +S1; +begin +S2; +begin +... + +begin +Sn-1; +Sn; +end; +end; +end: +end +#+end_src + +It is not immediately apparent that this sequencing can be expressed in a +purely applicative language. We can, however, take advantage of the implicit +sequencing of applicative order evaluation. Thus, for example, we may write a +two-statement sequence as follows: + +#+begin_src clojure +((fn[dummy] S2) S1) +#+end_src + +where =dummy= is an identifier not used in S2. From this it is +manifest that the value of S1 is ignored, and so is useful only for +side effects. (Note that we did not claim that S1 is expressed in a +purely applicative language, but only that the sequencing can be so +expressed.) From now on we will use the form =(BLOCK S1 S2)= as an +abbreviation for this expression, and =(BLOCK S1 S2...Sn-1 Sn)= as an +abbreviation for + +#+begin_src algol +(BLOCK S1 (BLOCK S2 (BLOCK ... (BLOCK Sn-1 Sn)...))) +#+end_src + +** 2.2. The G0 TO Statement + +A more general imperative structure is the compound statement with +labels and G0 T0s. Consider the following code fragment due to +Jacopini, taken from Knuth: [Knuth 74] +#+begin_src algol +begin +L1: if B1 then go to L2; S1; + if B2 then go to L2; S2; + go to L1; +L2: S3; +#+end_src + +It is perhaps surprising that this piece of code can be /syntactically/ +transformed into a purely applicative style. For example, in SCHEME we +could write: + +#+begin_src clojure +(in-ns 'categorical.imperative) +(let + [L1 (fn[] + (if B1 (L2) (do S1 + (if B2 (L2) (do S2 (L1)))))) + L2 (fn[] S3)] + (L1)) +#+end_src + +As with the D0 loop, this transformation depends critically on SCHEME's +treatment of tail-recursion and on lexical scoping of variables. The labels +are names of functions of no arguments. In order to ‘go to" the labeled code, +we merely call the function named by that label. + +** 2.3. Simple Assignment +Of course, all this sequencing of statements is useless unless the +statements have side effects. An important side effect is assignment. For +example, one often uses assignment to place intermediate results in a named +location (i.e. a variable) so that they may be used more than once later +without recomputing them: + +#+begin_src algol +begin +real a2, sqrtdisc; +a2 := 2*a; +sqrtdisc := sqrt(b^2 - 4*a*c); +root1 := (- b + sqrtdisc) / a2; +root2 := (- b - sqrtdisc) / a2; +print(root1); +print(root2); +print(root1 + root2); +end +#+end_src + +It is well known that such naming of intermediate results may be accomplished +by calling a function and binding the formal parameter variables to the +results: + +#+begin_src clojure +(fn [a2 sqrtdisc] + ((fn[root1 root2] + (do (println root1) + (println root2) + (println (+ root1 root2)))) + (/ (+ (- b) sqrtdisc) a2) + (/ (- (- b) sqrtdisc) a2)) + + (* 2 a) + (sqrt (- (* b b) (* 4 a c)))) +#+end_src +This technique can be extended to handle all simple variable assignments which +appear as statements in blocks, even if arbitrary G0 T0 statements also appear +in such blocks. {Note Mccarthywins} + +For example, here is a program which uses G0 TO statements in the form +presented before; it determines the parity of a non-negative integer by +counting it down until it reaches zero. + +#+begin_src algol +begin +L1: if a = 0 then begin parity := 0; go to L2; end; + a := a - 1; + if a = 0 then begin parity := 1; go to L2; end; + a := a - 1; + go to L1; +L2: print(parity); +#+end_src + +This can be expressed in SCHEME: + +#+begin_src clojure +(let + [L1 (fn [a parity] + (if (zero? a) (L2 a 0) + (L3 (dec a) parity))) + L3 (fn [a parity] + (if (zero? a) (L2 a 1) + (L1 (dec a) parity))) + L2 (fn [a parity] + (println parity))] + (L1 a parity)) +#+end_src + +The trick is to pass the set of variables which may be altered as arguments to +the label functions. {Note Flowgraph} It may be necessary to introduce new +labels (such as L3) so that an assignment may be transformed into the binding +for a function call. At worst, one may need as many labels as there are +statements (not counting the eliminated assignment and GO TO statements). + +** Compound Expressions ' +At this point we are almost in a position to model the most general +form of compound statement. In LISP, this is called the 'PROG feature". In +addition to statement sequencing and G0 T0 statements, one can return a /value/ +from a PROG by using the RETURN statement. + +Let us first consider the simplest compound statement, which in SCHEME +we call BLOCK. Recall that +=(BLOCK s1 sz)= is defined to be =((lambda (dummy) s2) s1)= + +Notice that this not only performs Sl before S2, but also returns the value of +S2. Furthermore, we defined =(BLOCK S1 S2 ... Sn)= so that it returns the value +of =Sn=. Thus BLOCK may be used as a compound expression, and models a LISP +PROGN, which is a PROG with no G0 T0 statements and an implicit RETURN of the +last ``statement'' (really an expression). + +Host LISP compilers compile D0 expressions by macro-expansion. We have +already seen one way to do this in SCHEME using only variable binding. A more +common technique is to expand the D0 into a PROG, using variable assignments +instead of bindings. Thus the iterative factorial program: + +#+begin_src clojure +(oarxnz FACT . +(LAMaoA (n) . +(D0 ((M N (- H 1)) +(Ans 1 (* M Ans))) +((- H 0) A"$)))) +#+end_src + +would expand into: + +#+begin_src clojure +(DEFINE FACT +. (LAMeoA (M) +(PRO6 (M Ans) +(sszro M n +Ans 1) ~ +LP (tr (- M 0) (RETURN Ans)) +(ssero M (- n 1) +Ans (' M Ans)) +(60 LP)))) +#+end_src + +where SSETQ is a simultaneous multiple assignment operator. (SSETQ is not a +SCHEME (or LISP) primitive. It can be defined in terms of a single assignment +operator, but we are more interested here in RETURN than in simultaneous +assignment. The SSETQ's will all be removed anyway and modelled by lambda +binding.) We can apply the same technique we used before to eliminate G0 T0 +statements and assignments from compound statements: + +#+begin_src clojure +(DEFINE FACT +(LAHBOA (I) +(LABELS ((L1 (LAManA (M Ans) +(LP n 1))) +(LP (LAMaoA (M Ans) +(IF (- M o) (nztunn Ans) +(£2 H An$)))) +(L2 (LAMaoA (M An ) +(LP (- " 1) (' H ﬂN$))))) +(L1 NIL NlL)))) +#+end_src clojure + +We still haven't done anything about RETURN. Let's see... +- ==> the value of (FACT 0) is the value of (Ll NIL NIL) +- ==> which is the value of (LP 0 1) +- ==> which is the value of (IF (= 0 0) (RETURN 1) (L2 0 1)) +- ==> which is the value of (RETURN 1) + +Notice that if RETURN were the /identity +function/ (LAMBDA (X) X), we would get the right answer. This is in fact a +general truth: if we Just replace a call to RETURN with its argument, then +our old transformation on compound statements extends to general compound +expressions, i.e. PROG. + +* Continuations +Up to now we have thought of SCHEME's LAMBDA expressions as functions, +and of a combination such as =(G (F X Y))= as meaning ``apply the function F to +the values of X and Y, and return a value so that the function G can be +applied and return a value ...'' But notice that we have seldom used LAMBDA +expressions as functions. Rather, we have used them as control structures and +environment modifiers. For example, consider the expression: +=(BLOCK (PRINT 3) (PRINT 4))= + +This is defined to be an abbreviation for: +=((LAMBDA (DUMMY) (PRINT 4)) (PRINT 3))= + +We do not care about the value of this BLOCK expression; it follows that we +do not care about the value of the =(LAMBDA (DUMMY) ...)=. We are not using +LAMBDA as a function at all. + +It is possible to write useful programs in terms of LAHBDA expressions +in which we never care about the value of /any/ lambda expression. We have +already demonstrated how one could represent any "FORTRAN" program in these +terms: all one needs is PROG (with G0 and SETQ), and PRINT to get the answers +out. The ultimate generalization of this imperative programing style is +/continuation-passing/. (Note Churchwins} . + +** Continuation-Passing Recursion +Consider the following alternative definition of FACT. It has an extra +argument, the continuation, which is a function to call with the answer, when +we have it, rather than return a value + +#+begin_src algol. +procedure Inc!(n, c); value n, c; +integer n: procedure c(integer value); +if n-=0 then c(l) else +begin +procedure l(!mp(d) value a: integer a; +c(mm); +Iacl(n-1, romp): +end; +#+end_src + +#+begin_src clojure +(in-ns 'categorical.imperative) +(defn fact-cps[n k] + (if (zero? n) (k 1) + (recur (dec n) (fn[a] (k (* n a)))))) +#+end_src clojure +It is fairly clumsy to use this version of‘ FACT in ALGOL; it is necessary to +do something like this: + +#+begin_src algol +begin +integer ann +procedure :emp(x); value 2:; integer x; +ans :1 x; +]'act(3. temp); +comment Now the variable "am" has 6; +end; +#+end_src + +Procedure =fact= does not return a value, nor does =temp=; we must use a side +effect to get the answer out. + +=FACT= is somewhat easier to use in SCHEME. We can call it like an +ordinary function in SCHEME if we supply an identity function as the second +argument. For example, =(FACT 3 (LAMBDA (X) X))= returns 6. Alternatively, we +could write =(FACT 3 (LAHBDA (X) (PRINT X)))=; we do not care about the value +of this, but about what gets printed. A third way to use the value is to +write +=(FACT 3 (LAMBDA (x) (SQRT x)))= +instead of +=(SQRT (FACT 3 (LAMBDA (x) x)))= + +In either of these cases we care about the value of the continuation given to +FACT. Thus we care about the value of FACT if and only if we care about the +value of its continuation! + +We can redefine other functions to take continuations in the same way. +For example, suppose we had arithmetic primitives which took continuations; to +prevent confusion, call the version of "+" which takes a continuation '++", +etc. Instead of writing +=(- (+ B Z)(* 4 A C))= +we can write +#+begin_src clojure +(in-ns 'categorical.imperative) +(defn enable-continuation "Makes f take a continuation as an additional argument."[f] + (fn[& args] ((fn[k] (k (apply f (reverse (rest (reverse args)))))) (last args)) )) +(def +& (enable-continuation +)) +(def *& (enable-continuation *)) +(def -& (enable-continuation -)) + + +(defn quadratic[a b c k] +(*& b b + (fn[x] (*& 4 a c + (fn[y] + (-& x y k)))))) +#+end_src + +where =k= is the continuation for the entire expression. + +This is an obscure way to write an algebraic expression, and we would +not advise writing code this way in general, but continuation-passing brings +out certain important features of the computation: + - [1] The operations to be performed appear in the order in which they are +performed. In fact, they /must/ be performed in this +order. Continuation- passing removes the need for the rule about +left-to-right argument evaluation{Note Evalorder) +- [2] In the usual applicative expression there are two implicit + temporary values: those of =(* B B)= and =(* 4 A C)=. The first of + these values must be preserved over the computation of the second, + whereas the second is used as soon as it is produced. These facts + are /manifest/ in the appearance and use of the variable X and Y in + the continuation-passing version. + +In short, the continuation-passing version specifies /exactly/ and + explicitly what steps are necessary to compute the value of the + expression. One can think of conventional functional application + for value as being an abbreviation for the more explicit + continuation-passing style. Alternatively, one can think of the + interpreter as supplying to each function an implicit default + continuation of one argument. This continuation will receive the + value of the function as its argument, and then carry on the + computation. In an interpreter this implicit continuation is + represented by the control stack mechanism for function returns. + Now consider what computational steps are implied by: + +=(LAMBDA (A B C ...) (F X Y Z ...))= when we "apply" the LAMBDA + expression we have some values to apply it to; we let the names A, + B, C ... refer to these values. We then determine the values of X, + Y, Z ... and pass these values (along with "the buck", + i.e. control!) to the lambda expression F (F is either a lambda + expression or a name for one). Passing control to F is an + unconditional transfer. (Note Jrsthack) {Note Hewitthack) Note that + we want values from X, Y, Z, ... If these are simple expressions, + such as variables, constants, or LAMBDA expressions, the evaluation + process is trivial, in that no temporary storage is required. In + pure continuation-passing style, all evaluations are trivial: no + combination is nested within another, and therefore no ‘hidden + temporaries" are required. But if X is a combination, say (G P Q), + then we want to think of G as a function, because we want a value + from it, and we will need an implicit temporary place to keep the + result while evaluating Y and Z. (An interpreter usually keeps these + temporary places in the control stack!) On the other hand, we do not + necessarily need a value from F. This is what we mean by tail- + recursion: F is called tail-recursively, whereas G is not. A better + name for tail-recursion would be "tail-transfer", since no real + recursion is implied. This is why we have made such a fuss about + tail-recursion: it can be used for transfer of control without + making any commitment about whether the expression expected to + return a value. Thus it can be used to model statement-like control + structures. Put another way, tail—recursion does not require a + control stack as nested recursion does. In our models of iteration + and imperative style all the LAMBDA expressions used for control (to + simulate GO statements, for example) are called in tail-recursive + fashion. + +** Escape Expressions +Reynolds [Reynolds 72] defines the construction += escape x in r = + +to evaluate the expression $r$ in an environment such that the variable $x$ is +bound to an escape function. If the escape function is never applied, then +the value of the escape expression is the value of $r$. If the escape function +is applied to an argument $a$, however, then evaluation of $r$ is aborted and the +escape expression returns $a$. {Note J-operator} (Reynolds points out that +this definition is not quite accurate, since the escape function may be called +even after the escape expression has returned a value; if this happens, it +“returns again"!) + +As an example of the use of an escape expression, consider this +procedure to compute the harmonic mean of an array of numbers. If any of the +numbers is zero, we want the answer to be zero. We have a function hannaunl +which will sum the reciprocals of numbers in an array, or call an escape +function with zero if any of the numbers is zero. (The implementation shown +here is awkward because ALGOL requires that a function return its value by +assignment.) +#+begin_src algol +begin +real procedure Imrmsum(a, n. escfun)§ p +real array at integer n; real procedure esciun(real); +begin +real sum; +sum := 0; +for i :1 0 until n-l do +begin +if a[i]=0 then cscfun(0); +sum :1 sum ¢ I/a[i]; +enm +harmsum SI sum; +end: . +real array b[0:99]: +print(escape x in I00/hm-m.mm(b, 100, x)); +end +#+end_src + +If harmsum exits normally, the number of elements is divided by the sum and +printed. Otherwise, zero is returned from the escape expression and printed +without the division ever occurring. +This program can be written in SCHEME using the built-in escape +operator =CATCH=: + +#+begin_src clojure +(in-ns 'categorical.imperative) +(defn harmonic-sum[coll escape] + ((fn [i sum] + (cond (= i (count coll) ) sum + (zero? (nth coll i)) (escape 0) + :else (recur (inc i) (+ sum (/ 1 (nth coll i)))))) + 0 0)) + +#+end_src + +This actually works, but elucidates very little about the nature of ESCAPE. +We can eliminate the use of CATCH by using continuation-passing. Let us do +for HARMSUM what we did earlier for FACT. Let it take an extra argument C, +which is called as a function on the result. + +#+begin_src clojure +(in-ns 'categorical.imperative) +(defn harmonic-sum-escape[coll escape k] + ((fn[i sum] + (cond (= i (count coll)) (k sum) + (zero? (nth coll i)) (escape 0) + (recur (inc i) (+ sum (/ 1 (nth coll i)))))) + 0 0)) + +(let [B (range 0 100) + after-the-catch println] + (harmonic-sum-escape B after-the-catch (fn[y] (after-the-catch (/ 100 y))))) + +#+end_src + +Notice that if we use ESCFUN, then C does not get called. In this way the +division is avoided. This example shows how ESCFUN may be considered to be an +"alternate continuation". + +** Dynamic Variable Scoping +In this section we will consider the problem of dynamically scoped, or +"fluid", variables. These do not exist in ALGOL, but are typical of many LISP +implementations, ECL, and APL. He will see that fluid variables may be +modelled in more than one way, and that one of these is closely related to +continuation—pass1ng. + +*** Free (Global) Variables +Consider the following program to compute square roots: + +#+begin_src clojure +(defn sqrt[x tolerance] + ( +(DEFINE soar +(LAHBDA (x EPSILON) +(Pace (ANS) +(stro ANS 1.0) +A (coup ((< (ADS (-s x (~s ANS ANS))) EPSILON) +(nzruau ANS))) ' +(sero ANS (//s (+5 x (//s x ANS)) 2.0)) +(60 A)))) . +This function takes two arguments: the radicand and the numerical tolerance +for the approximation. Now suppose we want to write a program QUAD to compute +solutions to a quadratic equation; p +(perms QUAD +(LAHDDA (A a c) +((LAHBDA (A2 sonmsc) ' +(LIST (/ (+ (- a) SQRTDISC) AZ) +(/ (- (- B) SORTOISC) AZ))) +(' 2 A) +- (SORT (- (9 D 2) (' 4 A C)) (tolerance>)))) +#+end_src +It is not clear what to write for (tolerance). One alternative is to pick +some tolerance for use by QUAD and write it as a constant in the code. This +is undesirable because it makes QUAD inflexible and hard to change. _Another +is to make QUAD take an extra argument and pass it to SQRT: +(DEFINE QUAD +(LANODA (A 8 C EPSILON) +(soar ... EPSILON) ...)) +This is undesirable because EPSILDN is not really part of the problem QUAD is +supposed to solve, and we don't want the user to have to provide it. +Furthermore, if QUAD were built into some larger function, and that into +another, all these functions would have to pass EPSILON down as an extra +argument. A third.possibi1ity would be to pass the SQRT function as an +argument to QUAD (don't laugh!), the theory being to bind EPSILON at the +appropriate level like this: A U‘ A +(QUAD 3 A 5 (LAMBDA (X) (SORT X ))) +where we define QUAD as: +(DEFINE QUAD +(LAMBDA (A a c soar) ...)) diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_0dd51e4ab8cb4eb616439d758fbec4ed120e8a15.png Binary file categorical/ltxpng/plausible_0dd51e4ab8cb4eb616439d758fbec4ed120e8a15.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_327b1246bf52889b71445e990ea38a8414f683cc.png Binary file categorical/ltxpng/plausible_327b1246bf52889b71445e990ea38a8414f683cc.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_34a4f0f7fd91bf7212f437fd9e0fae7c04dbdae1.png Binary file categorical/ltxpng/plausible_34a4f0f7fd91bf7212f437fd9e0fae7c04dbdae1.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_5072ddd55abeb45994e1f467e2041032d8a90350.png Binary file categorical/ltxpng/plausible_5072ddd55abeb45994e1f467e2041032d8a90350.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_5434f18d3d7ac2858cd4f914be8f5d2b49c79a90.png Binary file categorical/ltxpng/plausible_5434f18d3d7ac2858cd4f914be8f5d2b49c79a90.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_60c89d70f35bbeef7c8c5572941c3398e65f9eb1.png Binary file categorical/ltxpng/plausible_60c89d70f35bbeef7c8c5572941c3398e65f9eb1.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_671e7f9a1d6c0ed854ad469bb2e6446cf536f6ce.png Binary file categorical/ltxpng/plausible_671e7f9a1d6c0ed854ad469bb2e6446cf536f6ce.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_70e445784826d50e75fa759496de7dbd78e0c29e.png Binary file categorical/ltxpng/plausible_70e445784826d50e75fa759496de7dbd78e0c29e.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_9471153ffb8045fd30a14628197f6e1a31c4d2bf.png Binary file categorical/ltxpng/plausible_9471153ffb8045fd30a14628197f6e1a31c4d2bf.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_c0cfa0be1f5f4c93c4a07d2bf03f296a87ddb3f5.png Binary file categorical/ltxpng/plausible_c0cfa0be1f5f4c93c4a07d2bf03f296a87ddb3f5.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_d282b164f31f8cb4d3238ea37d09715da89272aa.png Binary file categorical/ltxpng/plausible_d282b164f31f8cb4d3238ea37d09715da89272aa.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_dd81a36f0fc30866fc59a06fba09770696fbcfce.png Binary file categorical/ltxpng/plausible_dd81a36f0fc30866fc59a06fba09770696fbcfce.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/ltxpng/plausible_df22ff07c07c3eae08269551b79e75e13ed321f9.png Binary file categorical/ltxpng/plausible_df22ff07c07c3eae08269551b79e75e13ed321f9.png has changed diff -r 8d8278e09888 -r b4de894a1e2e categorical/monad.clj --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/monad.clj Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,14 @@ + +(ns categorical.monad) +(use 'clojure.contrib.monads) + +(in-ns 'categorical.monad) + +;; To implement nondeterministic programs, we'll use a lazy seq to represent a value which may be any one of the members of seq. + +(defmonad amb + [ + m-result (fn[& vals] (cons 'amb vals)) + m-bind (fn[amb-seq f] (cons 'amb (map f (rest amb-seq)))) + ] +) diff -r 8d8278e09888 -r b4de894a1e2e categorical/monad.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/monad.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,251 @@ + + + + +Monads in CS and Math + + + + + + + + + + + +

+ + + +

+Monads in CS and Math +

+ +

+The purpose of a monad is to introduce a new data type into your +language and to integrate it with existing types; a monad augments a +pre-existing data type $A$, creating a new type $B$. +

+When describing monads, I will use the following terminology: +

Values of type $A$ are called ordinary values, because they belong + to the pre-existing data type. +
Values of type $B$ are called monadic values, because they + belong to the type introduced by the monad. +
A function $A\rightarrow B$ is called a monadic function; such + functions accept regular values as input and return monadic values. +

+ + +

+A monad consists of +two components for achieving this purpose: +

A function that converts ordinary values into monadic values: +
+ First, in order to create a new data type, each monad has an + upgrade function which systematically adds structure to values of + the existing data type, $A$. Because the upgrade function produces + enhanced values in a systematic way, the enhanced values constitute + a sort of new data type, $B$. +

+ + +

+ + +$\text{upgrade}:A\rightarrow B$ + +

+ + +

A function that enables monadic functions to accept monadic values: Second, each monad has a pipe function which is a + protocol for combining a monadic function and a monadic value to + produce a monadic value. The pipe function facilitates composition: + notice that a collection of monadic functions $A\rightarrow B$ is composable only if they have been modified by the + pipe function to become functions of type $B\rightarrow B$, otherwise their input and output types do not match. +

+ +

+ + +$\text{pipe}:B\times B^A\rightarrow B$ + +

+ + + + + + +

(ns categorical.monad)
+
+

+ + + + + +

1 Examples of monads +
- 1.1 The nondeterministic monad
+

+ +

1 Examples of monads

+ + +

+ +

1.1 The nondeterministic monad

+ + +

+We'll implement nondeterministic programming by defining a monad amb +that transforms ordinary deterministic functions into functions that ``ambiguously'' return zero or more +values. +

+ + + +

Date: 2011-07-11 03:07:07 EDT

Author: Robert McIntyre

Org version 7.6 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e categorical/monad.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/monad.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,76 @@ +#+TITLE: Monads in CS and Math +#+BEGIN_HTML +

+{{{TITLE}}} +

+#+END_HTML + +The purpose of a /monad/ is to introduce a new data type into your +language and to integrate it with existing types; a monad augments a +pre-existing data type $A$, creating a new type $B$. + +When describing monads, I will use the following terminology: +- Values of type $A$ are called *ordinary values*, because they belong + to the pre-existing data type. +- Values of type $B$ are called *monadic values*, because they + belong to the type introduced by the monad. +- A function $A\rightarrow B$ is called a *monadic function*; such + functions accept regular values as input and return monadic values. + + + +A monad consists of +two components for achieving this purpose: + +- A function that converts ordinary values into monadic values :: + First, in order to create a new data type, each monad has an + /upgrade function/ which systematically adds structure to values of + the existing data type, $A$. Because the upgrade function produces + enhanced values in a systematic way, the enhanced values constitute + a sort of new data type, $B$. + +#+begin_quote +$\text{upgrade}:A\rightarrow B$ +#+end_quote + +- A function that enables monadic functions to accept monadic values :: Second, each monad has a /pipe function/ which is a + protocol for combining a monadic function and a monadic value to + produce a monadic value. The pipe function facilitates composition: + notice that a collection of monadic functions $A\rightarrow B$ is composable only if they have been modified by the + pipe function to become functions of type $B\rightarrow B$, otherwise their input and output types do not match. +#+begin_quote +$\text{pipe}:B\times B^A\rightarrow B$ +#+end_quote + +#+srcname: intro +#+begin_src clojure :results silent +(ns categorical.monad) +(use 'clojure.contrib.monads) +#+end_src + + +* Examples of monads +** The nondeterministic monad + +We'll implement nondeterministic programming by defining a monad =amb= +that transforms ordinary deterministic functions into functions that ``ambiguously'' return zero or more +values. + +#+srcname: stuff +#+begin_src clojure :results silent +(in-ns 'categorical.monad) + +;; To implement nondeterministic programs, we'll use a lazy seq to represent a value which may be any one of the members of seq. + +(defmonad amb + [ + m-result (fn[& vals] (cons 'amb vals)) + m-bind (fn[amb-seq f] (cons 'amb (map f (rest amb-seq)))) + ] +) +#+end_src + +#+begin_src clojure :results silent :tangle monad.clj :noweb yes +<> +<> +#+end_src diff -r 8d8278e09888 -r b4de894a1e2e categorical/morphisms.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/morphisms.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,240 @@ + + + + +Categorical Programs + + + + + + + + + + + + +

+ + + + + +

Categorical Programs

+ + +

1 The Category of Types +
- 1.1 Building new data types out of existing ones
+

+ +

1 The Category of Types

+ +

Every computer language with data types (such as Integers, +Floats, or +Strings) corresponds with a certain category; the objects of the category are the data types +of the language, and there is one arrow between data types $A$ and $B$ +for every function in the language whose argument has type $A$ and whose return value +has type $B$. +

+ + + +

(ns categorical.morphisms)
+
+

+ + + + + +

+ +

1.1 Building new data types out of existing ones

+ + +

+We will now discuss two processes which combine existing +data types to yield composite data types. These higher-order processes will be +functorial, which means that they will recombine not only the objects of +our category (the data types), but the arrows as well (functions with +typed input and output). +

+The product functor combines two data types, $A$ and $B$, producing +a new type, denoted $A\times B$. Conceptually, values with the data type $A\times +B$ have two components, the first of which is a value of type $A$, the +second of which is a value of type $B$. Concretely, the data type $A\times B$ +can be implemented as pairs of values, the first of which is +of type $A$ and the second of which is of type $B$. With this +implementation in mind, you can see that the type $A$ and type $B$ +constituents of type $A\times B$ value can be recovered; the so-called +projection functors first and second recover them. A third +functor, the diagonal functor, completes our toolkit; given +value $a$ of type $A$, it produces the value $(a,a)$ of type $A\times A$. +

+ + + +

(defn product[a b] (list a b))
+(def projections (list first second))
+(defn diagonal[a] (product a a))
+
+

+ + + + +

+The coproduct functor combines two data types, $A$ and $B$, +producing a new type, denoted $A+B$. Conceptually, you construct the +$A+B$ datatype by creating labelled versions of the data types $A$ and +$B$ — we'll call them black-$A$ and white-$B$ — and +then creating a data type $A+B$ which contains all the black-$A$ +values and all the white-$B$ values. Concretely, you can implement the +$A+B$ data type using label-value pairs: the first value in the pair +is either the label `black' or the label `white'; the second +value must have type $A$ if the label is black, or type $B$ if the +label is white. +

+Now, the product functor came with projection functors first +and second that take $A\times B$ values and return (respectively) $A$ values and +$B$ values. The coproduct functor, conversely, comes with two +injection functors make-black and make-white that perform the +opposite function: they take (respectively) $A$ values and $B$ values, +promoting them to the $A+B$ datatype by making the $A$'s into +black-$A$'s and the $B$'s into white-$B$'s. +

+ + +

Date: 2011-07-01 12:28:39 CDT

Author: Dylan Holmes

Org version 7.5 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e categorical/morphisms.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/morphisms.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,81 @@ +#+TITLE: Categorical Programs +#+AUTHOR: Dylan Holmes + +#+STYLE: + +#+BEGIN_HTML +

{{{title}}}

+#+END_HTML + +* The Category of Types +Every computer language with data types (such as Integers, +Floats, or +Strings) corresponds with a certain category; the objects of the category are the data types +of the language, and there is one arrow between data types $A$ and $B$ +for every function in the language whose argument has type $A$ and whose return value +has type $B$. + +#+begin_src clojure :results silent +(ns categorical.morphisms) +#+end_src + +** Building new data types out of existing ones + +We will now discuss two processes which combine existing +data types to yield composite data types. These higher-order processes will be +/functorial/, which means that they will recombine not only the objects of +our category (the data types), but the arrows as well (functions with +typed input and output). + +The =product= functor combines two data types, $A$ and $B$, producing +a new type, denoted $A\times B$. Conceptually, values with the data type $A\times +B$ have two components, the first of which is a value of type $A$, the +second of which is a value of type $B$. Concretely, the data type $A\times B$ +can be implemented as pairs of values, the first of which is +of type $A$ and the second of which is of type $B$. With this +implementation in mind, you can see that the type $A$ and type $B$ +constituents of type $A\times B$ value can be recovered; the so-called +projection functors =first= and =second= recover them. A third +functor, the =diagonal= functor, completes our toolkit; given +value $a$ of type $A$, it produces the value $(a,a)$ of type $A\times A$. + +#+begin_src clojure :results silent +(defn product[a b] (list a b)) +(def projections (list first second)) +(defn diagonal[a] (product a a)) +#+end_src + +The =coproduct= functor combines two data types, $A$ and $B$, +producing a new type, denoted $A+B$. Conceptually, you construct the +$A+B$ datatype by creating labelled versions of the data types $A$ and +$B$ \mdash{} we'll call them black-$A$ and white-$B$ \mdash{} and +then creating a data type $A+B$ which contains all the black-$A$ +values and all the white-$B$ values. Concretely, you can implement the +$A+B$ data type using label-value pairs: the first value in the pair +is either the label `black' or the label `white'; the second +value must have type $A$ if the label is black, or type $B$ if the +label is white. + +Now, the =product= functor came with projection functors =first= +and =second= that take $A\times B$ values and return (respectively) $A$ values and +$B$ values. The =coproduct= functor, conversely, comes with two +injection functors =make-black= and =make-white= that perform the +opposite function: they take (respectively) $A$ values and $B$ values, +promoting them to the $A+B$ datatype by making the $A$'s into +black-$A$'s and the $B$'s into white-$B$'s. + +#the =coproduct= functor comes with two injection functors =make-black= +#and =make-white= that take (respectively) $A$ values and $B$ values, +#promoting them to the $A+B$ datatype by making them into black-$A$'s +#or white-$B$'s. + + +#The =coproduct= functor combines two data types, $A$ and $B$, +#producing a new type, denoted $A+B$. Conceptually, values with the +#data type $A+B$ can be either values with type $A$ that are designated /of the first +#kind/, or values with type $B$ that are designated /of the second +#kind/. Concretely, the data type $A+B$ can be implemented using pairs +#of values in conjunction with a labelling convention; the first value +#in the pair is of either type $A$ or $B$, and the second value is a label +#indicating whether it is of the first or second kind. Our labelling +#scheme uses the keywords =:black= and =:white= for this purpose. diff -r 8d8278e09888 -r b4de894a1e2e categorical/plausible.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/plausible.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,336 @@ + + + + +Categorification of Plausible Reasoning + + + + + + + + + + + +

+ + + + +

1 Deductive and inductive posets +
- 1.1 Definition
- 1.2 A canonical map from deductive posets to inductive posets
+
2 Assigning plausibilities to inductive posets +
+

+ +

1 Deductive and inductive posets

+ + + +

+ +

1.1 Definition

+ +

If you have a collection $P$ of logical propositions, you can order them by +implication: $a$ precedes $b$ if and only if $a$ implies +$b$. This makes $P$ into a poset. Since the ordering arose from +deductive implication, we'll call this a deductive poset. +

+If you have a deductive poset $P$, you can create a related poset $P^*$ +as follows: the underlying set is the same, and for any two +propositions $a$ and $b$ in $P$, $a$ precedes +$ab$ in $P^*$. We'll call this an inductive poset. +

+ +

1.2 A canonical map from deductive posets to inductive posets

+ +

Each poset corresponds with a poset-category, that is a category with +at most one arrow between any two objects. Considered as categories, +inductive and deuctive posets are related as follows: there is a map +$\mathscr{F}$ which sends each arrow $a\rightarrow b$ in $P$ to +the arrow $a\rightarrow ab$ in $P^*$. In fact, since $a$ implies +$b$ if and only if $a = ab$, $\mathscr{F}$ sends each arrow in $P$ to +an identity arrow in $P^*$ (specifically, it sends the arrow +$a\rightarrow b$ to the identity arrow $a\rightarrow a$). +

+ +

2 Assigning plausibilities to inductive posets

+ + +

+Inductive posets encode the relative (qualitative) plausibilities of its +propositions: there exists an arrow $x\rightarrow y$ only if $x$ +is at least as plausible as $y$. +

+ +

2.1 Consistent reasoning as a commutative diagram

+ +

Inductive categories enable the following neat trick: we can interpret +the objects of $P^*$ as states of given information and interpret +each arrow $a\rightarrow ab$ in $P^*$ as an inductive inference: the arrow +$a\rightarrow ab$ represents an inferential leap from the state of +knowledge where only $a$ is given to the state of knowledge where +both $a$ and $b$ are given— in this way, it represents +the process of inferring $b$ when given $a$, and we label the +arrow with $(b|a)$. +

+This trick has several important features that suggest its usefulness, +namely +

Composition of arrows corresponds to compound inference. +
In the special case of deductive inference, the inferential arrow is an + identity; the source and destination states of knowledge are the same. +
One aspect of the consistency requirement of Jaynes¹ takes the form of a + commutative square: $x\rightarrow ax \rightarrow abx$ = + $x\rightarrow bx \rightarrow abx$ is the categorified version of + $(AB|X)=(A|X)\cdot(B|AX)=(B|X)\cdot(A|BX)$. +
We can make plausibility assignments by enriching the inductive + category $P^*$ over some monoidal category, e.g. the set of real numbers + (considered as a category) with its usual multiplication. When we do, + the identity arrows of $P^*$ —corresponding to + deductive inferences— are assigned a value of certainty automatically. +

+ + +

+ +

2.2 ``Multiplicity'' is reciprocal probability

+ +

The natural numbers have a comparatively concrete origin: they are the +result of decategorifying the category of finite sets², or the +coequalizer of the arrows from a one-object category to a two-object +category with a single nonidentity arrow. Extensions of the set of +natural numbers— such as +the set of integers or rational numbers or real numbers— strike +me as being somewhat more abstract (however, see the Eudoxus +construction of the real numbers). +

+Jaynes points out that our existing choice of scale for probabilities +(i.e., the scale from 0 for impossibility to 1 for +certainty) has a degree of freedom: any monotonic function of +probability encodes the same information that probability does. Though +the resulting laws for compound probability and so on change in form +when probabilities are changed, they do not change in content. +

+With this in mind, it seems natural and permissible to use not probability but +reciprocal probability instead. This scale, which we might + call multiplicity, ranges from 1 (certainty) to +positive infinity (impossibility); higher numbers are ascribed to +less-plausible events. +

+In this way, the ``probability'' +associated with choosing one out of $n$ indistinguishable choices +becomes identified with $n$. +

+ +

2.3 Laws for multiplicity

+ +

Jaynes derives laws of probability; either his method or his results +can be used to obtain laws for multiplicities. +

product rule: The product rule is unchanged: $\xi(AB|X)=\xi(A|X)\cdot + \xi(B|AX) = \xi(B|X)\cdot \xi(A|BX)$ +
certainty: States of absolute certainty are assigned a multiplicity + of 1. States of absolute impossibility are assigned a + multiplicity of positive infinity. +
entropy: In terms of probability, entropy has the form $S=-\sum_i + p_i \ln{p_i} = \sum_i p_i (-\ln{p_i}) = \sum_i p_i \ln{(1/p_i)} + $. Hence, in terms of multiplicity, entropy + has the form $S = \sum_i \frac{\ln{\xi_i}}{\xi_i} $. + +
+ Another interesting quantity is $\exp{S}$, which behaves + multiplicitively rather than additively. $\exp{S} = + \prod_i \exp{\frac{\ln{\xi_i}}{\xi_i}} = + \left(\exp{\ln{\xi_i}}\right)^{1/\xi_i} = \prod_i \xi_i^{1/\xi_i} $ +

+ + + +

Footnotes:

¹ (IIIa) If a conclusion can be reasoned out in more than one way, then every possible way must lead to the same result. +

+ +

² As Baez explains. +

+ +

Date: 2011-07-09 14:19:42 EDT

Author: Dylan Holmes

Org version 7.6 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e categorical/plausible.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/plausible.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,114 @@ +#+TITLE: Categorification of Plausible Reasoning +#+AUTHOR: Dylan Holmes +#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js" +* COMMENT #+OPTIONS: LaTeX:dvipng + +* Deductive and inductive posets + +** Definition +If you have a collection $P$ of logical propositions, you can order them by +implication: $a$ precedes $b$ if and only if $a$ implies +$b$. This makes $P$ into a poset. Since the ordering arose from +deductive implication, we'll call this a /deductive poset/. + +If you have a deductive poset $P$, you can create a related poset $P^*$ +as follows: the underlying set is the same, and for any two +propositions $a$ and $b$ in $P$, $a$ precedes +$ab$ in $P^*$. We'll call this an /inductive poset/. + +** A canonical map from deductive posets to inductive posets +Each poset corresponds with a poset-category, that is a category with +at most one arrow between any two objects. Considered as categories, +inductive and deuctive posets are related as follows: there is a map +$\mathscr{F}$ which sends each arrow $a\rightarrow b$ in $P$ to +the arrow $a\rightarrow ab$ in $P^*$. In fact, since $a$ implies +$b$ if and only if $a = ab$, $\mathscr{F}$ sends each arrow in $P$ to +an identity arrow in $P^*$ (specifically, it sends the arrow +$a\rightarrow b$ to the identity arrow $a\rightarrow a$). + + +** Assigning plausibilities to inductive posets + +Inductive posets encode the relative (/qualitative/) plausibilities of its +propositions: there exists an arrow $x\rightarrow y$ only if $x$ +is at least as plausible as $y$. + +*** Consistent reasoning as a commutative diagram +Inductive categories enable the following neat trick: we can interpret +the objects of $P^*$ as states of given information and interpret +each arrow $a\rightarrow ab$ in $P^*$ as an inductive inference: the arrow +$a\rightarrow ab$ represents an inferential leap from the state of +knowledge where only $a$ is given to the state of knowledge where +both $a$ and $b$ are given\mdash{} in this way, it represents +the process of inferring $b$ when given $a$, and we label the +arrow with $(b|a)$. + +This trick has several important features that suggest its usefulness, +namely + - Composition of arrows corresponds to compound inference. + - In the special case of deductive inference, the inferential arrow is an + identity; the source and destination states of knowledge are the same. + - One aspect of the consistency requirement of Jaynes[fn:1] takes the form of a + commutative square: $x\rightarrow ax \rightarrow abx$ = + $x\rightarrow bx \rightarrow abx$ is the categorified version of + $(AB|X)=(A|X)\cdot(B|AX)=(B|X)\cdot(A|BX)$. + - We can make plausibility assignments by enriching the inductive + category $P^*$ over some monoidal category, e.g. the set of real numbers + (considered as a category) with its usual multiplication. /When we do/, + the identity arrows of $P^*$ \mdash{}corresponding to + deductive inferences\mdash{} are assigned a value of certainty automatically. + +[fn:1] /(IIIa) If a conclusion can be reasoned out in more than one +way, then every possible way must lead to the same result./ + + +*** Reciprocal probabilities +The natural numbers have a comparatively concrete origin: they are the +result of decategorifying the category of finite sets[fn:2], or the +coequalizer of the arrows from a one-object category to a two-object +category with a single nonidentity arrow. Extensions of the set of +natural numbers\mdash{} such as +the set of integers or rational numbers or real numbers\mdash{} strike +me as being somewhat more abstract. + +Jaynes points out that our existing choice of scale for probabilities +(i.e., the scale from 0 for impossibility to 1 for +certainty) has a degree of freedom: any monotonic function of +probability encodes the same information that probability does. + +With this in mind, it seems useful to use not /probability/ but +/reciprocal probability/ instead. This scale, which we might +tentatively call freeness, is a scale ranging 1 (certainty) to +positive infinity (impossibility). + +In this way, the ``probability'' +associated with choosing one out of $n$ indistinguishable choices +becomes identified with $n$. + +The entropy + +[fn:2] As Baez says. + + + +** self-questions + +What circumstances would make $\mathscr{F}$ an injection? + +What if $P=\{\top,\bot\}$? + + + +** COMMENT +Inductive and deductive posets are related as follows: there is a monotone +inclusion map $\mathscr{i}:P^*\hookrightarrow P$ which\mdash{} since $a$ +implies $b$ is equivalent to $a=ab$\mdash{} sends comparable +propositions in $P$ to the same proposition in $P^*$. Conversely, +only comparable propositions in $P$ are sent to the same proposition +in $P^*$. + + + +** Inductive posets and plausibility + +* Inverse Probability diff -r 8d8278e09888 -r b4de894a1e2e categorical/plausible.org_archive --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/plausible.org_archive Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,40 @@ +# -*- mode: org -*- + + +Archived entries from file /home/r/aurellem/src/categorical/plausible.org + +* Consistent reasoning as a commutative diagram + :PROPERTIES: + :ARCHIVE_TIME: 2011-07-09 Sat 01:00 + :ARCHIVE_FILE: ~/aurellem/src/categorical/plausible.org + :ARCHIVE_OLPATH: Deductive and inductive posets/Assigning plausibilities to inductive posets + :ARCHIVE_CATEGORY: plausible + :END: +Inductive categories enable the following neat trick: we can interpret +the objects of $P^*$ as states of given information and interpret +each arrow $a\rightarrow ab$ in $P^*$ as an inductive inference: the arrow +$a\rightarrow ab$ represents an inferential leap from the state of +knowledge where only $a$ is given to the state of knowledge where +both $a$ and $b$ are given\mdash{} in this way, it represents +the process of inferring $b$ when given $a$, and we label the +arrow with $(b|a)$. + +This trick has several important features that suggest its usefulness, +namely + - Composition of arrows corresponds to compound inference. + - In the special case of deductive inference, the inferential arrow is an + identity; the source and destination states of knowledge are the same. + - One aspect of the consistency requirement of Jaynes[fn:1] takes the form of a + commutative square: $x\rightarrow ax \rightarrow abx$ = + $x\rightarrow bx \rightarrow abx$ is the categorified version of + $(AB|X)=(A|X)\cdot(B|AX)=(B|X)\cdot(A|BX)$. + - We can make plausibility assignments by enriching the inductive + category $P^*$ over some monoidal category, e.g. the set of real numbers + (considered as a category) with its usual multiplication. /When we do/, + the identity arrows of $P^*$ \mdash{}corresponding to + deductive inferences\mdash{} are assigned a value of certainty automatically. + +[fn:1] /(IIIa) If a conclusion can be reasoned out in more than one +way, then every possible way must lead to the same result./ + + diff -r 8d8278e09888 -r b4de894a1e2e categorical/synthetic.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/synthetic.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,282 @@ + + + + +Synthetic Differential Geometry + + + + + + + + + + + + + +

+ + + + + +

Synthetic Differential Geometry

Written by Dylan Holmes

+ + + + + + +

+(My notes on Anders Kock's Synthetic Differential Geometry) +

+ +

1 Revisiting the real line +
+

+ +

1 Revisiting the real line

+ + +

+Lines, the kind which Euclid talked about, each constitute a commutative + ring: you choose any two points on the line to be 0 and 1, then add + and multiply as if you were dealing with real numbers $\mathbb{R}$. +

+Euclid moreover uses the axiom that for any two points, either they are the +same point or there is a unique line between them. Algebraically, +this amounts to saying that each line is not only a commutative ring +but a field, as well. This marks our first departure from euclidean +geometry, as our first axiom denies that each line is a field. +

+ + +

+ +

1.1 The first anti-euclidean axiom

+ +

A point in a ring is called nilpotent if its square is +zero. Normally (that is, in $\mathbb{R}^n$), only $0$ is +nilpotent. Here, as a consequence of the following axiom, there will +exist other elements that are nilpotent. These elements will +encapsulate our intuitive idea of “infinitesimally small” numbers. +

+ +
Axiom 1: Let $R$ be the line, considered as a commutative ring, and + let $D\subset R$ be the set of nilpotent elements on the line. Then for any + morphism $g:D\rightarrow R$, there exists a unique $b\in R$ such that +
+ + +$\forall d\in D, g(d) = g(0)+ b\cdot d$ + +
+Intuitively, this unique $b$ is the slope of the function $g$ near +zero. Because every morphism $g$ has exactly one such $b$, we have the +following results: +
+
+
The set $D$ of nilpotent elements contains more than + just 0. Indeed, suppose the contrary: if $D=\{0\}$, then for any $g$, every $b\in R$ has the + property described above;—$b$ isn't uniquely defined. +
+
Pick $b_1$ and $b_2$ in $R$. If every nilpotent $d$ satisfies $d\cdot + b_1 = d\cdot b_2$, then $b_1$ and $b_2$ are equal. +
+
+ + +

+ +

1.2 The first axiom $\ldots$ in terms of arrows

+ + +

+Define $\xi:R\times R\rightarrow R^D$ by $\xi:(a,b)\mapsto (d\mapsto +a+b\cdot d)$. The first axiom is equivalent to the statement +“ξ is invertible (i.e., a bijection)” +

+We give $R\times R$ the structure of an $R$-algebra by defining +multiplication: $ (a_1,b_1)\star(a_2,b_2) = (a_1\cdot a_2,\quad +a_1\cdot b_2 + a_2\cdot b_1)$. This is called dual-numbers multiplication, and is similar to muliplication of complex numbers. +

+ +

1.3 Ex

+ +

If $a$ and $b$ are nilpotent, then $ab$ is nilpotent. +
Even if $a$ and $b$ are nilpotent, the sum $a+b$ may not be. +
Even if $a+b$ is nilpotent, either summand $a$, $b$ may not be. +
+

+ + + + + +

Date: 2011-08-15 22:42:41 EDT

Author: Dylan Holmes

Org version 7.6 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e categorical/synthetic.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/categorical/synthetic.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,68 @@ +#+TITLE: Synthetic Differential Geometry +#+author: Dylan Holmes +#+EMAIL: rlm@mit.edu +#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js" +#+STYLE: +#+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t +#+SETUPFILE: ../templates/level-0.org +#+INCLUDE: ../templates/level-0.org +#+BABEL: :noweb yes :results silent + +(My notes on Anders Kock's /Synthetic Differential Geometry/) + +* Revisiting the real line + +*Lines*, the kind which Euclid talked about, each constitute a commutative + ring: you choose any two points on the line to be 0 and 1, then add + and multiply as if you were dealing with real numbers $\mathbb{R}$. + +Euclid moreover uses the axiom that for any two points, /either/ they are the +same point /or/ there is a unique line between them. Algebraically, +this amounts to saying that each line is not only a commutative ring +but a *field*, as well. This marks our first departure from euclidean +geometry, as our first axiom denies that each line is a field. + + +** The first anti-euclidean axiom +A point in a ring is called *nilpotent* if its square is +zero. Normally (that is, in $\mathbb{R}^n$), only $0$ is +nilpotent. Here, as a consequence of the following axiom, there will +exist other elements that are nilpotent. These elements will +encapsulate our intuitive idea of \ldquo{}infinitesimally small\rdquo{} numbers. + +#+begin_quote +*Axiom 1:* Let $R$ be the line, considered as a commutative ring, and + let $D\subset R$ be the set of nilpotent elements on the line. Then for any + morphism $g:D\rightarrow R$, there exists a unique $b\in R$ such that + +$\forall d\in D, g(d) = g(0)+ b\cdot d$ + +Intuitively, this unique $b$ is the slope of the function $g$ near +zero. Because every morphism $g$ has exactly one such $b$, we have the +following results: + +1. The set $D$ of nilpotent elements contains more than + just 0. Indeed, suppose the contrary: if $D=\{0\}$, then for any $g$, /every/ $b\in R$ has the + property described above;\mdash{}$b$ isn't uniquely defined. +2. Pick $b_1$ and $b_2$ in $R$. If every nilpotent $d$ satisfies $d\cdot + b_1 = d\cdot b_2$, then $b_1$ and $b_2$ are equal. + +** The first axiom $\ldots$ in terms of arrows + +Define $\xi:R\times R\rightarrow R^D$ by $\xi:(a,b)\mapsto (d\mapsto +a+b\cdot d)$. The first axiom is equivalent to the statement +\ldquo{}\xi is invertible (i.e., a bijection)\rdquo{} + +We give $R\times R$ the structure of an $R$-algebra by defining +multiplication: $ (a_1,b_1)\star(a_2,b_2) = (a_1\cdot a_2,\quad +a_1\cdot b_2 + a_2\cdot b_1)$. This is called *dual-numbers +multiplication*, and is similar to muliplication of complex numbers. + + +** Ex +1. If $a$ and $b$ are nilpotent, then $ab$ is nilpotent. +2. Even if $a$ and $b$ are nilpotent, the sum $a+b$ may not be. +3. Even if $a+b$ is nilpotent, either summand $a$, $b$ may not be. +4. + +#+end_quote diff -r 8d8278e09888 -r b4de894a1e2e chat/room.clj --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/chat/room.clj Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,34 @@ +(ns chat.room + (:use [clojure.contrib server-socket duck-streams str-utils]) + (:gen-class)) + + +;; CHAT ROOMS + +(def rooms (ref {})) + +(defn get-room! "Returns the chat room with the given name, creating it if necessary." [name] + (dosync + (or (@rooms name) + (let [new-room (agent '())] + (do (alter *rooms* assoc name new-room)))))) + +(defn say-in-room[room message] + (doseq [[_ output] room] + (binding [*out* output] + (println message)))) + + +;; USERS + +(defn user-join[room username output-channel] + (conj room [username output-channel])) +(defn user-leave[room username] + (remove #(= (% 0) username) room)) + + +(def *username* "Someone") +(defn- join-room [room] + (send room user-join *username* *out*) + (send room + ) diff -r 8d8278e09888 -r b4de894a1e2e clojure_magick/magick.clj --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/clojure_magick/magick.clj Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,151 @@ +(ns clojure-magick.magick) +(import [org.im4java.core ConvertCmd IMOperation Stream2BufferedImage] ) + +(defn read-image "Read the image in the specified file; returns a BufferedImage." [filename] + (let [op (doto (IMOperation.) + (.addImage (into-array String [filename])) + (.addImage (into-array String ["-"]))) + s2b (Stream2BufferedImage.)] + (doto (ConvertCmd.) + (.setOutputConsumer s2b) + (.run op (into-array Object []))) + (.getImage s2b)) +) + +(defn write-image "Write the image to the specified file." [img filename] + (let [op (doto (IMOperation.) + (.addImage) + (.addImage))] + (.run (ConvertCmd.) op (into-array Object [img filename]))) +) + + +(defn glue-vertical "Concatenate images in a sequence from top to bottom, producing a new image." [image & imgs] + (let [imgs (cons image imgs) + op (doto (IMOperation.) + (.addImage (count imgs)) + (.append) + (.addImage (into-array String ["-"]))) + s2b (Stream2BufferedImage.)] + (doto (ConvertCmd.) + (.setOutputConsumer s2b) + (.run op (into-array Object (vec imgs)))) + (.getImage s2b) +)) + +(defn glue-horizontal "Concatenate images in a sequence from left to right, producing a new image." [image & imgs] + (let [imgs (cons image imgs) + op (doto (IMOperation.) + (.addImage (count imgs)) + (.p_append) + (.addImage (into-array String ["-"]))) + s2b (Stream2BufferedImage.)] + (doto (ConvertCmd.) + (.setOutputConsumer s2b) + (.run op (into-array Object (vec imgs)))) + (.getImage s2b) + )) + + + +(defn resize "Resize the image, keeping the same aspect ratio." + ([img width height] + (let [op (doto (IMOperation.) + (.addImage) + (.resize width height) + (.addImage (into-array String ["-"]))) + s2b (Stream2BufferedImage.)] + (doto (ConvertCmd.) + (.setOutputConsumer s2b) + (.run op (into-array Object [img]))) + (.getImage s2b) + )) + ( + [img width] + (resize img width nil) + ) +) + + +(defn upside-down "Flip the image upside down." + [img] + (let [op (doto (IMOperation.) + (.addImage) + (.flip) + (.addImage (into-array String ["-"]))) + s2b (Stream2BufferedImage.)] + (doto (ConvertCmd.) + (.setOutputConsumer s2b) + (.run op (into-array Object [img]))) + (.getImage s2b) + ) +) + + +(defn scale-rotate-translate "Scale the image, then rotate the image about the axis coordinates, then translate the image by the given amount. Does not change the dimensions of the picture." + ([img axis-x axis-y scale-x scale-y degrees translate-x translate-y] + (let [arg-str (apply str (list axis-x "," axis-y " " scale-x "," scale-y " " degrees " " translate-x "," translate-y)) + op (doto (IMOperation.) + (.addImage) + (.distort "ScaleRotateTranslate" arg-str) + (.addImage (into-array String ["-"]))) + s2b (Stream2BufferedImage.)] + (doto (ConvertCmd.) + (.setOutputConsumer s2b) + (.run op (into-array Object [img]))) + (.getImage s2b))) + ([img axis-x axis-y scale degrees translate-x translate-y] + (recur axis-x axis-y scale scale degrees translate-x translate-y)) +) + + +(defn rotate "Rotate the image clockwise by the given amount." [img degrees] + (let [op (doto (IMOperation.) + (.addImage) + (.distort "ScaleRotateTranslate" (str degrees)) + (.addImage (into-array String ["-"]))) + s2b (Stream2BufferedImage.)] + (doto (ConvertCmd.) + (.setOutputConsumer s2b) + (.run op (into-array Object [img]))) + (.getImage s2b)) +) + + + + + + + + + + + + + + + + + + + + + + + + + + +(defn resize[] + (let [op (doto (IMOperation.) + (.addImage (into-array String ["/home/r/jiggly.gif"])) + (.resize 800 600) + (.addImage (into-array String ["/home/r/jiggly_sm.gif"])) + )] + (.run (ConvertCmd.) op (into-array Object [])))) + + +(defn run-test[] + (let [puff (read-image "/home/r/jiggly.gif")] + (write-image (glue-horizontal puff puff puff) "/home/r/multijiggly.gif") + )) diff -r 8d8278e09888 -r b4de894a1e2e mtg/bk.clj --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mtg/bk.clj Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,39 @@ +(ns mtg.frame) + +;; GENERALLY USEFUL FUNCTIONS + +(defn assay "Takes x and a series of pred-value pairs. Returns a list of vals for which the corresponding preds are true of x." [x & pred-vals] + (reduce #(if ((first %2) x) (conj %1 (second %2))) '() pred-vals) + ) +(defn alter-val "Applies f to the current value associated with each key, associating each key with the value returned." [m f & keys] + (map #(assoc m % (f (get m %))) keys)) + +(defn every-nth "Returns every nth member of coll. If n is not positive, returns an empty list." [n coll] + (if (<= n 0) '() + (take-while (comp not nil?) (map first (iterate #(nthnext % n) coll))))) + + + + + +;; FRAME MANIPULATION + + +(defn conj-key "Adds the xs to the seq associated with the given key." [map key & xs] + (assoc map key (apply conj (get map key []) xs))) + +(defn update "Takes a frame and a sequence of key-fn pairs. Applies f to the current value associated with key, updating the current value with the result. Frames generate and store a unique id for each call to update." + [frame & kfs] + (let [id (gensym "up_") + keys (every-nth 2 kfs) + fns (every-nth 2 (rest kfs))] + + ((reduce comp (map (fn[k f](fn[m](conj-key m k (list id f) )) keys fns)) + (conj-key frame :*bindings* (map (fn [k f](list id k)) keys fns)) + ) +)) + + + + +(def *frame* (atom {:*bindings* '()})) diff -r 8d8278e09888 -r b4de894a1e2e mtg/continuation.clj --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mtg/continuation.clj Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,159 @@ +(ns mtg.continuation) +;; convention: cps function names end in & + +(def not-nil? (comp not nil?)) +(defn cpt "Continuation-passing transform. Makes f take a continuation as an additional argument."[f] + (fn[& args] + ((last args) (apply f (butlast args))))) + + + +;; Card operations + + +;; Open predicates + ;(defmacro literal-to-keyword[x] `(keyword (quote ~x))) +(defrecord Oracle []) +(def *oracle* (atom (Oracle.))) + +(defn ask* [m k](get @m k (cpt(constantly nil)))) +(defn tell!* [m k v] (swap! m assoc k v)) + +(def ask(partial ask* *oracle*)) +;(def ? (fn [k & args] (apply (ask k) args))) + +(def tell! (partial tell!* *oracle*)) +(defn extend[k f] ((ask k) f)) + +;; PRELIMINARY DEFINITIONS +(defmacro defq "Defines a query." [name params & body] + `(tell! (keyword ~name) (fn ~params ~@body))) +(defmacro ? "Asks a query." [name & args] + `((ask (keyword name)) ~@args)) + +(defn true-preds "Returns a sequence of the preds for which (pred obj) returns true." + [obj & preds] + (map (comp (partial apply str) rest butlast str) + (filter #(? % obj) preds))) + +(tell! :the-players #{}) +(tell! :player? (fn[obj & _] (not-nil? (? the-players obj)))) + + + + + + + +(defq object? [obj & _] ;; 109.1 + ((comp not empty?)(true-preds obj :card? :copied-card? :token? :spell? :permanent? :emblem?))) + +(defq has-controller? [obj & _] + (if (and (not (? in-zone 'battlefield obj)) + (not (? on-stack obj))) false)) + +(defq take-turn [player & _] nil) + + +(tell! :characteristics (fn[obj & _] + (list :name :mana-cost :color :card-type :subtype :supertype :expansion-symbol :rules-text :abilities :power :toughness :loyalty :hand-modifier :life-modifier))) + + + + +(tell! :colors + (fn[obj & _] + (true-preds obj :red? :blue? :green? :white? :black?))) + + + + + + +;; GAME / TURN MECHANICS +(defn new-game "Start a new two-player game."[] + (tell! :the-players #{'PLAYERONE 'PLAYERTWO}) + (tell! :life-total (reduce #(assoc %1 (keyword %2) 20) {} (ask :the-players) )) +) + + + +;;(ask :blue) blue? = (fn[& args](fn[k] (k ...) )) +;; (cpt (constantly nil)) +;;(ask k) = (get self k (cpt(constantly nil))) + + + + + + +(defn reverse&[coll k] + (if (empty? coll) (k '()) + (recur (rest coll) (fn[x](k (conj x (first coll)))) + + ))) + + + + + + + + + + + + + + + + +(defn choose "Asks player to choose a member of the universe which matches all the given predicate; returns the chosen member." [player & preds] + + ) + + +(defn get* "Returns the value of key for the object, default or nil if key is not present." + ([obj key]) + ([obj key default])) + + + + + + + + + + + + + +;; shuffle decks +;; anyone may cut anyone else's deck +;; deck => library +;; decide who goes first +;; players decide whether to mulligan in APNAP +;; players mulligan simultaneously +;; ... finish mulligans +;; opening hand abilities +;; beginning game abilities +;; 2P game: first player skips draw phase of first turn + +;; pooled mana: produced-by, spendable-on + + + + +;; VALUES WITH MODIFICATION HISTORY +;(deftype trace[subject ]) + +;(defn trace [x]{:subject x :modifiers '()}) +;(defn modifiers[tr] (get :modifiers tr '())) +;(defn subject[tr] (:subject tr)) + +;(defn modify[tr f](assoc tr :modifiers (conj (modifiers tr) f))) +;(defn compute[tr] (:modifiers tr)) + +;; players + diff -r 8d8278e09888 -r b4de894a1e2e mtg/frame.clj --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mtg/frame.clj Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,50 @@ +(ns mtg.frame) + +;; GENERALLY USEFUL FUNCTIONS + +(defn assay "Takes x and a series of pred-value pairs. Returns a list of vals for which the corresponding preds are true of x." [x & pred-vals] + (reduce #(if ((first %2) x) (conj %1 (second %2))) '() pred-vals) + ) +(defn alter-val "Applies f to the current value associated with each key, associating each key with the value returned." [m f & keys] + (map #(assoc m % (f (get m %))) keys)) + +(defn every-nth "Returns every nth member of coll. If n is not positive, returns an empty list." [n coll] + (if (<= n 0) '() + (take-while (comp not nil?) (map first (iterate #(nthnext % n) coll))))) + + + + + +;; FRAME MANIPULATION + +(defn conj-key "Adds the xs to the seq associated with the given key." [map key & xs] + (assoc map key (apply conj (get map key []) xs))) + +(defn update "Takes a frame and a sequence of key-fn pairs. Applies f to the current value associated with key, updating the current value with the result. Frames generate and store a unique id for each call to update." + [frame & kfs] + (let [id (gensym "") + keys (every-nth 2 kfs) + fns (every-nth 2 (rest kfs))] + + ((reduce comp (map (fn[k f](fn[m](conj-key m k (list id f)))) keys fns)) + (conj-key frame :*bindings* (map (fn [k f](list id k)) keys fns)) + ) +)) + +(defn rollback "Undo the update with the given id." [frame id] + (let [affected-keys + (conj (map second (filter #(=(first %) id) (:*bindings* frame))) :*bindings*)] + (reduce (fn[frame key] + (alter-val (partial filter #(=(first %) id)) key) + ) frame affected-keys) + )) + + +(defn get-fn "Keys in a frame store lists of modifiers. Produces the end result of applying all the modifiers in order." [frame key] + (reduce #(%2) (constantly nil) (list (constantly 1))) +) + + +(def *frame* (atom {:*bindings* '()})) + diff -r 8d8278e09888 -r b4de894a1e2e mtg/mtg_cards.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mtg/mtg_cards.txt Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,42 @@ +:name "Mountain" +:color 'R +:type "Land" +:supertype "Basic" +:subtype "Mountain" +:rarity 'C + +:name "Island" +:color 'B +:type "Land" +:supertype "Basic" +:subtype "Island" +:rarity 'C + +:name "Forest" +:color 'G +:type "Land" +:supertype "Basic" +:subtype "Forest" +:rarity 'C + +:name "Plains" +:color 'W +:type "Land" +:supertype "Basic" +:subtype "Plains" +:rarity 'C + +:name "Swamp" +:color 'B +:type "Land" +:supertype "Basic" +:subtype "Swamp" +:rarity 'C + +:name "Knight Errant" +:color 'W +:cost "1W" +:type "Creature" +:subtype "Human Knight" +:p 2 +:t 2 diff -r 8d8278e09888 -r b4de894a1e2e org/science.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/org/science.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,44 @@ +#+title: Science Minus Science +#+author: Dylan Holmes +#+email: ocsenave@gmail.com +#+description: What's wrong with our current Science Education? +#+SETUPFILE: ../../aurellem/org/setup.org +#+INCLUDE: ../../aurellem/org/level-0.org + + +From what I've seen, today's science classrooms are remarkably +unscientific. Someone has decided that it is less important to teach +the empirical mindset than to impart our accumulated scientific +knowledge. Thus, because the field is so vast nowadays, teachers are +obliged to be frugal with the facts: they must prune tangential +subjects and pare whatever's left, watering down complicated results +into simplified half-truths. Needs must when the devil drives, of +course--but what is the end result? + +In modern science classrooms, we force-feed students a deluge of +unfamiliar scientific dogma which they must swallow in time to +regurgitate onto an exam. To accomplish this daunting task, they +cannot possibly stop to consider various alternatives which scientists +have methodically eliminated over the course of centuries; instead, +they must simply trust that science has done what it purports to have +done--or, faster, simply stamp out their conjectural, critical +instincts. + +By the end of such a course, students might be able to recite the +tenets of our current scientific creed and might employ those tenets +when answering carefully formulated questions. But even if, by chance, +our students get their facts straight, they will have acquired at most +only our pre-processed truths, and nothing of the empirical machinery +that produced them. In my opinion, such a lackluster result demands +that we re-evaluate our priorities. Surely the shibboleth of the +scientist is not his ability to recount the bleeding-edge depiction of +reality--after all, theories are transient and revolutions expected--but +rather his pervasive inquiries about the world and his methodical, +empirical approach to answering them? Indeed, don't we recognize the +scientist by his lack of allegiance to the status quo, by the way he +scrutinizes even his own theories with utmost irreverence? + +In valuing data absorption over methodical reason, we give our +students a fragmentary and moreover inexplicable impression of +reality. We must ask ourselves: how much of science is left in that? + diff -r 8d8278e09888 -r b4de894a1e2e proof.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/proof.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,151 @@ +#+TITLE: WFF 'n Proof +#+AUTHOR: Dylan Holmes +#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js" +#+STYLE: +#+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t +#+SETUPFILE: ../templates/level-0.org +#+INCLUDE: ../templates/level-0.org +#+BABEL: :noweb yes :results verbatim :cache yes + + +* Introduction +WFF 'n Proof is a competitive dice game designed to teach logic. The +game is played with several dice whose faces are marked with logical +symbols; the goal is to arrange the dice into /well-formed formulas/ (WFFs) \mdash{} +expressions that are syntactically correct\mdash{}then to manipulate +these formulas to form valid proofs. + +** About the dice +WFF and Proof has small cubes and large cubes which you use like +dice. The small cubes have lower-case letters =p=, =q=, =r=, =s=, =i=, +and =o= inscribed on their faces, whereas the big cubes have +upper-case letters =C=, =A=, =K=, =E=, =N= and =R= on their faces. The +lower-case letters stand for /propositions/ that can be either true or +false, whereas the +upper-case letters stand for certain logical connectives: + + +| Symbol | Meaning | +|--------+---------| +| =C= | Implies | +| =A= | Or | +| =K= | And | +| =E= | Equals | +|--------+---------| +| =N= | Not | + +** Well-formed formulas +An expression is a row of dice arranged left-to-right. An expression +is a well-formed formula (WFF) if and only if +- It is a =p=, =q=, =r=, or =s=; or +- It is an =N= followed by a WFF; or +- It is a =C=, =A=, =K=, or =E= followed by two WFFs. + +(Both the lower-case letters =i= and =o= and the upper-case letter =R= +are special, and are not used in forming WFFs) + +* Games you can play with WFF 'n Proof +** Shake a WFF +Shake a WFF is the simplest game; it teaches you to build and to +recognize well-formed formulas, skills that are fundamental to more advanced games. + +#This game is for two or three players. +Each player starts the game with two small cubes and one big cube; +these cubes constitute each player's hand. Divide the remaining cubes +evenly among the players so that each player has a pile of unused +cubes. Also allow space for completed WFFs. + +At the start of each turn, all of the players toss their dice +simultaneously; each player races to build the longest WFF possible +out of his own dice. As soon as a player is finished\mdash{}either he +has constructed the longest WFF possible or has found that it is impossible to make a WFF\mdash{}he +calls out \ldquo{}WFF!\rdquo{} or \ldquo{}No WFF!\rdquo{}. + +When a player calls out \ldquo{}WFF!\rdquo{} or \ldquo{}No +WFF!\rdquo{}, the other players must stop to validate the caller's +claim. When a player confirms that the claim is correct, that player +calls out \ldquo{}Check!\rdquo{}. Now, either all the players believe +that the claim is correct, or they do not. + + +- If all the players confirm that the claim is correct, the caller + puts his WFF in the area for completed WFFs, then replenishes his + supply of cubes by replacing each of the cubes in his WFF with a + cube from his unused-cubes pile (each big cube must be replaced with + a big cube, and each small cube must be replaced with a small cube), + then takes an additional cube from his pile and receives one + point. The turn then ends. + +- If, when checking the caller's claim, a player believes that the + claim was made in error, the player calls out + \ldquo{}Challenge!\rdquo{}, then points out the mistake: + - If the caller has presented an expression that is not a WFF, the + challenger shows that it is not a WFF. + - If the caller has presented a WFF that is not the longest possible, + the challenger shows how to make a longer WFF from the caller's + cubes. + - If it is possible to make a WFF from the caller's cubes but the + caller claims it is impossible, the challenger shows how to make a + WFF from the caller's cubes. + The other players verify this challenge. If the challenge was made + correctly, the challenging player wins a point and gains an extra + cube from his pile whereas the calling player must return a cube to his + pile. If the challenge was made wrongly, the challenger must return a + cube to his pile. In either case, the turn ends. + +After the turn ends, the player can play another turn. The game ends +when a certain number of cubes are in the area for completed WFFs; +this number can be decided before the game begins. + +One final note is the rule for cubes: +#+begin_quote +*Cube Rule:* A player's hand must have either the same number of big + cubes and small cubes, or one more small cube than big cube. When + taking cubes from or returning cubes to your unusued cubes pile, you + must make sure to follow this rule. +#+end_quote + + + + + + +* Misc + +#+srcname: proof +#+begin_src clojure +(ns wff.proof) + +(defn arity[x] + (first + (keep-indexed + (fn[i f](if (nil? (f x)) nil i)) + [#{'p 'q 'r 's} + #{'N} + #{'C 'A 'K 'E}]))) + +(defn wff? + [& symbols] + ((fn[stack symbols] + ;;(println stack symbols) + (if (or(empty? symbols)(empty? stack)) + (and(empty? symbols)(empty? stack)) + (let [n (arity (first symbols))] + (if (nil? n) + false + (recur (concat (rest stack) (repeat n '*)) (rest symbols)))))) + '(*) symbols)) + +#+end_src + + + + + + +* COMMENT tangle + +#+begin_src clojure :tangle ./proof.clj +<> +<> +#+end_src diff -r 8d8278e09888 -r b4de894a1e2e sicm/bk/utils.clj --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/bk/utils.clj Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,196 @@ + +(ns sicm.utils) + + ;***** GENERIC ARITHMETIC +(ns sicm.utils) +(in-ns 'sicm.utils) + +(defprotocol Arithmetic + (zero [this]) + (one [this])) + + +(extend-protocol Arithmetic + java.lang.Number + (zero [this] 0) + (one [this] 1)) + +(extend-protocol Arithmetic + clojure.lang.Seqable + (zero [this] (map zero this)) + (one [this] (map one this))) + + + ;***** TUPLES AND MATRICES +(in-ns 'sicm.utils) + +(defprotocol Spinning + (up? [this]) + (down? [this])) + +(defn spin + "Returns the spin of the Spinning s, either :up or :down" + [#^Spinning s] + (cond (up? s) :up (down? s) :down)) + + +(deftype Tuple + [spin coll] + clojure.lang.Seqable + (seq [this] (seq (.coll this))) + clojure.lang.Counted + (count [this] (count (.coll this)))) + +(extend-type Tuple + Spinning + (up? [this] (= ::up (.spin this))) + (down? [this] (= ::down (.spin this)))) + +(defmethod print-method Tuple + [o w] + (print-simple (str (if (up? o) 'u 'd) (.coll o)) w)) + + + +(defn up + "Create a new up-tuple containing the contents of coll." + [coll] + (Tuple. ::up coll)) + +(defn down + "Create a new down-tuple containing the contents of coll." + [coll] + (Tuple. ::down coll)) + + +(in-ns 'sicm.utils) + +(defn numbers? + "Returns true if all arguments are numbers, else false." + [& xs] + (every? number? xs)) + +(defn contractible? + "Returns true if the tuples a and b are compatible for contraction, + else false. Tuples are compatible if they have the same number of + components, they have opposite spins, and their elements are + pairwise-compatible." + [a b] + (and + (isa? (type a) Tuple) + (isa? (type b) Tuple) + (= (count a) (count b)) + (not= (spin a) (spin b)) + + (not-any? false? + (map #(or + (numbers? %1 %2) + (contractible? %1 %2)) + a b)))) + + + +(defn contract + "Contracts two tuples, returning the sum of the + products of the corresponding items. Contraction is recursive on + nested tuples." + [a b] + (if (not (contractible? a b)) + (throw + (Exception. "Not compatible for contraction.")) + (reduce + + (map + (fn [x y] + (if (numbers? x y) + (* x y) + (contract x y))) + a b)))) + + ;***** MATRICES +(in-ns 'sicm.utils) +(require 'incanter.core) ;; use incanter's fast matrices + +(defprotocol Matrix + (rows [this]) + (cols [this]) + (diagonal [this]) + (trace [this]) + (determinant [this])) + +(extend-protocol Matrix + incanter.Matrix + (rows [this] (map down this))) + + + + +(defn count-rows [matrix] + ((comp count rows) matrix)) + +(defn count-cols [matrix] + ((comp count cols) matrix)) + + +(defn matrix-by-rows + "Define a matrix by giving its rows." + [& rows] + (cond + (not (all-equal? (map count rows))) + (throw (Exception. "All rows in a matrix must have the same number of elements.")) + :else + (reify Matrix + (rows [this] (map down rows)) + (cols [this] (map up (apply map vector rows))) + (diagonal [this] (map-indexed (fn [i row] (nth row i) rows))) + (trace [this] + (if (not= (count-rows this) (count-cols this)) + (throw (Exception. + "Cannot take the trace of a non-square matrix.")) + (reduce + (diagonal this)))) + + (determinant [this] + (if (not= (count-rows this) (count-cols this)) + (throw (Exception. + "Cannot take the determinant of a non-square matrix.")) + (reduce * (diagonal this)))) + ))) + + +(defn matrix-by-cols + "Define a matrix by giving its columns." + [& cols] + (cond + (not (all-equal? (map count cols))) + (throw (Exception. "All columns in a matrix must have the same number of elements.")) + :else + (reify Matrix + (cols [this] (map up cols)) + (rows [this] (map down (apply map vector cols))) + (diagonal [this] (map-indexed (fn [i col] (nth col i) cols))) + (trace [this] + (if (not= (count-cols this) (count-rows this)) + (throw (Exception. + "Cannot take the trace of a non-square matrix.")) + (reduce + (diagonal this)))) + + (determinant [this] + (if (not= (count-cols this) (count-rows this)) + (throw (Exception. + "Cannot take the determinant of a non-square matrix.")) + (reduce * (map-indexed (fn [i col] (nth col i)) cols)))) + ))) + +(extend-protocol Matrix Tuple + (rows [this] (if (down? this) + (list this) + (map (comp up vector) this))) + + (cols [this] (if (up? this) + (list this) + (map (comp down vector) this)))) + +(defn matrix-multiply [A B] + (apply matrix-by-rows + (for [a (rows A)] + (for [b (cols B)] + (contract a b))))) diff -r 8d8278e09888 -r b4de894a1e2e sicm/bk/utils.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/bk/utils.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,1482 @@ + + + + +Building a Classical Mechanics Library in Clojure + + + + + + + + + + + + + +

+ + + + + +

Building a Classical Mechanics Library in Clojure

+ + + + + + + +

1 Generic Arithmetic
2 Useful Data Types +
- 2.1 Complex Numbers
- 2.2 Tuples and Tensors +
  - 2.2.1 Contraction
  - 2.2.2 Matrices
  +
- 2.3 Power Series
+
3 Basic Utilities +
- 3.1 Sequence manipulation
- 3.2 Ranges, Airity and Function Composition
+
4 Numerical Methods
5 Differentiation
6 Symbolic Manipulation

+ +

1 Generic Arithmetic

+ + + + + +

(ns sicm.utils)
+(in-ns 'sicm.utils)
+
+
+
+(defn all-equal? [coll]
+  (if (empty? (rest coll)) true
+      (and (= (first coll) (second coll))
+           (recur (rest coll)))))
+
+
+(defprotocol Arithmetic
+  (zero [this])
+  (one [this])
+  (negate [this])
+  (invert [this])
+  (conjugate [this]))
+
+(defn zero? [x]
+  (= x (zero x)))
+(defn one? [x]
+  (= x (one x)))
+
+
+
+(extend-protocol Arithmetic
+  java.lang.Number
+  (zero [x] 0)
+  (one [x] 1)
+  (negate [x] (- x))
+  (invert [x] (/ x))
+  )
+
+(extend-protocol Arithmetic
+  clojure.lang.IFn
+  (one [f] identity)
+  (negate [f] (comp negate f))
+  (invert [f] (comp invert f)))
+
+    
+(extend-protocol Arithmetic
+  clojure.lang.Seqable
+  (zero [this] (map zero this))
+  (one [this] (map one this))
+  (invert [this] (map invert this))
+  (negate [this] (map negate this)))
+
+
+(defn ordered-like
+  "Create a comparator using the sorted collection as an
+  example. Elements not in the sorted collection are sorted to the
+  end."
+  [sorted-coll]
+  (let [
+        sorted-coll? (set sorted-coll) 
+        ascending-pair?
+        (set(reduce concat
+                    (map-indexed
+                     (fn [n x]
+                       (map #(vector x %) (nthnext sorted-coll n)))
+                     sorted-coll)))]
+    (fn [x y]
+      (cond
+       (= x y) 0
+       (ascending-pair? [x y]) -1
+       (ascending-pair? [y x]) 1
+       (sorted-coll? x) -1
+       (sorted-coll? y) 1))))
+  
+
+
+(def type-precedence
+  (ordered-like [incanter.Matrix]))
+
+(defmulti add
+    (fn [x y]
+      (sort type-precedence [(type x)(type y)]))) 
+
+(defmulti multiply
+  (fn [x y]
+    (sort type-precedence [(type x) (type y)])))
+
+(defmethod add [java.lang.Number java.lang.Number] [x y] (+ x y))
+(defmethod multiply [java.lang.Number java.lang.Number] [x y] (* x y))
+
+(defmethod multiply [incanter.Matrix java.lang.Integer] [x y]
+           (let [args (sort #(type-precedence (type %1)(type %2)) [x y])
+                 matrix (first args)
+                 scalar (second args)]
+             (incanter.core/matrix (map (partial map (partial multiply scalar)) matrix))))
+           
+

+ + + + + +

+ +

2 Useful Data Types

+ + + +

+ +

2.1 Complex Numbers

+ + + + +

(in-ns 'sicm.utils)
+
+(defprotocol Complex
+  (real-part [z])
+  (imaginary-part [z])
+  (magnitude-squared [z])
+  (angle [z])
+  (conjugate [z])
+  (norm [z]))
+
+(defn complex-rectangular
+  "Define a complex number with the given real and imaginary
+  components."
+  [re im]
+  (reify Complex
+         (real-part [z] re)
+         (imaginary-part [z] im)
+         (magnitude-squared [z] (+ (* re re) (* im im)))
+         (angle [z] (java.lang.Math/atan2 im re))
+         (conjugate [z] (complex-rectangular re (- im)))
+
+         Arithmetic
+         (zero [z] (complex-rectangular 0 0))
+         (one [z] (complex-rectangular 1 0))
+         (negate [z] (complex-rectangular (- re) (- im)))
+         (invert [z] (complex-rectangular
+                      (/ re (magnitude-squared z))
+                      (/ (- im) (magnitude-squared z))))
+
+         Object
+         (toString [_]
+                   (if (and (zero? re) (zero? im)) (str 0)
+                       (str
+                        (if (not(zero? re))
+                          re)
+                        (if ((comp not zero?) im)
+                          (str
+                           (if (neg? im) "-" "+")
+                           (if ((comp not one?) (java.lang.Math/abs im))
+                           (java.lang.Math/abs im))
+                           "i")))))))
+
+(defn complex-polar
+  "Define a complex number with the given magnitude and angle."
+  [mag ang]
+  (reify Complex
+         (magnitude-squared [z] (* mag mag))
+         (angle [z] angle)
+         (real-part [z] (* mag (java.lang.Math/cos ang)))
+         (imaginary-part [z] (* mag (java.lang.Math/sin ang)))
+         (conjugate [z] (complex-polar mag (- ang)))
+
+         Arithmetic
+         (zero [z] (complex-polar 0 0))
+         (one [z] (complex-polar 1 0))
+         (negate [z] (complex-polar (- mag) ang))
+         (invert [z] (complex-polar (/ mag) (- ang)))
+
+         Object
+         (toString [_] (str mag " * e^(i" ang")"))
+         ))
+
+
+;; Numbers are complex quantities
+
+(extend-protocol Complex
+  java.lang.Number
+  (real-part [x] x)
+  (imaginary-part [x] 0)
+  (magnitude [x] x)
+  (angle [x] 0)
+  (conjugate [x] x))
+
+
+

+ + + + + + +

+ +

2.2 Tuples and Tensors

+ + +

+A tuple is a vector which is spinable—it can be either spin up or spin down. (Covariant, contravariant; dual vectors) +

+ + +

(in-ns 'sicm.utils)
+
+(defprotocol Spinning
+  (up? [this])
+  (down? [this]))
+
+(defn spin
+  "Returns the spin of the Spinning s, either :up or :down"
+  [#^Spinning s]
+  (cond (up? s) :up (down? s) :down))
+
+
+(deftype Tuple
+  [spin coll]
+  clojure.lang.Seqable
+  (seq [this] (seq (.coll this)))
+  clojure.lang.Counted
+  (count [this] (count (.coll this))))
+
+(extend-type Tuple
+  Spinning
+  (up? [this] (= ::up (.spin this)))
+  (down? [this] (= ::down (.spin this))))
+
+(defmethod print-method Tuple
+  [o w]
+  (print-simple (str (if (up? o) 'u 'd) (.coll o))  w))
+
+
+
+(defn up
+  "Create a new up-tuple containing the contents of coll."
+  [coll]
+  (Tuple. ::up coll))       
+
+(defn down
+  "Create a new down-tuple containing the contents of coll."
+  [coll]
+  (Tuple. ::down coll))
+
+
+

+ + + + + +

+ +

2.2.1 Contraction

+ +

Contraction is a binary operation that you can apply to compatible + tuples. Tuples are compatible for contraction if they have the same + length and opposite spins, and if the corresponding items in each + tuple are both numbers or both compatible tuples. +

+ + + +

(in-ns 'sicm.utils)
+
+(defn numbers?
+  "Returns true if all arguments are numbers, else false."
+  [& xs]
+  (every? number? xs))
+
+(defn contractible?
+  "Returns true if the tuples a and b are compatible for contraction,
+  else false. Tuples are compatible if they have the same number of
+  components, they have opposite spins, and their elements are
+  pairwise-compatible."
+  [a b]
+  (and
+   (isa? (type a) Tuple)
+   (isa? (type b) Tuple)
+   (= (count a) (count b))
+   (not= (spin a) (spin b))
+   
+   (not-any? false?
+             (map #(or
+                    (numbers? %1 %2)
+                    (contractible? %1 %2))
+                  a b))))
+
+
+
+(defn contract
+  "Contracts two tuples, returning the sum of the
+  products of the corresponding items. Contraction is recursive on
+  nested tuples."
+  [a b]
+  (if (not (contractible? a b))
+    (throw
+     (Exception. "Not compatible for contraction."))
+    (reduce +
+            (map
+             (fn [x y]
+               (if (numbers? x y)
+                 (* x y)
+                 (contract x y)))
+             a b))))
+
+

+ + + + +

+ +

2.2.2 Matrices

+ + + + +

(in-ns 'sicm.utils)
+(require 'incanter.core) ;; use incanter's fast matrices
+
+(defprotocol Matrix
+  (rows [matrix])
+  (cols [matrix])
+  (diagonal [matrix])
+  (trace [matrix])
+  (determinant [matrix])
+  (transpose [matrix])
+  (conjugate [matrix])
+)
+
+(extend-protocol Matrix
+  incanter.Matrix
+  (rows [rs] (map down (apply map vector (apply map vector rs))))
+  (cols [rs] (map up (apply map vector rs)))
+  (diagonal [matrix] (incanter.core/diag matrix) )
+  (determinant [matrix] (incanter.core/det matrix))
+  (trace [matrix] (incanter.core/trace matrix))
+  (transpose [matrix] (incanter.core/trans matrix))
+  )
+
+(defn count-rows [matrix]
+  ((comp count rows) matrix))
+
+(defn count-cols [matrix]
+  ((comp count cols) matrix))
+
+(defn square? [matrix]
+  (= (count-rows matrix) (count-cols matrix)))
+
+(defn identity-matrix
+  "Define a square matrix of size n-by-n with 1s along the diagonal and
+  0s everywhere else."
+  [n]
+  (incanter.core/identity-matrix n))
+
+
+
+
+
+(defn matrix-by-rows
+  "Define a matrix by giving its rows."
+  [& rows]
+  (if
+   (not (all-equal? (map count rows)))
+   (throw (Exception. "All rows in a matrix must have the same number of elements."))
+   (incanter.core/matrix (vec rows))))
+
+(defn matrix-by-cols
+  "Define a matrix by giving its columns"
+  [& cols]
+  (if (not (all-equal? (map count cols)))
+   (throw (Exception. "All columns in a matrix must have the same number of elements."))
+   (incanter.core/matrix (vec (apply map vector cols)))))
+
+(defn identity-matrix
+  "Define a square matrix of size n-by-n with 1s along the diagonal and
+  0s everywhere else."
+  [n]
+  (incanter.core/identity-matrix n))
+
+
+
+(extend-protocol Arithmetic
+  incanter.Matrix
+  (one [matrix]
+       (if (square? matrix)
+         (identity-matrix (count-rows matrix))
+         (throw (Exception. "Non-square matrices have no multiplicative unit."))))
+  (zero [matrix]
+        (apply matrix-by-rows (map zero (rows matrix))))
+  (negate [matrix]
+          (apply matrix-by-rows (map negate (rows matrix))))
+  (invert [matrix]
+          (incanter.core/solve matrix)))
+
+
+
+(defmulti coerce-to-matrix
+ "Converts x into a matrix, if possible."
+ type)
+
+(defmethod coerce-to-matrix incanter.Matrix [x] x)
+(defmethod coerce-to-matrix Tuple [x]
+           (if (apply numbers? (seq x))
+             (if (up? x)
+             (matrix-by-cols (seq x))
+             (matrix-by-rows (seq x)))
+             (throw (Exception. "Non-numerical tuple cannot be converted into a matrix.")))) 
+             
+  
+
+
+
+;; (defn matrix-by-cols
+;;   "Define a matrix by giving its columns."
+;;   [& cols]
+;;   (cond
+;;    (not (all-equal? (map count cols)))
+;;    (throw (Exception. "All columns in a matrix must have the same number of elements."))
+;;    :else
+;;    (reify Matrix
+;;        (cols [this] (map up cols))
+;;        (rows [this] (map down (apply map vector cols)))
+;;        (diagonal [this] (map-indexed (fn [i col] (nth col i) cols)))
+;;        (trace [this]
+;;               (if (not= (count-cols this) (count-rows this))
+;;                 (throw (Exception.
+;;                         "Cannot take the trace of a non-square matrix."))
+;;                 (reduce + (diagonal this))))
+          
+;;        (determinant [this]
+;;                     (if (not= (count-cols this) (count-rows this))
+;;                       (throw (Exception.
+;;                               "Cannot take the determinant of a non-square matrix."))
+;;                       (reduce * (map-indexed (fn [i col] (nth col i)) cols))))
+;;        )))
+
+(extend-protocol Matrix  Tuple
+                 (rows [this] (if (down? this)
+                                (list this)
+                                (map (comp up vector) this)))
+
+                 (cols [this] (if (up? this)
+                                (list this)
+                                (map (comp down vector) this))
+                 ))
+  
+(defn matrix-multiply
+  "Returns the matrix resulting from the matrix multiplication of the given arguments."
+  ([A] (coerce-to-matrix A))
+  ([A B] (incanter.core/mmult (coerce-to-matrix A) (coerce-to-matrix B)))
+  ([M1 M2 & Ms] (reduce matrix-multiply (matrix-multiply M1 M2) Ms)))
+
+

+ + + + + +

+ +

2.3 Power Series

+ + + + +

(in-ns 'sicm.utils)
+(use 'clojure.contrib.def)
+
+
+
+(defn series-fn
+  "The function corresponding to the given power series."
+  [series]
+  (fn [x]
+    (reduce +
+            (map-indexed  (fn[n x] (* (float (nth series n)) (float(java.lang.Math/pow (float x) n)) ))
+                          (range 20)))))
+
+(deftype PowerSeries
+  [coll]
+  clojure.lang.Seqable
+  (seq [this] (seq (.coll this)))
+  
+  clojure.lang.Indexed
+  (nth [this n] (nth (.coll this) n 0))
+  (nth [this n not-found] (nth (.coll this) n not-found))
+
+  ;; clojure.lang.IFn
+  ;; (call [this] (throw(Exception.)))
+  ;; (invoke [this & args] args
+  ;;      (let [f 
+  ;; )
+  ;; (run [this] (throw(Exception.)))
+  )
+
+(defn power-series
+  "Returns a power series with the items of the coll as its
+  coefficients. Trailing zeros are added to the end of coll."
+  [coeffs]
+  (PowerSeries. coeffs))
+
+(defn power-series-indexed
+  "Returns a power series consisting of the result of mapping f to the non-negative integers."
+  [f]
+  (PowerSeries. (map f (range))))
+
+
+(defn-memo nth-partial-sum
+  ([series n]
+     (if (zero? n) (first series)
+         (+ (nth series n)
+            (nth-partial-sum series (dec n))))))
+
+(defn partial-sums [series]
+  (lazy-seq (map nth-partial-sum (range))))
+
+
+
+
+(def cos-series
+     (power-series-indexed
+      (fn[n]
+        (if (odd? n) 0
+            (/
+             (reduce *
+                     (reduce * (repeat (/ n 2) -1))
+                     (range 1 (inc n)))
+             )))))
+
+(def sin-series
+     (power-series-indexed
+      (fn[n]
+        (if (even? n) 0
+            (/
+             (reduce *
+                     (reduce * (repeat (/ (dec n) 2) -1))
+                     (range 1 (inc n)))
+             )))))
+
+

+ + + + + +

+ +

3 Basic Utilities

+ + + +

+ +

3.1 Sequence manipulation

+ + + + + +

(ns sicm.utils)
+
+(defn do-up
+  "Apply f to each number from low to high, presumably for
+  side-effects."
+  [f low high]
+  (doseq [i (range low high)] (f i)))
+   
+(defn do-down
+  "Apply f to each number from high to low, presumably for
+   side-effects."
+  [f high low]
+  (doseq [i (range high low -1)] (f i)))
+
+
+(defn all-equal? [coll]
+  (if (empty? (rest coll)) true
+      (and (= (first coll) (second coll))
+           (recur (rest coll))))))
+
+(defn multiplier
+  "Returns a function that 'multiplies' the members of a collection,
+returning unit if called on an empty collection."
+  [multiply unit]
+  (fn [coll] ((partial reduce multiply unit) coll)))
+
+(defn divider
+  "Returns a function that 'divides' the first element of a collection
+by the 'product' of the rest of the collection."
+  [divide multiply invert unit]
+  (fn [coll]
+    (apply
+     (fn
+       ([] unit)
+       ([x] (invert x))
+       ([x y] (divide x y))
+       ([x y & zs] (divide x (reduce multiply y zs))))
+     coll)))
+
+(defn left-circular-shift
+  "Remove the first element of coll, adding it to the end of coll."
+  [coll]
+  (concat (rest coll) (take 1 coll)))
+
+(defn right-circular-shift
+  "Remove the last element of coll, adding it to the front of coll."
+  [coll]
+  (cons (last coll) (butlast coll)))
+

+ + + + + + + +

+ +

3.2 Ranges, Airity and Function Composition

+ + + + +

(in-ns 'sicm.utils)
+(def infinity Double/POSITIVE_INFINITY)
+(defn infinite? [x] (Double/isInfinite x))
+(def finite? (comp not infinite?)) 
+
+(defn arity-min
+  "Returns the smallest number of arguments f can take."
+  [f]
+  (apply
+   min
+   (map (comp alength #(.getParameterTypes %))
+        (filter (comp (partial = "invoke") #(.getName %))
+                (.getDeclaredMethods (class f))))))
+  
+(defn arity-max
+  "Returns the largest number of arguments f can take, possibly
+  Infinity."
+  [f]
+  (let [methods (.getDeclaredMethods (class f))]
+    (if (not-any? (partial = "doInvoke") (map #(.getName %) methods))
+      (apply max
+             (map (comp alength #(.getParameterTypes %))
+                  (filter (comp (partial = "invoke") #(.getName %))  methods)))
+      infinity)))
+
+
+(def ^{:arglists '([f])
+       :doc "Returns a two-element list containing the minimum and
+     maximum number of args that f can take."}
+     arity-interval
+     (juxt arity-min arity-max))
+
+
+
+;; --- intervals
+
+(defn intersect
+  "Returns the interval of overlap between interval-1 and interval-2"
+  [interval-1 interval-2]
+  (if (or (empty? interval-1) (empty? interval-2)) []
+      (let [left (max (first interval-1) (first interval-2))
+            right (min (second interval-1) (second interval-2))]
+        (if (> left right) []
+            [left right]))))
+
+(defn same-endpoints?
+  "Returns true if the left endpoint is the same as the right
+  endpoint."
+  [interval]
+  (= (first interval) (second interval)))
+
+(defn naturals?
+  "Returns true if the left endpoint is 0 and the right endpoint is
+infinite."
+  [interval]
+  (and (zero? (first interval))
+       (infinite? (second interval))))
+ 
+
+(defn fan-in
+  "Returns a function that pipes its input to each of the gs, then
+  applies f to the list of results. Consequently, f must be able to
+  take a number of arguments equal to the number of gs."
+  [f & gs]
+  (fn [& args]
+    (apply f (apply (apply juxt gs) args))))
+
+(defn fan-in
+  "Returns a function that pipes its input to each of the gs, then applies f to the list of results. The resulting function takes any number of arguments, but will fail if given arguments that are incompatible with any of the gs."
+  [f & gs]
+  (comp (partial apply f) (apply juxt gs)))
+
+
+
+(defmacro airty-blah-sad [f n more?]
+  (let [syms (vec (map (comp gensym (partial str "x")) (range n)))
+         optional (gensym "xs")]
+    (if more?
+      `(fn ~(conj syms '& optional)
+         (apply ~f ~@syms ~optional))
+      `(fn ~syms (~f ~@syms)))))
+
+(defmacro airt-whaa* [f n more?]
+  `(airty-blah-sad ~f ~n ~more?))
+
+
+
+
+(defn fan-in*
+  "Returns a function that pipes its input to each of the gs, then
+  applies f to the list of results. Unlike fan-in, fan-in* strictly
+  enforces arity: it will fail if the gs do not have compatible
+  arities."
+  [f & gs]
+  (let [arity-in (reduce intersect (map arity-interval gs))
+        left (first arity-in)
+        right (second arity-in)
+        composite (fan-in f gs)
+        ]
+    (cond
+     (empty? arity-in)
+     (throw (Exception. "Cannot compose functions with incompatible arities."))
+
+     (not
+      (or (= left right)
+          (and (finite? left)
+               (= right infinity))))
+         
+     (throw (Exception.
+             "Compose can only handle arities of the form [n n] or [n infinity]"))
+     :else
+     (airty-blah-sad composite left (= right infinity)))))
+    
+
+
+(defn compose-n "Compose any number of functions together."
+  ([] identity)
+  ([f] f)
+  ([f & fs]
+  (let [fns (cons f fs)]
+    (compose-bin (reduce fan-in (butlast fs)) (last fs))))
+)
+
+
+
+
+           
+
+(defn iterated
+  ([f n id] (reduce comp id (repeat n f)))
+  ([f n] (reduce comp identity (repeat n f))))
+
+(defn iterate-until-stable
+  "Repeatedly applies f to x, returning the first result that is close
+enough to its predecessor."
+  [f close-enough? x]
+  (second (swank.util/find-first
+           (partial apply close-enough?)
+           (partition 2 1 (iterate f x)))))
+
+(defn lexical< [x y]
+  (neg? (compare (str x) (str y))))
+
+
+;; do-up
+;; do-down
+(def make-pairwise-test comparator)
+;;all-equal?
+(def accumulation multiplier)
+(def inverse-accumulation divider)
+;;left-circular-shift
+;;right-circular-shift
+(def exactly-n? same-endpoints?)
+(def any-number? naturals?)
+;; TODO compose
+;; TODO compose-n
+;; identity
+(def compose-2 fan-in)
+(def compose-bin fan-in*)
+(def any? (constantly true))
+(def none? (constantly false))
+(def constant constantly)
+(def joint-arity intersect)
+(def a-reduce reduce)
+;; filter
+(def make-map (partial partial map) )
+(def bracket juxt)
+;; TODO apply-to-all
+;; TODO nary-combine
+;; TODO binary-combine
+;; TODO unary-combine
+;; iterated
+;; iterate-until-stable
+(def make-function-of-vector (partial partial map))
+(def make-function-of-arguments (fn [f] (fn [& args] (f args))))
+(def alphaless lexical<)
+
+

+ + + + + + + + + + + + +

+ +

4 Numerical Methods

+ + + + +

(in-ns 'sicm.utils)
+(import java.lang.Math)
+(use 'clojure.contrib.def)
+
+;; ---- USEFUL CONSTANTS
+
+(defn machine-epsilon
+  "Smallest value representable on your machine, as determined by
+successively dividing a number in half until consecutive results are
+indistinguishable."
+  []
+  (ffirst
+   (drop-while
+    (comp not zero? second)
+    (partition 2 1
+               (iterate (partial * 0.5) 1)))))
+
+
+(def pi (Math/PI))
+(def two-pi (* 2 pi))
+
+(def eulers-gamma 0.5772156649015328606065)
+
+(def phi (/ (inc (Math/sqrt 5)) 2))
+
+(def ln2 (Math/log 2))
+(def ln10 (Math/log 10))
+(def exp10 #(Math/pow 10 %))
+(def exp2 #(Math/pow 2 %))
+
+
+;;
+
+;; ---- ANGLES AND TRIGONOMETRY
+
+(defn angle-restrictor
+  "Returns a function that ensures that angles lie in the specified interval of length two-pi."
+  [max-angle]
+  (let [min-angle (- max-angle two-pi)]
+    (fn [x]
+      (if (and
+           (<= min-angle x)
+           (< x max-angle))
+        x
+        (let [corrected-x (- x (* two-pi (Math/floor (/ x two-pi))))]
+          (if (< corrected-x max-angle)
+            corrected-x
+            (- corrected-x two-pi)))))))
+
+(defn angle-restrict-pi
+  "Coerces angles to lie in the interval from -pi to pi."
+  [angle]
+  ((angle-restrictor pi) angle))
+
+(defn angle-restrict-two-pi
+  "Coerces angles to lie in the interval from zero to two-pi"
+  [angle]
+  ((angle-restrictor two-pi) angle))
+
+
+
+
+(defn invert [x] (/ x))
+(defn negate [x] (- x))
+
+(defn exp [x] (Math/exp x))
+
+(defn sin [x] (Math/sin x))
+(defn cos [x] (Math/cos x))
+(defn tan [x] (Math/tan x))
+
+(def sec (comp invert cos))
+(def csc (comp invert sin))
+
+(defn sinh [x] (Math/sinh x))
+(defn cosh [x] (Math/cosh x))
+(defn tanh [x] (Math/tanh x))
+
+(def sech (comp invert cosh))
+(def csch (comp invert sinh))
+
+
+;; ------------
+
+(defn factorial
+  "Computes the factorial of the nonnegative integer n."
+  [n]
+  (if (neg? n)
+    (throw (Exception. "Cannot compute the factorial of a negative number."))
+    (reduce * 1 (range 1 (inc n)))))
+
+(defn exact-quotient [n d] (/ n d))
+
+(defn binomial-coefficient
+  "Computes the number of different ways to choose m elements from n."
+  [n m]
+  (assert (<= 0 m n))
+  (let [difference (- n m)]
+    (exact-quotient
+     (reduce * (range n (max difference m) -1 ))
+     (factorial (min difference m)))))
+
+(defn-memo stirling-1
+  "Stirling numbers of the first kind: the number of permutations of n
+  elements with exactly m permutation cycles. "
+  [n k]
+  ;(assert (<= 1 k n))
+  (if (zero? n)
+    (if (zero? k) 1 0)
+    (+ (stirling-1 (dec n) (dec k))
+       (* (dec n) (stirling-1 (dec n) k)))))
+
+(defn-memo stirling-2 ;;count-partitions
+  "Stirling numbers of the second kind: the number of ways to partition a set of n elements into k subsets."
+  [n k]
+  (cond
+   (= k 1) 1
+   (= k n) 1
+   :else (+ (stirling-2 (dec n) (dec k))
+            (* k (stirling-2 (dec n) k)))))
+
+(defn harmonic-number [n]
+  (/ (stirling-1 (inc n) 2)
+     (factorial n)))
+
+
+(defn sum
+  [f low high]
+  (reduce + (map f (range low (inc high)))))
+
+

+ + + + + + + + + + + + + + + +

+ +

5 Differentiation

+ + +

+We compute derivatives by passing special differential objects $[x, +dx]$ through functions. Roughly speaking, applying a function $f$ to a +differential object $[x, dx]$ should produce a new differential +object $[f(x),\,Df(x)\cdot dx]$. +

+ + +$[x,\,dx]\xrightarrow{\quad f \quad}[f(x),\,Df(x)\cdot dx]$ +

+Notice that you can obtain the derivative of $f$ from this +differential object, as it is the coefficient of the $dx$ term. Also, +as you apply successive functions using this rule, you get the +chain-rule answer you expect: +

+ + +$[f(x),\,Df(x)\cdot dx]\xrightarrow{\quad g\quad} [gf(x),\, +Dgf(x)\cdot Df(x) \cdot dx ]$ + +

+In order to generalize to multiple variables and multiple derivatives, +we use a power series of differentials, a sortred infinite sequence which +contains all terms like $dx\cdot dy$, $dx^2\cdot dy$, etc. +

+ + + + +

(in-ns 'sicm.utils)
+
+(defn differential-seq
+  "Constructs a sequence of differential terms from a numerical
+coefficient and a list of keys for variables. If no coefficient is supplied, uses 1."
+  ([variables] (differential-seq* 1 variables))
+  ([coefficient variables & cvs]
+     (if (number? coefficient)
+       (conj (assoc {} (apply sorted-set variables) coefficient)
+             (if (empty? cvs)
+               nil
+               (apply differential-seq* cvs)))
+       (apply differential-seq* 1 coefficient 1 variables cvs)  
+       )))
+
+
+
+(defn differential-add
+  "Add two differential sequences by combining like terms."
+  [dseq1 dseq2]
+  (merge-with + dseq1 dseq2))
+
+(defn differential-multiply
+  "Multiply two differential sequences. The square of any differential variable is zero since differential variables are infinitesimally small."
+  [dseq1 dseq2]
+  (reduce
+   (fn [m [[vars1 coeff1] [vars2 coeff2]]]
+     (if (empty? (clojure.set/intersection vars1 vars2))
+       (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))
+       m))
+   {}
+  (clojure.contrib.combinatorics/cartesian-product
+   dseq1
+   dseq2)))
+
+
+
+(defn big-part
+  "Returns the 'finite' part of the differential sequence."
+  [dseq]
+  (let
+      [keys (sort-by count (keys dseq))
+       pivot-key (first (last keys))]
+
+    (apply hash-map
+    (reduce concat
+            (filter (comp pivot-key first) dseq)
+            ))))
+
+
+(defn small-part
+  "Returns the 'infinitesimal' part of the differential sequence."
+  [dseq]
+    (let
+      [keys (sort-by count (keys dseq))
+       pivot-key (first (last keys))]
+
+    (apply hash-map
+    (reduce concat
+            (remove (comp pivot-key first) dseq)
+            ))))
+
+
+
+
+
+
+
+
+
+;; ;; A differential term consists of a numerical coefficient and a
+;; ;; sorted  
+;; (defrecord DifferentialTerm [coefficient variables])
+;; (defmethod print-method DifferentialTerm
+;;    [o w]
+;;     (print-simple
+;;      (apply str (.coefficient o)(map (comp (partial str "d") name) (.variables o)))
+;;      w))
+
+
+;; (defn differential-seq
+;;   "Constructs a sequence of differential terms from a numerical
+;; coefficient and a list of keywords for variables. If no coefficient is
+;; supplied, uses 1."
+;;   ([variables] (differential-seq 1 variables))
+;;   ([coefficient variables]
+;;      (list
+;;       (DifferentialTerm. coefficient (apply sorted-set variables))))
+;;   ([coefficient variables & cvs]
+;;      (sort-by
+;;       #(vec(.variables %))
+;;       (concat (differential-seq coefficient variables) (apply differential-seq cvs)))))
+
+

+ + + + + + + + + + + + + + + + +

+ +

6 Symbolic Manipulation

+ + + + + +

(in-ns 'sicm.utils)
+
+(deftype Symbolic [type expression]
+  Object
+  (equals [this that]
+          (cond
+           (= (.expression this) (.expression that)) true
+           :else
+           (Symbolic.
+            java.lang.Boolean
+            (list '= (.expression this) (.expression that)))
+          )))
+
+
+
+
+(deftype AbstractSet [glyph membership-test])
+(defn member? [abstract-set x]
+  ((.membership-test abstract-set) x))
+
+;; ------------ Some important AbstractSets
+
+
+(def Real
+     (AbstractSet.
+      'R
+      (fn[x](number? x))))
+
+
+;; ------------ Create new AbstractSets from existing ones
+
+(defn abstract-product
+  "Gives the cartesian product of abstract sets."
+  ([sets]
+       (if (= (count sets) 1) (first sets)
+           (AbstractSet.
+            (symbol
+             (apply str
+                    (interpose 'x (map #(.glyph %) sets))))
+           (fn [x]
+             (and
+              (coll? x)
+              (= (count sets) (count x))
+              (reduce #(and %1 %2)
+                      true
+                      (map #(member? %1 %2) sets x)))))))
+  ([abstract-set n]
+     (abstract-product (repeat n abstract-set))))
+
+
+              
+(defn abstract-up
+  "Returns the abstract set of all up-tuples whose items belong to the
+  corresponding abstract sets in coll."
+  ([coll]
+     (AbstractSet.
+      (symbol (str "u["
+                   (apply str
+                          (interpose " "
+                                     (map #(.glyph %) coll)))
+                   "]"))
+      (fn [x]
+        (and
+         (satisfies? Spinning x)
+         (up? x)
+         (= (count coll) (count x))
+         (reduce
+          #(and %1 %2)
+          true
+          (map #(member? %1 %2) coll x))))))
+  ([abstract-set n]
+     (abstract-up (repeat n abstract-set))))
+
+
+(defn abstract-down
+  "Returns the abstract set of all down-tuples whose items belong to the
+  corresponding abstract sets in coll."
+  ([coll]
+     (AbstractSet.
+      (symbol (str "d["
+                   (apply str
+                          (interpose " "
+                                     (map #(.glyph %) coll)))
+                   "]"))
+      (fn [x]
+        (and
+         (satisfies? Spinning x)
+         (down? x)
+         (= (count coll) (count x))
+         (reduce
+          #(and %1 %2)
+          true
+          (map #(member? %1 %2) coll x))))))
+  ([abstract-set n]
+     (abstract-down (repeat n abstract-set))))
+
+              
+
+
+
+                                        ;-------ABSTRACT FUNCTIONS
+(defrecord AbstractFn
+  [#^AbstractSet domain #^AbstractSet codomain])
+
+
+(defmethod print-method AbstractFn
+  [o w]
+  (print-simple
+   (str
+    "f:"
+   (.glyph (:domain o))
+   "-->"
+   (.glyph (:codomain o))) w))
+

+ + + + + + + +

Date: 2011-08-09 18:41:37 EDT

Author: Robert McIntyre & Dylan Holmes

Org version 7.6 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e sicm/bk/utils.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/bk/utils.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,1262 @@ +#+TITLE:Building a Classical Mechanics Library in Clojure +#+author: Robert McIntyre & Dylan Holmes +#+EMAIL: rlm@mit.edu +#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js" +#+STYLE: +#+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t +#+SETUPFILE: ../templates/level-0.org +#+INCLUDE: ../templates/level-0.org +#+BABEL: :noweb yes :results silent + +* Generic Arithmetic + +#+srcname: generic-arithmetic +#+begin_src clojure +(ns sicm.utils) +(in-ns 'sicm.utils) + + + +(defn all-equal? [coll] + (if (empty? (rest coll)) true + (and (= (first coll) (second coll)) + (recur (rest coll))))) + + +(defprotocol Arithmetic + (zero [this]) + (one [this]) + (negate [this]) + (invert [this]) + (conjugate [this])) + +(defn zero? [x] + (= x (zero x))) +(defn one? [x] + (= x (one x))) + + + +(extend-protocol Arithmetic + java.lang.Number + (zero [x] 0) + (one [x] 1) + (negate [x] (- x)) + (invert [x] (/ x)) + ) + +(extend-protocol Arithmetic + clojure.lang.IFn + (one [f] identity) + (negate [f] (comp negate f)) + (invert [f] (comp invert f))) + + +(extend-protocol Arithmetic + clojure.lang.Seqable + (zero [this] (map zero this)) + (one [this] (map one this)) + (invert [this] (map invert this)) + (negate [this] (map negate this))) + + +(defn ordered-like + "Create a comparator using the sorted collection as an + example. Elements not in the sorted collection are sorted to the + end." + [sorted-coll] + (let [ + sorted-coll? (set sorted-coll) + ascending-pair? + (set(reduce concat + (map-indexed + (fn [n x] + (map #(vector x %) (nthnext sorted-coll n))) + sorted-coll)))] + (fn [x y] + (cond + (= x y) 0 + (ascending-pair? [x y]) -1 + (ascending-pair? [y x]) 1 + (sorted-coll? x) -1 + (sorted-coll? y) 1)))) + + + +(def type-precedence + (ordered-like [incanter.Matrix])) + +(defmulti add + (fn [x y] + (sort type-precedence [(type x)(type y)]))) + +(defmulti multiply + (fn [x y] + (sort type-precedence [(type x) (type y)]))) + +(defmethod add [java.lang.Number java.lang.Number] [x y] (+ x y)) +(defmethod multiply [java.lang.Number java.lang.Number] [x y] (* x y)) + +(defmethod multiply [incanter.Matrix java.lang.Integer] [x y] + (let [args (sort #(type-precedence (type %1)(type %2)) [x y]) + matrix (first args) + scalar (second args)] + (incanter.core/matrix (map (partial map (partial multiply scalar)) matrix)))) + +#+end_src + + +* Useful Data Types + +** Complex Numbers +#+srcname: complex-numbers +#+begin_src clojure +(in-ns 'sicm.utils) + +(defprotocol Complex + (real-part [z]) + (imaginary-part [z]) + (magnitude-squared [z]) + (angle [z]) + (conjugate [z]) + (norm [z])) + +(defn complex-rectangular + "Define a complex number with the given real and imaginary + components." + [re im] + (reify Complex + (real-part [z] re) + (imaginary-part [z] im) + (magnitude-squared [z] (+ (* re re) (* im im))) + (angle [z] (java.lang.Math/atan2 im re)) + (conjugate [z] (complex-rectangular re (- im))) + + Arithmetic + (zero [z] (complex-rectangular 0 0)) + (one [z] (complex-rectangular 1 0)) + (negate [z] (complex-rectangular (- re) (- im))) + (invert [z] (complex-rectangular + (/ re (magnitude-squared z)) + (/ (- im) (magnitude-squared z)))) + + Object + (toString [_] + (if (and (zero? re) (zero? im)) (str 0) + (str + (if (not(zero? re)) + re) + (if ((comp not zero?) im) + (str + (if (neg? im) "-" "+") + (if ((comp not one?) (java.lang.Math/abs im)) + (java.lang.Math/abs im)) + "i"))))))) + +(defn complex-polar + "Define a complex number with the given magnitude and angle." + [mag ang] + (reify Complex + (magnitude-squared [z] (* mag mag)) + (angle [z] angle) + (real-part [z] (* mag (java.lang.Math/cos ang))) + (imaginary-part [z] (* mag (java.lang.Math/sin ang))) + (conjugate [z] (complex-polar mag (- ang))) + + Arithmetic + (zero [z] (complex-polar 0 0)) + (one [z] (complex-polar 1 0)) + (negate [z] (complex-polar (- mag) ang)) + (invert [z] (complex-polar (/ mag) (- ang))) + + Object + (toString [_] (str mag " * e^(i" ang")")) + )) + + +;; Numbers are complex quantities + +(extend-protocol Complex + java.lang.Number + (real-part [x] x) + (imaginary-part [x] 0) + (magnitude [x] x) + (angle [x] 0) + (conjugate [x] x)) + + +#+end_src + + + +** Tuples and Tensors + +A tuple is a vector which is spinable\mdash{}it can be either /spin +up/ or /spin down/. (Covariant, contravariant; dual vectors) +#+srcname: tuples +#+begin_src clojure +(in-ns 'sicm.utils) + +(defprotocol Spinning + (up? [this]) + (down? [this])) + +(defn spin + "Returns the spin of the Spinning s, either :up or :down" + [#^Spinning s] + (cond (up? s) :up (down? s) :down)) + + +(deftype Tuple + [spin coll] + clojure.lang.Seqable + (seq [this] (seq (.coll this))) + clojure.lang.Counted + (count [this] (count (.coll this)))) + +(extend-type Tuple + Spinning + (up? [this] (= ::up (.spin this))) + (down? [this] (= ::down (.spin this)))) + +(defmethod print-method Tuple + [o w] + (print-simple (str (if (up? o) 'u 'd) (.coll o)) w)) + + + +(defn up + "Create a new up-tuple containing the contents of coll." + [coll] + (Tuple. ::up coll)) + +(defn down + "Create a new down-tuple containing the contents of coll." + [coll] + (Tuple. ::down coll)) + + +#+end_src + +*** Contraction +Contraction is a binary operation that you can apply to compatible + tuples. Tuples are compatible for contraction if they have the same + length and opposite spins, and if the corresponding items in each + tuple are both numbers or both compatible tuples. + +#+srcname: tuples-2 +#+begin_src clojure +(in-ns 'sicm.utils) + +(defn numbers? + "Returns true if all arguments are numbers, else false." + [& xs] + (every? number? xs)) + +(defn contractible? + "Returns true if the tuples a and b are compatible for contraction, + else false. Tuples are compatible if they have the same number of + components, they have opposite spins, and their elements are + pairwise-compatible." + [a b] + (and + (isa? (type a) Tuple) + (isa? (type b) Tuple) + (= (count a) (count b)) + (not= (spin a) (spin b)) + + (not-any? false? + (map #(or + (numbers? %1 %2) + (contractible? %1 %2)) + a b)))) + + + +(defn contract + "Contracts two tuples, returning the sum of the + products of the corresponding items. Contraction is recursive on + nested tuples." + [a b] + (if (not (contractible? a b)) + (throw + (Exception. "Not compatible for contraction.")) + (reduce + + (map + (fn [x y] + (if (numbers? x y) + (* x y) + (contract x y))) + a b)))) + +#+end_src + +*** Matrices +#+srcname: matrices +#+begin_src clojure +(in-ns 'sicm.utils) +(require 'incanter.core) ;; use incanter's fast matrices + +(defprotocol Matrix + (rows [matrix]) + (cols [matrix]) + (diagonal [matrix]) + (trace [matrix]) + (determinant [matrix]) + (transpose [matrix]) + (conjugate [matrix]) +) + +(extend-protocol Matrix + incanter.Matrix + (rows [rs] (map down (apply map vector (apply map vector rs)))) + (cols [rs] (map up (apply map vector rs))) + (diagonal [matrix] (incanter.core/diag matrix) ) + (determinant [matrix] (incanter.core/det matrix)) + (trace [matrix] (incanter.core/trace matrix)) + (transpose [matrix] (incanter.core/trans matrix)) + ) + +(defn count-rows [matrix] + ((comp count rows) matrix)) + +(defn count-cols [matrix] + ((comp count cols) matrix)) + +(defn square? [matrix] + (= (count-rows matrix) (count-cols matrix))) + +(defn identity-matrix + "Define a square matrix of size n-by-n with 1s along the diagonal and + 0s everywhere else." + [n] + (incanter.core/identity-matrix n)) + + + + + +(defn matrix-by-rows + "Define a matrix by giving its rows." + [& rows] + (if + (not (all-equal? (map count rows))) + (throw (Exception. "All rows in a matrix must have the same number of elements.")) + (incanter.core/matrix (vec rows)))) + +(defn matrix-by-cols + "Define a matrix by giving its columns" + [& cols] + (if (not (all-equal? (map count cols))) + (throw (Exception. "All columns in a matrix must have the same number of elements.")) + (incanter.core/matrix (vec (apply map vector cols))))) + +(defn identity-matrix + "Define a square matrix of size n-by-n with 1s along the diagonal and + 0s everywhere else." + [n] + (incanter.core/identity-matrix n)) + + + +(extend-protocol Arithmetic + incanter.Matrix + (one [matrix] + (if (square? matrix) + (identity-matrix (count-rows matrix)) + (throw (Exception. "Non-square matrices have no multiplicative unit.")))) + (zero [matrix] + (apply matrix-by-rows (map zero (rows matrix)))) + (negate [matrix] + (apply matrix-by-rows (map negate (rows matrix)))) + (invert [matrix] + (incanter.core/solve matrix))) + + + +(defmulti coerce-to-matrix + "Converts x into a matrix, if possible." + type) + +(defmethod coerce-to-matrix incanter.Matrix [x] x) +(defmethod coerce-to-matrix Tuple [x] + (if (apply numbers? (seq x)) + (if (up? x) + (matrix-by-cols (seq x)) + (matrix-by-rows (seq x))) + (throw (Exception. "Non-numerical tuple cannot be converted into a matrix.")))) + + + + + +;; (defn matrix-by-cols +;; "Define a matrix by giving its columns." +;; [& cols] +;; (cond +;; (not (all-equal? (map count cols))) +;; (throw (Exception. "All columns in a matrix must have the same number of elements.")) +;; :else +;; (reify Matrix +;; (cols [this] (map up cols)) +;; (rows [this] (map down (apply map vector cols))) +;; (diagonal [this] (map-indexed (fn [i col] (nth col i) cols))) +;; (trace [this] +;; (if (not= (count-cols this) (count-rows this)) +;; (throw (Exception. +;; "Cannot take the trace of a non-square matrix.")) +;; (reduce + (diagonal this)))) + +;; (determinant [this] +;; (if (not= (count-cols this) (count-rows this)) +;; (throw (Exception. +;; "Cannot take the determinant of a non-square matrix.")) +;; (reduce * (map-indexed (fn [i col] (nth col i)) cols)))) +;; ))) + +(extend-protocol Matrix Tuple + (rows [this] (if (down? this) + (list this) + (map (comp up vector) this))) + + (cols [this] (if (up? this) + (list this) + (map (comp down vector) this)) + )) + +(defn matrix-multiply + "Returns the matrix resulting from the matrix multiplication of the given arguments." + ([A] (coerce-to-matrix A)) + ([A B] (incanter.core/mmult (coerce-to-matrix A) (coerce-to-matrix B))) + ([M1 M2 & Ms] (reduce matrix-multiply (matrix-multiply M1 M2) Ms))) + +#+end_src + + +** Power Series +#+srcname power-series +#+begin_src clojure +(in-ns 'sicm.utils) +(use 'clojure.contrib.def) + + + +(defn series-fn + "The function corresponding to the given power series." + [series] + (fn [x] + (reduce + + (map-indexed (fn[n x] (* (float (nth series n)) (float(java.lang.Math/pow (float x) n)) )) + (range 20))))) + +(deftype PowerSeries + [coll] + clojure.lang.Seqable + (seq [this] (seq (.coll this))) + + clojure.lang.Indexed + (nth [this n] (nth (.coll this) n 0)) + (nth [this n not-found] (nth (.coll this) n not-found)) + + ;; clojure.lang.IFn + ;; (call [this] (throw(Exception.))) + ;; (invoke [this & args] args + ;; (let [f + ;; ) + ;; (run [this] (throw(Exception.))) + ) + +(defn power-series + "Returns a power series with the items of the coll as its + coefficients. Trailing zeros are added to the end of coll." + [coeffs] + (PowerSeries. coeffs)) + +(defn power-series-indexed + "Returns a power series consisting of the result of mapping f to the non-negative integers." + [f] + (PowerSeries. (map f (range)))) + + +(defn-memo nth-partial-sum + ([series n] + (if (zero? n) (first series) + (+ (nth series n) + (nth-partial-sum series (dec n)))))) + +(defn partial-sums [series] + (lazy-seq (map nth-partial-sum (range)))) + + + + +(def cos-series + (power-series-indexed + (fn[n] + (if (odd? n) 0 + (/ + (reduce * + (reduce * (repeat (/ n 2) -1)) + (range 1 (inc n))) + ))))) + +(def sin-series + (power-series-indexed + (fn[n] + (if (even? n) 0 + (/ + (reduce * + (reduce * (repeat (/ (dec n) 2) -1)) + (range 1 (inc n))) + ))))) + +#+end_src + + +* Basic Utilities + +** Sequence manipulation + +#+srcname: seq-manipulation +#+begin_src clojure +(ns sicm.utils) + +(defn do-up + "Apply f to each number from low to high, presumably for + side-effects." + [f low high] + (doseq [i (range low high)] (f i))) + +(defn do-down + "Apply f to each number from high to low, presumably for + side-effects." + [f high low] + (doseq [i (range high low -1)] (f i))) + + +(defn all-equal? [coll] + (if (empty? (rest coll)) true + (and (= (first coll) (second coll)) + (recur (rest coll)))))) + +(defn multiplier + "Returns a function that 'multiplies' the members of a collection, +returning unit if called on an empty collection." + [multiply unit] + (fn [coll] ((partial reduce multiply unit) coll))) + +(defn divider + "Returns a function that 'divides' the first element of a collection +by the 'product' of the rest of the collection." + [divide multiply invert unit] + (fn [coll] + (apply + (fn + ([] unit) + ([x] (invert x)) + ([x y] (divide x y)) + ([x y & zs] (divide x (reduce multiply y zs)))) + coll))) + +(defn left-circular-shift + "Remove the first element of coll, adding it to the end of coll." + [coll] + (concat (rest coll) (take 1 coll))) + +(defn right-circular-shift + "Remove the last element of coll, adding it to the front of coll." + [coll] + (cons (last coll) (butlast coll))) +#+end_src + + + + +** Ranges, Airity and Function Composition +#+srcname: arity +#+begin_src clojure +(in-ns 'sicm.utils) +(def infinity Double/POSITIVE_INFINITY) +(defn infinite? [x] (Double/isInfinite x)) +(def finite? (comp not infinite?)) + +(defn arity-min + "Returns the smallest number of arguments f can take." + [f] + (apply + min + (map (comp alength #(.getParameterTypes %)) + (filter (comp (partial = "invoke") #(.getName %)) + (.getDeclaredMethods (class f)))))) + +(defn arity-max + "Returns the largest number of arguments f can take, possibly + Infinity." + [f] + (let [methods (.getDeclaredMethods (class f))] + (if (not-any? (partial = "doInvoke") (map #(.getName %) methods)) + (apply max + (map (comp alength #(.getParameterTypes %)) + (filter (comp (partial = "invoke") #(.getName %)) methods))) + infinity))) + + +(def ^{:arglists '([f]) + :doc "Returns a two-element list containing the minimum and + maximum number of args that f can take."} + arity-interval + (juxt arity-min arity-max)) + + + +;; --- intervals + +(defn intersect + "Returns the interval of overlap between interval-1 and interval-2" + [interval-1 interval-2] + (if (or (empty? interval-1) (empty? interval-2)) [] + (let [left (max (first interval-1) (first interval-2)) + right (min (second interval-1) (second interval-2))] + (if (> left right) [] + [left right])))) + +(defn same-endpoints? + "Returns true if the left endpoint is the same as the right + endpoint." + [interval] + (= (first interval) (second interval))) + +(defn naturals? + "Returns true if the left endpoint is 0 and the right endpoint is +infinite." + [interval] + (and (zero? (first interval)) + (infinite? (second interval)))) + + +(defn fan-in + "Returns a function that pipes its input to each of the gs, then + applies f to the list of results. Consequently, f must be able to + take a number of arguments equal to the number of gs." + [f & gs] + (fn [& args] + (apply f (apply (apply juxt gs) args)))) + +(defn fan-in + "Returns a function that pipes its input to each of the gs, then applies f to the list of results. The resulting function takes any number of arguments, but will fail if given arguments that are incompatible with any of the gs." + [f & gs] + (comp (partial apply f) (apply juxt gs))) + + + +(defmacro airty-blah-sad [f n more?] + (let [syms (vec (map (comp gensym (partial str "x")) (range n))) + optional (gensym "xs")] + (if more? + `(fn ~(conj syms '& optional) + (apply ~f ~@syms ~optional)) + `(fn ~syms (~f ~@syms))))) + +(defmacro airt-whaa* [f n more?] + `(airty-blah-sad ~f ~n ~more?)) + + + + +(defn fan-in* + "Returns a function that pipes its input to each of the gs, then + applies f to the list of results. Unlike fan-in, fan-in* strictly + enforces arity: it will fail if the gs do not have compatible + arities." + [f & gs] + (let [arity-in (reduce intersect (map arity-interval gs)) + left (first arity-in) + right (second arity-in) + composite (fan-in f gs) + ] + (cond + (empty? arity-in) + (throw (Exception. "Cannot compose functions with incompatible arities.")) + + (not + (or (= left right) + (and (finite? left) + (= right infinity)))) + + (throw (Exception. + "Compose can only handle arities of the form [n n] or [n infinity]")) + :else + (airty-blah-sad composite left (= right infinity))))) + + + +(defn compose-n "Compose any number of functions together." + ([] identity) + ([f] f) + ([f & fs] + (let [fns (cons f fs)] + (compose-bin (reduce fan-in (butlast fs)) (last fs)))) +) + + + + + + +(defn iterated + ([f n id] (reduce comp id (repeat n f))) + ([f n] (reduce comp identity (repeat n f)))) + +(defn iterate-until-stable + "Repeatedly applies f to x, returning the first result that is close +enough to its predecessor." + [f close-enough? x] + (second (swank.util/find-first + (partial apply close-enough?) + (partition 2 1 (iterate f x))))) + +(defn lexical< [x y] + (neg? (compare (str x) (str y)))) + + +;; do-up +;; do-down +(def make-pairwise-test comparator) +;;all-equal? +(def accumulation multiplier) +(def inverse-accumulation divider) +;;left-circular-shift +;;right-circular-shift +(def exactly-n? same-endpoints?) +(def any-number? naturals?) +;; TODO compose +;; TODO compose-n +;; identity +(def compose-2 fan-in) +(def compose-bin fan-in*) +(def any? (constantly true)) +(def none? (constantly false)) +(def constant constantly) +(def joint-arity intersect) +(def a-reduce reduce) +;; filter +(def make-map (partial partial map) ) +(def bracket juxt) +;; TODO apply-to-all +;; TODO nary-combine +;; TODO binary-combine +;; TODO unary-combine +;; iterated +;; iterate-until-stable +(def make-function-of-vector (partial partial map)) +(def make-function-of-arguments (fn [f] (fn [& args] (f args)))) +(def alphaless lexical<) + +#+end_src + + + + + + + +: + +* Numerical Methods +#+srcname: numerical-methods +#+begin_src clojure +(in-ns 'sicm.utils) +(import java.lang.Math) +(use 'clojure.contrib.def) + +;; ---- USEFUL CONSTANTS + +(defn machine-epsilon + "Smallest value representable on your machine, as determined by +successively dividing a number in half until consecutive results are +indistinguishable." + [] + (ffirst + (drop-while + (comp not zero? second) + (partition 2 1 + (iterate (partial * 0.5) 1))))) + + +(def pi (Math/PI)) +(def two-pi (* 2 pi)) + +(def eulers-gamma 0.5772156649015328606065) + +(def phi (/ (inc (Math/sqrt 5)) 2)) + +(def ln2 (Math/log 2)) +(def ln10 (Math/log 10)) +(def exp10 #(Math/pow 10 %)) +(def exp2 #(Math/pow 2 %)) + + +;; + +;; ---- ANGLES AND TRIGONOMETRY + +(defn angle-restrictor + "Returns a function that ensures that angles lie in the specified interval of length two-pi." + [max-angle] + (let [min-angle (- max-angle two-pi)] + (fn [x] + (if (and + (<= min-angle x) + (< x max-angle)) + x + (let [corrected-x (- x (* two-pi (Math/floor (/ x two-pi))))] + (if (< corrected-x max-angle) + corrected-x + (- corrected-x two-pi))))))) + +(defn angle-restrict-pi + "Coerces angles to lie in the interval from -pi to pi." + [angle] + ((angle-restrictor pi) angle)) + +(defn angle-restrict-two-pi + "Coerces angles to lie in the interval from zero to two-pi" + [angle] + ((angle-restrictor two-pi) angle)) + + + + +(defn invert [x] (/ x)) +(defn negate [x] (- x)) + +(defn exp [x] (Math/exp x)) + +(defn sin [x] (Math/sin x)) +(defn cos [x] (Math/cos x)) +(defn tan [x] (Math/tan x)) + +(def sec (comp invert cos)) +(def csc (comp invert sin)) + +(defn sinh [x] (Math/sinh x)) +(defn cosh [x] (Math/cosh x)) +(defn tanh [x] (Math/tanh x)) + +(def sech (comp invert cosh)) +(def csch (comp invert sinh)) + + +;; ------------ + +(defn factorial + "Computes the factorial of the nonnegative integer n." + [n] + (if (neg? n) + (throw (Exception. "Cannot compute the factorial of a negative number.")) + (reduce * 1 (range 1 (inc n))))) + +(defn exact-quotient [n d] (/ n d)) + +(defn binomial-coefficient + "Computes the number of different ways to choose m elements from n." + [n m] + (assert (<= 0 m n)) + (let [difference (- n m)] + (exact-quotient + (reduce * (range n (max difference m) -1 )) + (factorial (min difference m))))) + +(defn-memo stirling-1 + "Stirling numbers of the first kind: the number of permutations of n + elements with exactly m permutation cycles. " + [n k] + ;(assert (<= 1 k n)) + (if (zero? n) + (if (zero? k) 1 0) + (+ (stirling-1 (dec n) (dec k)) + (* (dec n) (stirling-1 (dec n) k))))) + +(defn-memo stirling-2 ;;count-partitions + "Stirling numbers of the second kind: the number of ways to partition a set of n elements into k subsets." + [n k] + (cond + (= k 1) 1 + (= k n) 1 + :else (+ (stirling-2 (dec n) (dec k)) + (* k (stirling-2 (dec n) k))))) + +(defn harmonic-number [n] + (/ (stirling-1 (inc n) 2) + (factorial n))) + + +(defn sum + [f low high] + (reduce + (map f (range low (inc high))))) + +#+end_src + + + + + + + + + + + + +* Differentiation + +We compute derivatives by passing special *differential objects* $[x, +dx]$ through functions. Roughly speaking, applying a function $f$ to a +differential object $[x, dx]$ should produce a new differential +object $[f(x),\,Df(x)\cdot dx]$. + +$[x,\,dx]\xrightarrow{\quad f \quad}[f(x),\,Df(x)\cdot dx]$ +Notice that you can obtain the derivative of $f$ from this +differential object, as it is the coefficient of the $dx$ term. Also, +as you apply successive functions using this rule, you get the +chain-rule answer you expect: + +$[f(x),\,Df(x)\cdot dx]\xrightarrow{\quad g\quad} [gf(x),\, +Dgf(x)\cdot Df(x) \cdot dx ]$ + +In order to generalize to multiple variables and multiple derivatives, +we use a *power series of differentials*, a sortred infinite sequence which +contains all terms like $dx\cdot dy$, $dx^2\cdot dy$, etc. + + +#+srcname:differential +#+begin_src clojure +(in-ns 'sicm.utils) +(use 'clojure.contrib.combinatorics) +(use 'clojure.contrib.generic.arithmetic) + +(defprotocol DifferentialTerm + "Protocol for an infinitesimal quantity." + (coefficient [this]) + (partials [this])) + +(extend-protocol DifferentialTerm + java.lang.Number + (coefficient [this] this) + (partials [this] [])) + +(deftype DifferentialSeq + [terms] + clojure.lang.IPersistentCollection + (cons [this x] + (throw (Exception. x)) + (DifferentialSeq. (cons x terms))) + ;; (seq [this] (seq terms)) + (count [this] (count terms)) + (empty [this] (empty? terms))) + + + + + + + + + + + +(defn coerce-to-differential-seq [x] + (cond + (= (type x) DifferentialSeq) x + (satisfies? DifferentialTerm x) (DifferentialSeq. x))) + + +(defn differential-term + "Returns a differential term with the given coefficient and + partials. Coefficient may be any arithmetic object; partials must + be a list of non-negative integers." + [coefficient partials] + (reify DifferentialTerm + (partials [_] (set partials)) + (coefficient [_] coefficient))) + + + +(defn differential-seq* + ([coefficient partials] + (DifferentialSeq. [(differential-term coefficient partials)])) + ([coefficient partials & cps] + (if cps + + + + +(defn differential-seq + "Constructs a sequence of differential terms from a numerical +coefficient and a list of keys for variables. If no coefficient is supplied, uses 1." + ([variables] (differential-seq 1 variables)) + ([coefficient variables & cvs] + (if (number? coefficient) + (conj (assoc {} (apply sorted-set variables) coefficient) + (if (empty? cvs) + nil + (apply differential-seq cvs))) + (apply differential-seq 1 coefficient 1 variables cvs) + ))) + + +(defn differential-add + "Add two differential sequences by combining like terms." + [dseq1 dseq2] + (merge-with + dseq1 dseq2)) + +(defn differential-multiply + "Multiply two differential sequences. The square of any differential variable is zero since differential variables are infinitesimally small." + [dseq1 dseq2] + (reduce + (fn [m [[vars1 coeff1] [vars2 coeff2]]] + (if (empty? (clojure.set/intersection vars1 vars2)) + (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2)) + m)) + {} + (clojure.contrib.combinatorics/cartesian-product + dseq1 + dseq2))) + + +(defn big-part + "Returns the part of the differential sequence that is finite, + i.e. not infinitely small." + [dseq] + (let + [keys (sort-by count (keys dseq)) + smallest-var (last(last keys))] + (apply hash-map + (reduce concat + (remove (comp smallest-var first) dseq))))) + +(defn small-part + "Returns the part of the differential sequence that is + infinitesimal." + [dseq] + (let + [keys (sort-by count (keys dseq)) + smallest-var (last(last keys))] + (apply hash-map + (reduce concat + (filter (comp smallest-var first) dseq))))) + +(defn big-part + "Returns the 'finite' part of the differential sequence." + [dseq] + (let + [keys (sort-by count (keys dseq)) + smallest-var (last (last keys))] + (apply hash-map + (reduce concat + (filter (remove smallest-var first) dseq) + )))) + + +(defn small-part + "Returns the 'infinitesimal' part of the differential sequence." + [dseq] + (let + [keys (sort-by count (keys dseq)) + smallest-var (last (last keys))] + + (apply hash-map + (reduce concat + (filter (comp smallest-var first) dseq) + )))) + +(defn linear-approximator + "Returns the function that corresponds to unary-fn but which accepts and returns differential objects." + ([unary-f dfdx] + (fn F [dseq] + (differential-add + (unary-f (big-part dseq)) + (differential-multiply + (dfdx (big-part dseq)) + (small-part dseq)))))) + + + + + + +;; ;; A differential term consists of a numerical coefficient and a +;; ;; sorted +;; (defrecord DifferentialTerm [coefficient variables]) +;; (defmethod print-method DifferentialTerm +;; [o w] +;; (print-simple +;; (apply str (.coefficient o)(map (comp (partial str "d") name) (.variables o))) +;; w)) + + +;; (defn differential-seq +;; "Constructs a sequence of differential terms from a numerical +;; coefficient and a list of keywords for variables. If no coefficient is +;; supplied, uses 1." +;; ([variables] (differential-seq 1 variables)) +;; ([coefficient variables] +;; (list +;; (DifferentialTerm. coefficient (apply sorted-set variables)))) +;; ([coefficient variables & cvs] +;; (sort-by +;; #(vec(.variables %)) +;; (concat (differential-seq coefficient variables) (apply differential-seq cvs))))) + +#+end_src + + + + + + + + + + + + + +* Symbolic Manipulation + +#+srcname:symbolic +#+begin_src clojure +(in-ns 'sicm.utils) + +(deftype Symbolic [type expression] + Object + (equals [this that] + (cond + (= (.expression this) (.expression that)) true + :else + (Symbolic. + java.lang.Boolean + (list '= (.expression this) (.expression that))) + ))) + + + + +(deftype AbstractSet [glyph membership-test]) +(defn member? [abstract-set x] + ((.membership-test abstract-set) x)) + +;; ------------ Some important AbstractSets + + +(def Real + (AbstractSet. + 'R + (fn[x](number? x)))) + + +;; ------------ Create new AbstractSets from existing ones + +(defn abstract-product + "Gives the cartesian product of abstract sets." + ([sets] + (if (= (count sets) 1) (first sets) + (AbstractSet. + (symbol + (apply str + (interpose 'x (map #(.glyph %) sets)))) + (fn [x] + (and + (coll? x) + (= (count sets) (count x)) + (reduce #(and %1 %2) + true + (map #(member? %1 %2) sets x))))))) + ([abstract-set n] + (abstract-product (repeat n abstract-set)))) + + + +(defn abstract-up + "Returns the abstract set of all up-tuples whose items belong to the + corresponding abstract sets in coll." + ([coll] + (AbstractSet. + (symbol (str "u[" + (apply str + (interpose " " + (map #(.glyph %) coll))) + "]")) + (fn [x] + (and + (satisfies? Spinning x) + (up? x) + (= (count coll) (count x)) + (reduce + #(and %1 %2) + true + (map #(member? %1 %2) coll x)))))) + ([abstract-set n] + (abstract-up (repeat n abstract-set)))) + + +(defn abstract-down + "Returns the abstract set of all down-tuples whose items belong to the + corresponding abstract sets in coll." + ([coll] + (AbstractSet. + (symbol (str "d[" + (apply str + (interpose " " + (map #(.glyph %) coll))) + "]")) + (fn [x] + (and + (satisfies? Spinning x) + (down? x) + (= (count coll) (count x)) + (reduce + #(and %1 %2) + true + (map #(member? %1 %2) coll x)))))) + ([abstract-set n] + (abstract-down (repeat n abstract-set)))) + + + + + + ;-------ABSTRACT FUNCTIONS +(defrecord AbstractFn + [#^AbstractSet domain #^AbstractSet codomain]) + + +(defmethod print-method AbstractFn + [o w] + (print-simple + (str + "f:" + (.glyph (:domain o)) + "-->" + (.glyph (:codomain o))) w)) +#+end_src + + +* COMMENT +#+begin_src clojure :tangle utils.clj +(ns sicm.utils) + + ;***** GENERIC ARITHMETIC +<> + + + ;***** TUPLES AND MATRICES +<> +<> + ;***** MATRICES +<> + +#+end_src + + + diff -r 8d8278e09888 -r b4de894a1e2e sicm/deriv.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/deriv.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,240 @@ + + + + +An Unambiguous Notation for Derivatives + + + + + + + + + + + + + +

+ + + + + +

An Unambiguous Notation for Derivatives

+ + + + + + + +

1 Calculus of Infinitesimals +
- 1.1 Differential Objects
- 1.2 Interactions obey the chain rule
+

+ +

1 Calculus of Infinitesimals

+ + +

+ +

1.1 Differential Objects

+ + +

+A differential object is a pair $[x,\,dx]$ consisting of a variable +and an infinitely small increment of it. Differential objects can +interact with functions, producing a new differential object as a +result; this interaction is for calculating derivatives of functions. +

+Differential objects are for +calculating derivatives of functions: the derivative of $f$ with +respect to $x$ +

+You can “apply” +functions to differential objects; the result is: +

+ + +$[x,dx]\xrightarrow{\quad f \quad}[f(x), Df(x)\cdot dx].$ + +

+Loosely speaking, the interaction of $f$ and a differential object +of $x$ is a differential object of $f$. +

+ +

1.2 Interactions obey the chain rule

+ + +

+The interaction of $f$ and the differential object $[x, dx]$ is +a differential object $[f(x), Df(x)\cdot dx]$. Because of the rule for +interactions, if you apply another function $g$, you get the +chain-rule answer you expect: +

+ + +$[f(x), Df(x)\cdot dx]\xrightarrow{\quad g\quad}\left[g(f(x)),\, +Dg(f(x))\cdot Df(x)\cdot dx\right]$ + + + + + + + + +

Date: 2011-08-08 02:49:24 EDT

Author: Dylan Holmes

Org version 7.6 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e sicm/deriv.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/deriv.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,53 @@ +#+TITLE:An Unambiguous Notation for Derivatives +#+author: Dylan Holmes +#+EMAIL: rlm@mit.edu +#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js" +#+STYLE: +#+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t +#+SETUPFILE: ../templates/level-0.org +#+INCLUDE: ../templates/level-0.org +#+BABEL: :noweb yes + +* Calculus of Infinitesimals +** Differential Objects + +A *differential object* is a pair $[x,\,dx]$ consisting of a variable +and an infinitely small increment of it. We want differential objects +to enable us to compute derivatives of functions. + +Differential objects are for +calculating derivatives of functions: the derivative of $f$ with +respect to $x$ + +You can \ldquo{}apply\rdquo{} +functions to differential objects; the result is: + +$[x,dx]\xrightarrow{\quad f \quad}[f(x), Df(x)\cdot dx].$ + +Loosely speaking, the interaction of $f$ and a differential object +of $x$ is a differential object of $f$. + +#As a linguistic convention, we'll call this interaction /applying f +#to the differential object/. This is not to be confused with the +#=apply= function in Clojure. + +** Interactions obey the chain rule + +The interaction of $f$ and the differential object $[x, dx]$ is +a differential object $[f(x), Df(x)\cdot dx]$. Because of the rule for +interactions, if you apply another function $g$, you get the +chain-rule answer you expect: + +$[f(x), Df(x)\cdot dx]\xrightarrow{\quad g\quad}\left[g(f(x)),\, +Dg(f(x))\cdot Df(x)\cdot dx\right]$ + + +#+begin_src clojure :tangle deriv.clj + +#+end_src + +#+results: +: nil + + + diff -r 8d8278e09888 -r b4de894a1e2e sicm/notes.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/notes.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,141 @@ + + + + +Notes for SICM + + + + + + + + + + + + +

+ + + +

1 Tuples
2 Generic arithmetic

+ +

1 Tuples

+ +

Tuples are a new data type: sequences with spin. A tuple can be either spin-up or spin-down. +
A pair of compatible tuples can be contracted into a single number. +
- Tuples are compatible if they have the same length and opposite + spin, and if their corresponding pairs of items are either both + numbers or both compatible tuples. +
- To contract tuples, take the sum of the products of corresponding + pairs of items. (To take the product of compatible tuples, + contract them.) +
+ +
+

+ + +

+ +

2 Generic arithmetic

+ + +

Date: 2011-08-10 14:07:53 EDT

Author: Robert McIntyre

Org version 7.6 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e sicm/notes.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/notes.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,25 @@ +#+TITLE: Notes for SICM + +* Complex numbers + + +* Multidimensional Structures +** Tuples +- Tuples are a new data type: sequences with /spin/. A tuple can be either spin-up or spin-down. +- A pair of compatible tuples can be *contracted* into a single number. + - Tuples are compatible if they have the same length and opposite + spin, and if their corresponding pairs of items are either both + numbers or both compatible tuples. + - To contract tuples, take the sum of the products of corresponding + pairs of items. (To take the product of compatible tuples, + contract them.) + +*Generic Arithmetic:* To multiply two tuples, check to see if they are + compatible. If they are, contract them. If they aren't, multiply each + item of the second tuple by the first tuple. + +To multiply a tuple by a number, multiply each element of the tuple by +the number. + +** Matrices +Uses incanter matrices. diff -r 8d8278e09888 -r b4de894a1e2e sicm/utils.clj --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/utils.clj Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,496 @@ + +(ns sicm.utils) + +(in-ns 'sicm.utils) + +;; Let some objects have spin + +(defprotocol Spinning + (up? [this]) + (down? [this])) + +(defn spin + "Returns the spin of the Spinning s, either :up or :down" + [#^Spinning s] + (cond (up? s) :up (down? s) :down)) + + +;; DEFINITION: A tuple is a sequence with spin + +(deftype Tuple + [spin coll] + + clojure.lang.Seqable + (seq [this] (seq (.coll this))) + + clojure.lang.Counted + (count [this] (count (.coll this))) + + Spinning + (up? [this] (= ::up (.spin this))) + (down? [this] (= ::down (.spin this)))) + +(defmethod print-method Tuple + [o w] + (print-simple + (if (up? o) + (str "u" (.coll o)) + (str "d" (vec(.coll o)))) + w)) + +(def tuple? #(= (type %) Tuple)) + +;; CONSTRUCTORS + +(defn up + "Create a new up-tuple containing the contents of coll." + [coll] + (Tuple. ::up coll)) + +(defn down + "Create a new down-tuple containing the contents of coll." + [coll] + (Tuple. ::down coll)) + +(defn same-spin + "Creates a tuple which has the same spin as tuple and which contains +the contents of coll." + [tuple coll] + (if (up? tuple) + (up coll) + (down coll))) + +(defn opposite-spin + "Create a tuple which has opposite spin to tuple and which contains +the contents of coll." + [tuple coll] + (if (up? tuple) + (down coll) + (up coll))) +(in-ns 'sicm.utils) +(require 'incanter.core) ;; use incanter's fast matrices + + +(defn all-equal? [coll] + (if (empty? (rest coll)) true + (and (= (first coll) (second coll)) + (recur (rest coll))))) + + +(defprotocol Matrix + (rows [matrix]) + (cols [matrix]) + (diagonal [matrix]) + (trace [matrix]) + (determinant [matrix]) + (transpose [matrix]) + (conjugate [matrix]) +) + +(extend-protocol Matrix incanter.Matrix + (rows [rs] (map down (apply map vector (apply map vector rs)))) + (cols [rs] (map up (apply map vector rs))) + (diagonal [matrix] (incanter.core/diag matrix) ) + (determinant [matrix] (incanter.core/det matrix)) + (trace [matrix] (incanter.core/trace matrix)) + (transpose [matrix] (incanter.core/trans matrix))) + +(defn count-rows [matrix] + ((comp count rows) matrix)) + +(defn count-cols [matrix] + ((comp count cols) matrix)) + +(defn square? [matrix] + (= (count-rows matrix) (count-cols matrix))) + +(defn identity-matrix + "Define a square matrix of size n-by-n with 1s along the diagonal and + 0s everywhere else." + [n] + (incanter.core/identity-matrix n)) + + +(defn matrix-by-rows + "Define a matrix by giving its rows." + [& rows] + (if + (not (all-equal? (map count rows))) + (throw (Exception. "All rows in a matrix must have the same number of elements.")) + (incanter.core/matrix (vec rows)))) + +(defn matrix-by-cols + "Define a matrix by giving its columns" + [& cols] + (if (not (all-equal? (map count cols))) + (throw (Exception. "All columns in a matrix must have the same number of elements.")) + (incanter.core/matrix (vec (apply map vector cols))))) + +(defn identity-matrix + "Define a square matrix of size n-by-n with 1s along the diagonal and + 0s everywhere else." + [n] + (incanter.core/identity-matrix n)) + +(in-ns 'sicm.utils) +(use 'clojure.contrib.generic.arithmetic + 'clojure.contrib.generic.collection + 'clojure.contrib.generic.functor + 'clojure.contrib.generic.math-functions) + +(defn numbers? + "Returns true if all arguments are numbers, else false." + [& xs] + (every? number? xs)) + +(defn tuple-surgery + "Applies the function f to the items of tuple and the additional + arguments, if any. Returns a Tuple of the same type as tuple." + [tuple f & xs] + ((if (up? tuple) up down) + (apply f (seq tuple) xs))) + + + +;;; CONTRACTION collapses two compatible tuples into a number. + +(defn contractible? + "Returns true if the tuples a and b are compatible for contraction, + else false. Tuples are compatible if they have the same number of + components, they have opposite spins, and their elements are + pairwise-compatible." + [a b] + (and + (isa? (type a) Tuple) + (isa? (type b) Tuple) + (= (count a) (count b)) + (not= (spin a) (spin b)) + + (not-any? false? + (map #(or + (numbers? %1 %2) + (contractible? %1 %2)) + a b)))) + +(defn contract + "Contracts two tuples, returning the sum of the + products of the corresponding items. Contraction is recursive on + nested tuples." + [a b] + (if (not (contractible? a b)) + (throw + (Exception. "Not compatible for contraction.")) + (reduce + + (map + (fn [x y] + (if (numbers? x y) + (* x y) + (contract x y))) + a b)))) + + + + + +(defmethod conj Tuple + [tuple & xs] + (tuple-surgery tuple #(apply conj % xs))) + +(defmethod fmap Tuple + [f tuple] + (tuple-surgery tuple (partial map f))) + + + +;; TODO: define Scalar, and add it to the hierarchy above Number and Complex + + +(defmethod * [Tuple Tuple] ; tuple*tuple + [a b] + (if (contractible? a b) + (contract a b) + (map (partial * a) b))) + + +(defmethod * [java.lang.Number Tuple] ;; scalar * tuple + [a x] (fmap (partial * a) x)) + +(defmethod * [Tuple java.lang.Number] + [x a] (* a x)) + +(defmethod * [java.lang.Number incanter.Matrix] ;; scalar * matrix + [x M] (incanter.core/mult x M)) + +(defmethod * [incanter.Matrix java.lang.Number] + [M x] (* x M)) + +(defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix + [M1 M2] + (incanter.core/mmult M1 M2)) + +(defmethod * [incanter.Matrix Tuple] ;; matrix * tuple + [M v] + (if (and (apply numbers? v) (up? v)) + (* M (matrix-by-cols v)) + (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*")) + )) + +(defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix + [v M] + (if (and (apply numbers? v) (down? v)) + (* (matrix-by-rows v) M) + (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*")) + )) + + +(defmethod exp incanter.Matrix + [M] + (incanter.core/exp M)) + + +(in-ns 'sicm.utils) +(use 'clojure.contrib.seq + 'clojure.contrib.generic.arithmetic + 'clojure.contrib.generic.collection + 'clojure.contrib.generic.math-functions) + +;;∂ + +;; DEFINITION : Differential Term + +;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The +;; coefficient of the quantity is returned by the 'coefficient' method, +;; while the sequence of differential parameters is returned by the +;; method 'partials'. + +;; Instead of using (potentially ambiguous) letters to denote +;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2]. + +;; The coefficient can be any arithmetic object; the +;; partials must be a nonrepeating sorted sequence of nonnegative +;; integers. + +(deftype DifferentialTerm [coefficient partials]) + +(defn differential-term + "Make a differential term from a coefficient and list of partials." + [coefficient partials] + (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials)) + (DifferentialTerm. coefficient (set partials)) + (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers.")))) + + +;; DEFINITION : Differential Sequence +;; A differential sequence is a sequence of differential terms, all with different partials. +;; Internally, it is a map from the partials of each term to their coefficients. + +(deftype DifferentialSeq + [terms] + ;;clojure.lang.IPersistentMap + clojure.lang.Associative + (assoc [this key val] + (DifferentialSeq. + (cons (differential-term val key) terms))) + (cons [this x] + (DifferentialSeq. (cons x terms))) + (containsKey [this key] + (not(nil? (find-first #(= (.partials %) key) terms)))) + (count [this] (count (.terms this))) + (empty [this] (DifferentialSeq. [])) + (entryAt [this key] + ((juxt #(.partials %) #(.coefficient %)) + (find-first #(= (.partials %) key) terms))) + (seq [this] (seq (.terms this)))) + +(def differential? #(= (type %) DifferentialSeq)) + +(defn zeroth-order? + "Returns true if the differential sequence has at most a constant term." + [dseq] + (and + (differential? dseq) + (every? + #(= #{} %) + (keys (.terms dseq))))) + +(defmethod fmap DifferentialSeq + [f dseq] + (DifferentialSeq. + (fmap f (.terms dseq)))) + + + + +;; BUILDING DIFFERENTIAL OBJECTS + +(defn differential-seq + "Define a differential sequence by specifying an alternating +sequence of coefficients and lists of partials." + ([coefficient partials] + (DifferentialSeq. {(set partials) coefficient})) + ([coefficient partials & cps] + (if (odd? (count cps)) + (throw (Exception. "differential-seq requires an even number of terms.")) + (DifferentialSeq. + (reduce + #(assoc %1 (set (second %2)) (first %2)) + {(set partials) coefficient} + (partition 2 cps)))))) + + + +(defn big-part + "Returns the part of the differential sequence that is finite, + i.e. not infinitely small. If the sequence is zeroth-order, returns + the coefficient of the zeroth-order term instead. " + [dseq] + (if (zeroth-order? dseq) (get (.terms dseq) #{}) + (let [m (.terms dseq) + keys (sort-by count (keys m)) + smallest-var (last (last keys))] + (DifferentialSeq. + (reduce + #(assoc %1 (first %2) (second %2)) + {} + (remove #((first %) smallest-var) m)))))) + + +(defn small-part + "Returns the part of the differential sequence that infinitely + small. If the sequence is zeroth-order, returns zero." + [dseq] + (if (zeroth-order? dseq) 0 + (let [m (.terms dseq) + keys (sort-by count (keys m)) + smallest-var (last (last keys))] + (DifferentialSeq. + (reduce + #(assoc %1 (first %2) (second %2)) {} + (filter #((first %) smallest-var) m)))))) + + + +(defn cartesian-product [set1 set2] + (reduce concat + (for [x set1] + (for [y set2] + [x y])))) + +(defn nth-subset [n] + (if (zero? n) [] + (let [lg2 #(/ (log %) (log 2)) + k (int(java.lang.Math/floor (lg2 n))) + ] + (cons k + (nth-subset (- n (pow 2 k))))))) + +(def all-partials + (lazy-seq (map nth-subset (range)))) + + +(defn differential-multiply + "Multiply two differential sequences. The square of any differential + variable is zero since differential variables are infinitesimally + small." + [dseq1 dseq2] + (DifferentialSeq. + (reduce + (fn [m [[vars1 coeff1] [vars2 coeff2]]] + (if (not (empty? (clojure.set/intersection vars1 vars2))) + m + (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2)))) + {} + (cartesian-product (.terms dseq1) (.terms dseq2))))) + + + +(defmethod * [DifferentialSeq DifferentialSeq] + [dseq1 dseq2] + (differential-multiply dseq1 dseq2)) + +(defmethod + [DifferentialSeq DifferentialSeq] + [dseq1 dseq2] + (DifferentialSeq. + (merge-with + (.terms dseq1) (.terms dseq2)))) + +(defmethod * [java.lang.Number DifferentialSeq] + [x dseq] + (fmap (partial * x) dseq)) + +(defmethod * [DifferentialSeq java.lang.Number] + [dseq x] + (fmap (partial * x) dseq)) + +(defmethod + [java.lang.Number DifferentialSeq] + [x dseq] + (+ (differential-seq x []) dseq)) +(defmethod + [DifferentialSeq java.lang.Number] + [dseq x] + (+ dseq (differential-seq x []))) + +(defmethod - DifferentialSeq + [x] + (fmap - x)) + + +;; DERIVATIVES + + + +(defn linear-approximator + "Returns an operator that linearly approximates the given function." + ([f df|dx] + (fn [x] + (let [big-part (big-part x) + small-part (small-part x)] + ;; f(x+dx) ~= f(x) + f'(x)dx + (+ (f big-part) + (* (df|dx big-part) small-part) + )))) + + ([f df|dx df|dy] + (fn [x y] + (let [X (big-part x) + Y (big-part y) + DX (small-part x) + DY (small-part y)] + (+ (f X Y) + (* DX (f df|dx X Y)) + (* DY (f df|dy X Y))))))) + + + + + +(defn D[f] + (fn[x] (f (+ x (differential-seq 1 [0] 1 [1] 1 [2]))))) + +(defn d[partials f] + (fn [x] + (get + (.terms ((D f)x)) + (set partials) + 0 + ))) + +(defmethod exp DifferentialSeq [x] + ((linear-approximator exp exp) x)) + +(defmethod sin DifferentialSeq + [x] + ((linear-approximator sin cos) x)) + +(defmethod cos DifferentialSeq + [x] + ((linear-approximator cos #(- (sin %))) x)) + +(defmethod log DifferentialSeq + [x] + ((linear-approximator log (fn [x] (/ x)) ) x)) + +(defmethod / [DifferentialSeq DifferentialSeq] + [x y] + ((linear-approximator / + (fn [x y] (/ 1 y)) + (fn [x y] (- (/ x (* y y))))) + x y)) diff -r 8d8278e09888 -r b4de894a1e2e sicm/utils.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/utils.html Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,702 @@ + + + + +Building a Classical Mechanics Library in Clojure + + + + + + + + + + + + + +

+ + + + + +

Building a Classical Mechanics Library in Clojure

+ + + + + + + +

1 Useful data types +
- 1.1 Complex numbers
- 1.2 Power series
- 1.3 Tuples and tensors +
  - 1.3.1 Tuples are “sequences with spin”
  - 1.3.2 Matrices
  +
- 1.4 Generic Operations
+
2 Operators and Differentiation +
- 2.1 Operators
- 2.2 Differential Terms and Sequences
+

+ +

1 Useful data types

+ + + +

+ +

1.1 Complex numbers

+ + +

+ +

1.2 Power series

+ + +

+ +

1.3 Tuples and tensors

+ + + +

+ +

1.3.1 Tuples are “sequences with spin”

+ + + + + +

(in-ns 'sicm.utils)
+
+;; Let some objects have spin
+
+(defprotocol Spinning
+  (up? [this])
+  (down? [this]))
+
+(defn spin
+  "Returns the spin of the Spinning s, either :up or :down"
+  [#^Spinning s]
+  (cond (up? s) :up (down? s) :down))
+
+
+;; DEFINITION: A tuple is a sequence with spin
+
+(deftype Tuple
+  [spin coll]
+
+  clojure.lang.Seqable
+  (seq [this] (seq (.coll this)))
+
+  clojure.lang.Counted
+  (count [this] (count (.coll this)))
+
+  Spinning
+  (up? [this] (= ::up (.spin this)))
+  (down? [this] (= ::down (.spin this))))
+
+(defmethod print-method Tuple
+  [o w]
+  (print-simple
+   (if (up? o)
+     (str "u" (.coll o))
+     (str "d" (vec(.coll o))))
+   w))
+
+
+;; CONSTRUCTORS
+
+(defn up
+  "Create a new up-tuple containing the contents of coll."
+  [coll]
+  (Tuple. ::up coll))       
+
+(defn down
+  "Create a new down-tuple containing the contents of coll."
+  [coll]
+  (Tuple. ::down coll))
+

+ + + + + + +

+ +

1.3.2 Matrices

+ + + + +

(in-ns 'sicm.utils)
+(require 'incanter.core) ;; use incanter's fast matrices
+
+
+(defn all-equal? [coll]
+  (if (empty? (rest coll)) true
+      (and (= (first coll) (second coll))
+           (recur (rest coll)))))
+
+
+(defprotocol Matrix
+  (rows [matrix])
+  (cols [matrix])
+  (diagonal [matrix])
+  (trace [matrix])
+  (determinant [matrix])
+  (transpose [matrix])
+  (conjugate [matrix])
+)
+
+(extend-protocol Matrix incanter.Matrix
+  (rows [rs] (map down (apply map vector (apply map vector rs))))
+  (cols [rs] (map up (apply map vector rs)))
+  (diagonal [matrix] (incanter.core/diag matrix) )
+  (determinant [matrix] (incanter.core/det matrix))
+  (trace [matrix] (incanter.core/trace matrix))
+  (transpose [matrix] (incanter.core/trans matrix)))
+
+(defn count-rows [matrix]
+  ((comp count rows) matrix))
+
+(defn count-cols [matrix]
+  ((comp count cols) matrix))
+
+(defn square? [matrix]
+  (= (count-rows matrix) (count-cols matrix)))
+
+(defn identity-matrix
+  "Define a square matrix of size n-by-n with 1s along the diagonal and
+  0s everywhere else."
+  [n]
+  (incanter.core/identity-matrix n))
+
+
+(defn matrix-by-rows
+  "Define a matrix by giving its rows."
+  [& rows]
+  (if
+   (not (all-equal? (map count rows)))
+   (throw (Exception. "All rows in a matrix must have the same number of elements."))
+   (incanter.core/matrix (vec rows))))
+
+(defn matrix-by-cols
+  "Define a matrix by giving its columns"
+  [& cols]
+  (if (not (all-equal? (map count cols)))
+   (throw (Exception. "All columns in a matrix must have the same number of elements."))
+   (incanter.core/matrix (vec (apply map vector cols)))))
+
+(defn identity-matrix
+  "Define a square matrix of size n-by-n with 1s along the diagonal and
+  0s everywhere else."
+  [n]
+  (incanter.core/identity-matrix n))
+
+

+ + + + +

+ +

1.4 Generic Operations

+ + + + +

(in-ns 'sicm.utils)
+(use 'clojure.contrib.generic.arithmetic
+     'clojure.contrib.generic.collection
+     'clojure.contrib.generic.functor
+     'clojure.contrib.generic.math-functions)
+
+(defn numbers?
+  "Returns true if all arguments are numbers, else false."
+  [& xs]
+  (every? number? xs))
+
+(defn tuple-surgery
+  "Applies the function f to the items of tuple and the additional
+  arguments, if any. Returns a Tuple of the same type as tuple."
+  [tuple f & xs]
+  ((if (up? tuple) up down)
+   (apply f (seq tuple) xs)))
+
+
+
+;;; CONTRACTION collapses two compatible tuples into a number.
+
+(defn contractible?
+  "Returns true if the tuples a and b are compatible for contraction,
+  else false. Tuples are compatible if they have the same number of
+  components, they have opposite spins, and their elements are
+  pairwise-compatible."
+  [a b]
+  (and
+   (isa? (type a) Tuple)
+   (isa? (type b) Tuple)
+   (= (count a) (count b))
+   (not= (spin a) (spin b))
+   
+   (not-any? false?
+             (map #(or
+                    (numbers? %1 %2)
+                    (contractible? %1 %2))
+                  a b))))
+
+(defn contract
+  "Contracts two tuples, returning the sum of the
+  products of the corresponding items. Contraction is recursive on
+  nested tuples."
+  [a b]
+  (if (not (contractible? a b))
+    (throw
+     (Exception. "Not compatible for contraction."))
+    (reduce +
+            (map
+             (fn [x y]
+               (if (numbers? x y)
+                 (* x y)
+                 (contract x y)))
+             a b))))
+
+
+
+
+
+(defmethod conj Tuple
+  [tuple & xs]
+  (tuple-surgery tuple #(apply conj % xs)))
+
+(defmethod fmap Tuple
+  [f tuple]
+  (tuple-surgery tuple (partial map f)))
+
+
+
+;; TODO: define Scalar, and add it to the hierarchy above Number and Complex
+
+                                        
+(defmethod * [Tuple Tuple]             ; tuple*tuple
+  [a b]
+  (if (contractible? a b)
+    (contract a b)
+    (map (partial * a) b)))
+
+
+(defmethod * [java.lang.Number Tuple]  ;; scalar *  tuple
+  [a x] (fmap (partial * a) x))
+
+(defmethod * [Tuple java.lang.Number]
+  [x a] (* a x))
+
+(defmethod * [java.lang.Number incanter.Matrix] ;; scalar *  matrix
+  [x M] (incanter.core/mult x M))
+
+(defmethod * [incanter.Matrix java.lang.Number]
+  [M x] (* x M))
+
+(defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
+  [M1 M2]
+  (incanter.core/mmult M1 M2))
+
+(defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
+  [M v]
+  (if (and (apply numbers? v) (up? v)) 
+    (* M (matrix-by-cols v))
+    (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
+    ))
+
+(defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
+  [v M]
+  (if (and (apply numbers? v) (down? v))
+    (* (matrix-by-rows v) M)
+    (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
+    ))
+
+
+(defmethod exp incanter.Matrix
+  [M]
+  (incanter.core/exp M))
+
+

+ + + + +

+ +

2 Operators and Differentiation

+ + + + +

+ +

2.1 Operators

+ + +

+ +

2.2 Differential Terms and Sequences

+ + + + +

(in-ns 'sicm.utils)
+(use 'clojure.contrib.seq
+     'clojure.contrib.generic.arithmetic
+     'clojure.contrib.generic.collection
+     'clojure.contrib.generic.math-functions)
+
+;;∂
+
+;; DEFINITION : Differential Term
+
+;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
+;; coefficient of the quantity is returned by the 'coefficient' method,
+;; while the sequence of differential parameters is returned by the
+;; method 'partials'.
+
+;; Instead of using (potentially ambiguous) letters to denote
+;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].
+
+;; The coefficient can be any arithmetic object; the
+;; partials must be a nonrepeating sorted sequence of nonnegative
+;; integers.
+
+(deftype DifferentialTerm [coefficient partials])
+
+(defn differential-term
+  "Make a differential term given pairs of coefficients and partials."
+  [coefficient partials]
+  (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials)) 
+    (DifferentialTerm. coefficient (set partials))
+    (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))
+
+
+;; DEFINITION : Differential Sequence
+;; A differential sequence is a sequence of differential terms, all with different partials.
+;; Internally, it is a map from the partials of each term to their coefficients.
+
+(deftype DifferentialSeq
+  [terms]
+  ;;clojure.lang.IPersistentMap
+  clojure.lang.Associative
+  (assoc [this key val]
+    (DifferentialSeq.
+     (cons (differential-term val key) terms)))
+  (cons [this x]
+        (DifferentialSeq. (cons x terms)))
+  (containsKey [this key]
+               (not(nil? (find-first #(= (.partials %) key) terms))))
+  (count [this] (count (.terms this)))
+  (empty [this] (DifferentialSeq. []))
+  (entryAt [this key]
+           ((juxt #(.partials %) #(.coefficient %))
+            (find-first #(= (.partials %) key) terms)))
+  (seq [this] (seq (.terms this))))
+
+
+(defmethod fmap DifferentialSeq
+  [f dseq]
+  (DifferentialSeq. (fmap f (.terms dseq))))
+
+
+;; BUILDING DIFFERENTIAL OBJECTS
+
+(defn differential-seq
+    "Define a differential sequence by specifying coefficient/partials
+pairs, which are used to create differential terms to populate the
+sequence."
+  ([coefficient partials]
+     (DifferentialSeq. {(set partials) coefficient}))
+  ([coefficient partials & cps]
+     (if (odd? (count cps))
+       (throw (Exception. "differential-seq requires an even number of terms."))
+       (DifferentialSeq.
+        (reduce
+         #(assoc %1 (set (second %2)) (first %2))
+         {(set partials) coefficient}
+         (partition 2 cps))))))
+  
+
+
+(defn big-part
+  "Returns the part of the differential sequence that is finite,
+  i.e. not infinitely small, as a differential sequence. If given a
+  non-differential object, returns that object."
+  [dseq]
+  (if (not= (type dseq) DifferentialSeq) dseq
+      (let [m (.terms dseq)
+            keys (sort-by count (keys m))
+            smallest-var (last (last keys))]
+        (DifferentialSeq.
+         (reduce
+          #(assoc %1 (first %2) (second %2))
+          {}
+          (remove #((first %) smallest-var) m))))))
+
+
+(defn small-part
+  "Returns the part of the differential sequence that infinitely
+  small, as a differential sequence. If given a non-differential
+  object, returns that object."
+  [dseq]
+  (if (not= (type dseq) DifferentialSeq) 0
+      (let [m (.terms dseq)
+            keys (sort-by count (keys m))
+            smallest-var (last (last keys))]
+        (DifferentialSeq.
+         (reduce
+          #(assoc %1 (first %2) (second %2)) {}
+          (filter #((first %) smallest-var) m))))))
+
+
+
+(defn cartesian-product [set1 set2]
+  (reduce concat
+          (for [x set1]
+            (for [y set2]
+              [x y]))))
+
+(defn differential-multiply
+  "Multiply two differential sequences. The square of any differential
+  variable is zero since differential variables are infinitesimally
+  small."
+  [dseq1 dseq2]
+  (DifferentialSeq.
+   (reduce
+    (fn [m [[vars1 coeff1] [vars2 coeff2]]]
+      (if (not (empty? (clojure.set/intersection vars1 vars2)))
+        m
+        (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
+    {}
+    (cartesian-product (.terms dseq1) (.terms dseq2)))))
+
+
+
+(defmethod * [DifferentialSeq DifferentialSeq]
+  [dseq1 dseq2]
+   (differential-multiply dseq1 dseq2))
+
+(defmethod + [DifferentialSeq DifferentialSeq]
+  [dseq1 dseq2]
+  (DifferentialSeq.
+   (merge-with + (.terms dseq1) (.terms dseq2))))
+
+(defmethod * [java.lang.Number DifferentialSeq]
+  [x dseq]
+  (fmap (partial * x) dseq))
+
+(defmethod * [DifferentialSeq java.lang.Number]
+  [dseq x]
+  (fmap (partial * x) dseq))
+;; DERIVATIVES
+
+
+
+(defn linear-approximator
+  "Returns an operator that linearly approximates the given function."
+  ([f df|dx]
+     (fn [x]
+       (let [big-part (big-part x)
+             small-part (small-part x)]
+         (+ (fmap f big-part)
+            (* (fmap df|dx big-part) small-part)))))
+  ([f df|dx df|dy]
+     (fn [x y]
+       (let [X (big-part x)
+             Y (big-part y)
+             DX (small-part x)
+             DY (small-part y)]
+         (+ (f X Y)
+            (+ (* DX (df|dx X Y))
+               (* DY (df|dy X Y)))))))
+  )
+
+
+
+
+
+(defn D[f]
+  (fn[x]
+    (f (differential-seq x [] 1 [0]))))
+
+(defn d[f partials]
+  (fn [x]
+    (get 
+     (.terms ((D f)x))
+     (set partials)
+     0
+    )))
+
+(defmethod exp DifferentialSeq [x]
+           ((linear-approximator exp exp) x))
+
+(defmethod sin DifferentialSeq
+  [x]
+  ((linear-approximator sin cos) x))
+
+(defmethod cos DifferentialSeq
+  [x]
+  ((linear-approximator cos #(- (sin %))) x))
+     
+
+
+

+ + + + + + + + + +

Date: 2011-08-11 04:10:36 EDT

Author: Dylan Holmes

Org version 7.6 with Emacs version 23

+Validate XHTML 1.0 +

+ + diff -r 8d8278e09888 -r b4de894a1e2e sicm/utils.org --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sicm/utils.org Fri Oct 28 00:03:05 2011 -0700 @@ -0,0 +1,690 @@ +#+TITLE: Building a Classical Mechanics Library in Clojure +#+author: Dylan Holmes +#+EMAIL: rlm@mit.edu +#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js" +#+STYLE: +#+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t +#+SETUPFILE: ../templates/level-0.org +#+INCLUDE: ../templates/level-0.org +#+BABEL: :noweb yes :results silent + +* Useful data types + +** Complex numbers + +** Power series + +** Tuples and tensors + +*** Tuples are \ldquo{}sequences with spin\rdquo{} + +#+srcname: tuples +#+begin_src clojure +(in-ns 'sicm.utils) + +;; Let some objects have spin + +(defprotocol Spinning + (up? [this]) + (down? [this])) + +(defn spin + "Returns the spin of the Spinning s, either :up or :down" + [#^Spinning s] + (cond (up? s) :up (down? s) :down)) + + +;; DEFINITION: A tuple is a sequence with spin + +(deftype Tuple + [spin coll] + + clojure.lang.Seqable + (seq [this] (seq (.coll this))) + + clojure.lang.Counted + (count [this] (count (.coll this))) + + Spinning + (up? [this] (= ::up (.spin this))) + (down? [this] (= ::down (.spin this)))) + +(defmethod print-method Tuple + [o w] + (print-simple + (if (up? o) + (str "u" (.coll o)) + (str "d" (vec(.coll o)))) + w)) + +(def tuple? #(= (type %) Tuple)) + +;; CONSTRUCTORS + +(defn up + "Create a new up-tuple containing the contents of coll." + [coll] + (Tuple. ::up coll)) + +(defn down + "Create a new down-tuple containing the contents of coll." + [coll] + (Tuple. ::down coll)) + +(defn same-spin + "Creates a tuple which has the same spin as tuple and which contains +the contents of coll." + [tuple coll] + (if (up? tuple) + (up coll) + (down coll))) + +(defn opposite-spin + "Create a tuple which has opposite spin to tuple and which contains +the contents of coll." + [tuple coll] + (if (up? tuple) + (down coll) + (up coll))) +#+end_src + + + +*** Matrices +#+srcname:matrices +#+begin_src clojure +(in-ns 'sicm.utils) +(require 'incanter.core) ;; use incanter's fast matrices + + +(defn all-equal? [coll] + (if (empty? (rest coll)) true + (and (= (first coll) (second coll)) + (recur (rest coll))))) + + +(defprotocol Matrix + (rows [matrix]) + (cols [matrix]) + (diagonal [matrix]) + (trace [matrix]) + (determinant [matrix]) + (transpose [matrix]) + (conjugate [matrix]) +) + +(extend-protocol Matrix incanter.Matrix + (rows [rs] (map down (apply map vector (apply map vector rs)))) + (cols [rs] (map up (apply map vector rs))) + (diagonal [matrix] (incanter.core/diag matrix) ) + (determinant [matrix] (incanter.core/det matrix)) + (trace [matrix] (incanter.core/trace matrix)) + (transpose [matrix] (incanter.core/trans matrix))) + +(defn count-rows [matrix] + ((comp count rows) matrix)) + +(defn count-cols [matrix] + ((comp count cols) matrix)) + +(defn square? [matrix] + (= (count-rows matrix) (count-cols matrix))) + +(defn identity-matrix + "Define a square matrix of size n-by-n with 1s along the diagonal and + 0s everywhere else." + [n] + (incanter.core/identity-matrix n)) + + +(defn matrix-by-rows + "Define a matrix by giving its rows." + [& rows] + (if + (not (all-equal? (map count rows))) + (throw (Exception. "All rows in a matrix must have the same number of elements.")) + (incanter.core/matrix (vec rows)))) + +(defn matrix-by-cols + "Define a matrix by giving its columns" + [& cols] + (if (not (all-equal? (map count cols))) + (throw (Exception. "All columns in a matrix must have the same number of elements.")) + (incanter.core/matrix (vec (apply map vector cols))))) + +(defn identity-matrix + "Define a square matrix of size n-by-n with 1s along the diagonal and + 0s everywhere else." + [n] + (incanter.core/identity-matrix n)) + +#+end_src + +** Generic Operations +#+srcname:arith-tuple +#+begin_src clojure +(in-ns 'sicm.utils) +(use 'clojure.contrib.generic.arithmetic + 'clojure.contrib.generic.collection + 'clojure.contrib.generic.functor + 'clojure.contrib.generic.math-functions) + +(defn numbers? + "Returns true if all arguments are numbers, else false." + [& xs] + (every? number? xs)) + +(defn tuple-surgery + "Applies the function f to the items of tuple and the additional + arguments, if any. Returns a Tuple of the same type as tuple." + [tuple f & xs] + ((if (up? tuple) up down) + (apply f (seq tuple) xs))) + + + +;;; CONTRACTION collapses two compatible tuples into a number. + +(defn contractible? + "Returns true if the tuples a and b are compatible for contraction, + else false. Tuples are compatible if they have the same number of + components, they have opposite spins, and their elements are + pairwise-compatible." + [a b] + (and + (isa? (type a) Tuple) + (isa? (type b) Tuple) + (= (count a) (count b)) + (not= (spin a) (spin b)) + + (not-any? false? + (map #(or + (numbers? %1 %2) + (contractible? %1 %2)) + a b)))) + +(defn contract + "Contracts two tuples, returning the sum of the + products of the corresponding items. Contraction is recursive on + nested tuples." + [a b] + (if (not (contractible? a b)) + (throw + (Exception. "Not compatible for contraction.")) + (reduce + + (map + (fn [x y] + (if (numbers? x y) + (* x y) + (contract x y))) + a b)))) + + + + + +(defmethod conj Tuple + [tuple & xs] + (tuple-surgery tuple #(apply conj % xs))) + +(defmethod fmap Tuple + [f tuple] + (tuple-surgery tuple (partial map f))) + + + +;; TODO: define Scalar, and add it to the hierarchy above Number and Complex + + +(defmethod * [Tuple Tuple] ; tuple*tuple + [a b] + (if (contractible? a b) + (contract a b) + (map (partial * a) b))) + + +(defmethod * [java.lang.Number Tuple] ;; scalar * tuple + [a x] (fmap (partial * a) x)) + +(defmethod * [Tuple java.lang.Number] + [x a] (* a x)) + +(defmethod * [java.lang.Number incanter.Matrix] ;; scalar * matrix + [x M] (incanter.core/mult x M)) + +(defmethod * [incanter.Matrix java.lang.Number] + [M x] (* x M)) + +(defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix + [M1 M2] + (incanter.core/mmult M1 M2)) + +(defmethod * [incanter.Matrix Tuple] ;; matrix * tuple + [M v] + (if (and (apply numbers? v) (up? v)) + (* M (matrix-by-cols v)) + (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*")) + )) + +(defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix + [v M] + (if (and (apply numbers? v) (down? v)) + (* (matrix-by-rows v) M) + (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*")) + )) + + +(defmethod exp incanter.Matrix + [M] + (incanter.core/exp M)) + +#+end_src + +* Operators and Differentiation +** Operators +#+scrname: operators +#+begin_src clojure +(in-ns 'sicm.utils) +(use 'clojure.contrib.seq + 'clojure.contrib.generic.arithmetic + 'clojure.contrib.generic.collection + 'clojure.contrib.generic.math-functions) + +(defmethod + [clojure.lang.IFn clojure.lang.IFn] + [f g] + (fn [& args] + (+ (apply f args) (apply g args)))) + +(defmethod * [clojure.lang.IFn clojure.lang.IFn] + [f g] + (fn [& args] + (* (apply f args) (apply g args)))) + +(defmethod / [clojure.lang.IFn java.lang.Number] + [f x] + (fn [& args] + (/ (apply f args) x))) + + +(defmethod - [clojure.lang.IFn] + [f] + (fn [& args] + (- (apply f args)))) + +(defmethod - [clojure.lang.IFn clojure.lang.IFn] + [f g] + (fn [& args] + (- (apply f args) (apply g args)))) + +(defmethod pow [clojure.lang.IFn java.lang.Number] + [f x] + (fn [& args] + (pow (apply f args) x))) + + +(defmethod + [java.lang.Number clojure.lang.IFn] + [x f] + (fn [& args] + (+ x (apply f args)))) + +(defmethod * [java.lang.Number clojure.lang.IFn] + [x f] + (fn [& args] + (* x (apply f args)))) + +(defmethod * [clojure.lang.IFn java.lang.Number] + [f x] + (* x f)) +(defmethod + [clojure.lang.IFn java.lang.Number] + [f x] + (+ x f)) + +#+end_src + +** Differential Terms and Sequences +#+srcname: differential +#+begin_src clojure +(in-ns 'sicm.utils) +(use 'clojure.contrib.seq + 'clojure.contrib.generic.arithmetic + 'clojure.contrib.generic.collection + 'clojure.contrib.generic.math-functions) + +;;∂ + +;; DEFINITION : Differential Term + +;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The +;; coefficient of the quantity is returned by the 'coefficient' method, +;; while the sequence of differential parameters is returned by the +;; method 'partials'. + +;; Instead of using (potentially ambiguous) letters to denote +;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2]. + +;; The coefficient can be any arithmetic object; the +;; partials must be a nonrepeating sorted sequence of nonnegative +;; integers. + +;; (deftype DifferentialTerm [coefficient partials]) + +;; (defn differential-term +;; "Make a differential term from a coefficient and list of partials." +;; [coefficient partials] +;; (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials)) +;; (DifferentialTerm. coefficient (set partials)) +;; (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers.")))) + + +;; DEFINITION : Differential Sequence +;; A differential sequence is a sequence of differential terms, all with different partials. +;; Internally, it is a map from the partials of each term to their coefficients. + +(deftype DifferentialSeq + [terms] + ;;clojure.lang.IPersistentMap + clojure.lang.Associative + (assoc [this key val] + (DifferentialSeq. + (cons (differential-term val key) terms))) + (cons [this x] + (DifferentialSeq. (cons x terms))) + (containsKey [this key] + (not(nil? (find-first #(= (.partials %) key) terms)))) + (count [this] (count (.terms this))) + (empty [this] (DifferentialSeq. [])) + (entryAt [this key] + ((juxt #(.partials %) #(.coefficient %)) + (find-first #(= (.partials %) key) terms))) + (seq [this] (seq (.terms this)))) + +(def differential? #(= (type %) DifferentialSeq)) + +(defn zeroth-order? + "Returns true if the differential sequence has at most a constant term." + [dseq] + (and + (differential? dseq) + (every? + #(= #{} %) + (keys (.terms dseq))))) + +(defmethod fmap DifferentialSeq + [f dseq] + (DifferentialSeq. + (fmap f (.terms dseq)))) + + + + +;; BUILDING DIFFERENTIAL OBJECTS + +(defn differential-seq + "Define a differential sequence by specifying an alternating +sequence of coefficients and lists of partials." + ([coefficient partials] + (DifferentialSeq. {(set partials) coefficient})) + ([coefficient partials & cps] + (if (odd? (count cps)) + (throw (Exception. "differential-seq requires an even number of terms.")) + (DifferentialSeq. + (reduce + #(assoc %1 (set (second %2)) (first %2)) + {(set partials) coefficient} + (partition 2 cps)))))) + + + +(defn big-part + "Returns the part of the differential sequence that is finite, + i.e. not infinitely small. If the sequence is zeroth-order, returns + the coefficient of the zeroth-order term instead. " + [dseq] + (if (zeroth-order? dseq) (get (.terms dseq) #{}) + (let [m (.terms dseq) + keys (sort-by count (keys m)) + smallest-var (last (last keys))] + (DifferentialSeq. + (reduce + #(assoc %1 (first %2) (second %2)) + {} + (remove #((first %) smallest-var) m)))))) + + +(defn small-part + "Returns the part of the differential sequence that infinitely + small. If the sequence is zeroth-order, returns zero." + [dseq] + (if (zeroth-order? dseq) 0 + (let [m (.terms dseq) + keys (sort-by count (keys m)) + smallest-var (last (last keys))] + (DifferentialSeq. + (reduce + #(assoc %1 (first %2) (second %2)) {} + (filter #((first %) smallest-var) m)))))) + + + +(defn cartesian-product [set1 set2] + (reduce concat + (for [x set1] + (for [y set2] + [x y])))) + +(defn nth-subset [n] + (if (zero? n) [] + (let [lg2 #(/ (log %) (log 2)) + k (int(java.lang.Math/floor (lg2 n))) + ] + (cons k + (nth-subset (- n (pow 2 k))))))) + +(def all-partials + (lazy-seq (map nth-subset (range)))) + + +(defn differential-multiply + "Multiply two differential sequences. The square of any differential + variable is zero since differential variables are infinitesimally + small." + [dseq1 dseq2] + (DifferentialSeq. + (reduce + (fn [m [[vars1 coeff1] [vars2 coeff2]]] + (if (not (empty? (clojure.set/intersection vars1 vars2))) + m + (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2)))) + {} + (cartesian-product (.terms dseq1) (.terms dseq2))))) + + + +(defmethod * [DifferentialSeq DifferentialSeq] + [dseq1 dseq2] + (differential-multiply dseq1 dseq2)) + +(defmethod + [DifferentialSeq DifferentialSeq] + [dseq1 dseq2] + (DifferentialSeq. + (merge-with + (.terms dseq1) (.terms dseq2)))) + +(defmethod * [java.lang.Number DifferentialSeq] + [x dseq] + (fmap (partial * x) dseq)) + +(defmethod * [DifferentialSeq java.lang.Number] + [dseq x] + (fmap (partial * x) dseq)) + +(defmethod + [java.lang.Number DifferentialSeq] + [x dseq] + (+ (differential-seq x []) dseq)) +(defmethod + [DifferentialSeq java.lang.Number] + [dseq x] + (+ dseq (differential-seq x []))) + +(defmethod - DifferentialSeq + [x] + (fmap - x)) + + +;; DERIVATIVES + + + +(defn linear-approximator + "Returns an operator that linearly approximates the given function." + ([f df|dx] + (fn [x] + (let [big-part (big-part x) + small-part (small-part x)] + ;; f(x+dx) ~= f(x) + f'(x)dx + (+ (f big-part) + (* (df|dx big-part) small-part) + )))) + + ([f df|dx df|dy] + (fn [x y] + (let [X (big-part x) + Y (big-part y) + DX (small-part x) + DY (small-part y)] + (+ (f X Y) + (* DX (f df|dx X Y)) + (* DY (f df|dy X Y))))))) + + + + + +(defn D[f] + (fn[x] (f (+ x (differential-seq 1 [0] 1 [1] 1 [2]))))) + +(defn d[partials f] + (fn [x] + (get + (.terms ((D f)x)) + (set partials) + 0 + ))) + +(defmethod exp DifferentialSeq [x] + ((linear-approximator exp exp) x)) + +(defmethod sin DifferentialSeq + [x] + ((linear-approximator sin cos) x)) + +(defmethod cos DifferentialSeq + [x] + ((linear-approximator cos #(- (sin %))) x)) + +(defmethod log DifferentialSeq + [x] + ((linear-approximator log (fn [x] (/ x)) ) x)) + +(defmethod / [DifferentialSeq DifferentialSeq] + [x y] + ((linear-approximator / + (fn [x y] (/ 1 y)) + (fn [x y] (- (/ x (* y y))))) + x y)) + +#+end_src + + + +* Derivatives Revisited +#+begin_src clojure +(in-ns 'sicm.utils) +(use 'clojure.contrib.seq + 'clojure.contrib.generic.arithmetic + 'clojure.contrib.generic.collection + 'clojure.contrib.generic.math-functions) + +(defn replace-in + "Replaces the nth item in coll with the given item. If n is larger + than the size of coll, adds n to the end of the collection." + [coll n item] + (concat + (take n coll) + [item] + (drop (inc n) coll))) + +(defn euclidean-structure [f partials] + (fn sd[g v] + (cond + (tuple? v) + (opposite-spin + v + (map + (fn [n] + (sd (fn [xn] + (g + (same-spin v + (replace-in v n xn)))) + (nth v n))) + (range (count v))))))) + + + + +#+end_src + + +* Symbolic Quantities + +#+srcname: symbolic +#+begin_src clojure +(in-ns 'sicm.utils) +(use 'clojure.contrib.generic.arithmetic + 'clojure.contrib.generic.math-functions) + +(deftype Symbolic + [type + expression]) + +(defn print-expression [s] + (print (.expression s))) + +(defn symbolic-number + [symbol] + (Symbolic. java.lang.Number symbol)) + +(defn simplify-expression [s] + (let [e (.expression s)] + (cond + (not(list? e)) e + (= (first e) '+) +) + + + +(defmethod + [Symbolic Symbolic] + [x y] + (Symbolic. (.type x) (list '+ (.expression x) (.expression y)))) + +(defmethod + [Symbolic java.lang.Number] + [s x] + (if (zero? x) + s + (Symbolic. (.type s) (list '+ (.expression s) x)))) + +(defmethod sin Symbolic + [s] + (Symbolic. (.type s) (list 'sin (.expression s)))) + +#+end_src + +* COMMENT To-be-tangled Source +#+begin_src clojure :tangle utils.clj +(ns sicm.utils) + +<> +<> +<> + +<> +#+end_src +

Table of Contents

1 Abstract

2 Introduction

2.1 Simple Loops

2.1.1 Simple Recursion

2.1.2 Iteration

3 Imperative Programming

3.1 Compound Statements

3.2 2.2. The G0 TO Statement

3.3 2.3. Simple Assignment

3.4 Compound Expressions '

4 Continuations

4.1 Continuation-Passing Recursion

4.2 Escape Expressions

Footnotes:

+Monads in CS and Math +

Table of Contents

1 Examples of monads

1.1 The nondeterministic monad

+{{{TITLE}}} +

Categorical Programs

Table of Contents

1 The Category of Types

1.1 Building new data types out of existing ones

{{{title}}}

Table of Contents

1 Deductive and inductive posets

1.1 Definition

1.2 A canonical map from deductive posets to inductive posets

2 Assigning plausibilities to inductive posets

2.1 Consistent reasoning as a commutative diagram

2.2 ``Multiplicity'' is reciprocal probability

2.3 Laws for multiplicity

Footnotes:

aurellem ☉

Synthetic Differential Geometry

Table of Contents

1 Revisiting the real line

1.1 The first anti-euclidean axiom

1.2 The first axiom \(\ldots\) in terms of arrows

1.3 Ex

aurellem ☉

Building a Classical Mechanics Library in Clojure

Table of Contents

1 Generic Arithmetic

2 Useful Data Types

2.1 Complex Numbers

2.2 Tuples and Tensors

2.2.1 Contraction

2.2.2 Matrices

2.3 Power Series

3 Basic Utilities

3.1 Sequence manipulation

3.2 Ranges, Airity and Function Composition

4 Numerical Methods

5 Differentiation

6 Symbolic Manipulation

aurellem ☉

An Unambiguous Notation for Derivatives

Table of Contents

1 Calculus of Infinitesimals

1.1 Differential Objects

1.2 Interactions obey the chain rule

Table of Contents

1 Tuples

2 Generic arithmetic

aurellem ☉

Building a Classical Mechanics Library in Clojure

Table of Contents

1 Useful data types

1.1 Complex numbers

1.2 Power series

1.3 Tuples and tensors

1.3.1 Tuples are “sequences with spin”

1.3.2 Matrices

1.4 Generic Operations

2 Operators and Differentiation

2.1 Operators

2.2 Differential Terms and Sequences