-
-
Notifications
You must be signed in to change notification settings - Fork 68
/
indexed.cljc
367 lines (328 loc) · 14.3 KB
/
indexed.cljc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
;;
;; Copyright © 2021 Sam Ritchie.
;; This work is based on the Scmutils system of MIT/GNU Scheme:
;; Copyright © 2002 Massachusetts Institute of Technology
;;
;; This is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 3 of the License, or (at
;; your option) any later version.
;;
;; This software is distributed in the hope that it will be useful, but
;; WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
;; General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this code; if not, see <http://www.gnu.org/licenses/>.
;;
(ns sicmutils.calculus.indexed
"This namespace implements minimal support for indexed objects and typed
functions."
(:refer-clojure :exclude [+ - * /])
(:require [sicmutils.calculus.basis :as b]
[sicmutils.calculus.form-field :as ff]
[sicmutils.calculus.manifold :as m]
[sicmutils.calculus.vector-field :as vf]
[sicmutils.function :as f]
[sicmutils.generic :as g :refer [+ - * /]]
[sicmutils.operator :as o]
[sicmutils.structure :as s]
[sicmutils.util :as u]
[sicmutils.util.aggregate :as ua]
[sicmutils.util.permute :as permute]
[sicmutils.value :as v]))
;; ## Minimal support for Indexed Objects
;;
;; Indices here are the components of tensors relative to a basis.
(defn- meta-k
"Takes a function or operator `f` and a metadata (or context) key `k` and
attempts to fetch it from the metadata (or context). Returns `default` if `k`
has no entry."
([f k]
(meta-k f k nil))
([f k default]
(if (o/operator? f)
(k (o/context f) default)
(k (meta f) default))))
(defn- with-kv
"Returns a copy of `f` with the `k`, `v` pair added to its metadata (if a
function) or context (if an operator)."
[f k v]
(if (o/operator? f)
(o/with-context f (assoc (o/context f) k v))
(vary-meta f assoc k v)))
(defn argument-types
"Given an operator or function `f`, returns its registered vector of argument
types, or `[]` if none exist.
argument types are, for example,
```clojure
[::ff/oneform-field ::vf/vector-field ::vf/vector-field]
```
for a `Christoffel-2`, which takes one oneform field and two vector fields."
[f]
(meta-k f :arguments []))
(defn ^:no-doc has-argument-types?
"Returns true if `f` has any argument types registered, false otherwise."
[f]
(boolean
(seq (argument-types f))))
(defn with-argument-types
"Given some operator or function `f`, returns a copy of `f` with the supplied
argument types `types` registered in its metadata (if a function) or
context (if an operator).
Retrieve these types with [[argument-types]]."
[f types]
(with-kv f :arguments (into [] types)))
(defn index-types
"Given an operator or function `f`, returns its registered vector of index
types, or `[]` if none exist.
index types are, for example,
```clojure
['up 'down 'down]
```
for a `Christoffel-2`, which takes one oneform field and two vector fields."
[f]
(meta-k f :index-types []))
(defn ^:no-doc has-index-types?
"Returns true if `f` has any index types registered, false otherwise."
[f]
(boolean
(seq (index-types f))))
(defn with-index-types
"Given some operator or function `f`, returns a copy of `f` with the supplied
index types `types` registered in its metadata (if a function) or
context (if an operator).
Retrieve these types with [[index-types]]."
[f types]
(with-kv f :index-types (into [] types)))
;; An argument-typed function of type (n . m) takes n oneform fields and m
;; vector-fields, in that order, and produces a function on a manifold point.
;;
;; An indexed function of type (n . m) takes n+m indices and gives a function on
;; a manifold point.
;;
;; For each argument-typed function and basis, there is an indexed function that
;; gives the function resulting from applying the argument-typed function to the
;; basis elements specified by the corresponding indices.
(defn- valid-arg-types?
"Returns true if `ts` is a well-formed-enough sequence of argument types to use
for generating an indexed function via [[typed->indexed]], false otherwise.
Validates that:
- The sequence of types `ts` is not empty
- every entry in `ts` derives from `::vf/vector-field` or `::ff/oneform-field`
- form fields come before vector fields."
[ts]
(letfn [(one-ff? [t]
(isa? t ::ff/oneform-field))
(vf? [t]
(isa? t ::vf/vector-field))]
(and (seq ts)
(every? (some-fn one-ff? vf?) ts)
(every? vf? (drop-while one-ff? ts)))))
(defn typed->indexed [f basis]
(let [arg-types (argument-types f)]
(assert (valid-arg-types? arg-types)
(str "Invalid arg types: " arg-types))
(let [vector-basis (b/basis->vector-basis basis)
oneform-basis (b/basis->oneform-basis basis)
idx-types (map (fn [t]
(if (isa? t ::vf/vector-field)
'down
'up))
arg-types)]
(-> (fn indexed [indices]
(assert (= (count indices)
(count arg-types))
(str "Indices count doesn't match expected argument types."
" Indices:" indices
", arg-types: " arg-types))
(let [args (mapv (fn [t idx]
(if (isa? t ::vf/vector-field)
(get vector-basis idx)
(get oneform-basis idx)))
arg-types
indices)]
(apply f args)))
(with-index-types idx-types)))))
(defn- valid-index-types?
"Returns true if `ts` is a well-formed-enough sequence of index types to use for
generating a typed function via [[indexed->typed]], false otherwise.
Validates that:
- The sequence of types `ts` is not empty
- every entry in `ts` is either the symbol `'up` or `'down`
- all `'up` entries (corresponding to oneform fields) come before `'down`
entries (corresponding to vector fields)"
[ts]
(boolean
(and (seq ts)
(every? #{'up 'down} ts)
(every? #{'down} (drop-while #{'up} ts)))))
(defn- validate-typed-args!
"Returns true if:
- every argument in `args` has a corresponding index type in `index-types`
- every `'up` in `index-types` is aligned with a [[form-field/oneform-field?]]
argument in `args`
- every `'down` in `index-types` is aligned with a [[vector-field/vector-field?]]
argument in `args`
false otherwise.
`index-types` is assumed to have passed a [[valid-index-types?]] check."
[index-types args]
(assert (= (count index-types)
(count args))
(str "The number "
(count index-types) " of index-types doesn't match the number "
(count args) " of arguments."))
(assert (every? true?
(map (fn [index-type arg]
(or (and (= index-type 'up)
(ff/oneform-field? arg))
(and (= index-type 'down)
(vf/vector-field? arg))))
index-types
args))
(str "Args do not match index-types 'up must be paired with oneform-fields and 'down with vector fields."
" Args:" (pr-str args)
", indices: " (pr-str index-types))))
(defn indexed->typed [indexed basis]
(let [index-types (index-types indexed)]
(assert (valid-index-types? index-types)
(str "Invalid index types: " index-types))
(let [vector-basis (b/basis->vector-basis basis)
oneform-basis (b/basis->oneform-basis basis)
n (b/basis->dimension basis)
arg-types (mapv {'up ::ff/oneform-field
'down ::vf/vector-field}
index-types)]
(-> (fn typed [& args]
(validate-typed-args! index-types args)
(let [n-seq (reverse (range n))
combos (permute/cartesian-product
(for [x args]
(map (fn [i arg]
[i (if (vf/vector-field? arg)
((get oneform-basis i) arg)
(arg (get vector-basis i)))])
n-seq
(repeat x))))]
(ua/generic-sum
(for [combo combos
:let [indices (map first combo)
product-args (map peek combo)]]
(apply *
(indexed indices)
(reverse product-args))))))
(with-argument-types arg-types)))))
(defn outer-product [T1 T2]
(let [i1 (index-types T1)
i2 (index-types T2)]
(assert (seq i1) "No index types registered for T1!")
(assert (seq i2) "No index types registered for T2!")
(let [{nu1 'up nd1 'down} (frequencies i1)
{nu2 'up nd2 'down} (frequencies i2)
nup (+ nu1 nu2)
ndp (+ nd1 nd2)
np (+ nup ndp)
n1 (+ nup nd1)]
(letfn [(product [args]
(assert (= (count args) np)
(str "Wrong number of args to outer-product: "
(count args)
", expected: " np))
(let [argv (into [] args)]
(* (T1 (into (subvec argv 0 nu1)
(subvec argv nup n1)))
(T2 (into (subvec argv nu1 nup)
(subvec argv n1 np))))))]
(with-index-types product
(concat (repeat nup 'up)
(repeat ndp 'down)))))))
(letfn [(insertv [coll i v]
(let [l (subvec coll 0 i)
r (subvec coll i)]
(apply conj l v r)))]
(defn contract [T u d n]
(let [i-types (index-types T)]
(assert (seq i-types) "No index types registered for T!")
(let [{nu 'up nd 'down} (frequencies i-types)]
(assert (and (<= 0 u) (< u nu)
(<= 0 d) (< d nd))
(str "Contraction indices " u ", " d " not in the correct range. "
"Each must be >= 0 and < the respective number of "
"'up and 'down instances in the index types registered with T. "
"These were " nu " and " nd "."))
(let [nuc (dec nu)
ndc (dec nd)]
(-> (fn contraction [args]
(let [argv (into [] args)]
(ua/generic-sum
(fn [i]
(T (concat
(insertv (subvec argv 0 nuc) u i)
(insertv (subvec argv nuc) d i))))
0 n)))
(with-index-types
(concat (repeat nuc 'up)
(repeat ndc 'down)))))))))
(defn typed->structure [T basis]
(let [vector-basis (b/basis->vector-basis basis)
oneform-basis (b/basis->oneform-basis basis)]
(letfn [(lp [arg-types argv]
(if (empty? arg-types)
(apply T argv)
(let [[t & ts] arg-types]
(s/mapr (fn [e]
(lp ts (conj argv e)))
(cond (isa? t ::vf/vector-field)
vector-basis
(isa? t ::ff/oneform-field)
oneform-basis
:else (u/illegal
(str "Invalid argument type: " (pr-str t)
". Every arg must be a vector field or oneform field.")))))))]
(lp (argument-types T) []))))
(defn structure->typed [coeff-functions basis]
(let [vector-basis (b/basis->vector-basis basis)
oneform-basis (b/basis->oneform-basis basis)
arg-types (loop [cf coeff-functions
acc []]
(if-not (s/structure? cf)
acc
(let [shape (s/opposite-orientation
(s/orientation cf))
t (case shape
::s/up ::vf/vector-field
::s/down ::ff/oneform-field)]
(recur (nth cf 0)
(conj acc t)))))]
(-> (fn indexed-fn [& args]
(assert (= (count args)
(count arg-types))
(str "The number of args " (count args)
" did not match the expected arity " (count arg-types) ". "
"Please supply args corresponding to the expected types " arg-types "."))
(assert (every? true? (map (fn [arg arg-type]
(isa? (v/kind arg) arg-type))
args arg-types))
(str "Invalid arguments: " args ". "
"Please supply args corresponding to the expected types " arg-types "."))
(letfn [(lp [args arg-types]
(if (empty? args)
m/one-manifold-function
(let [arg (first args)
arg-type (first arg-types)]
(cond (isa? arg-type ::vf/vector-field)
(s/mapr (fn [etilde]
(* (etilde arg)
(lp (rest args)
(rest arg-types))))
oneform-basis)
(isa? arg-type ::ff/oneform-field)
(s/mapr (fn [e]
(* (arg e)
(lp (rest args)
(rest arg-types))))
vector-basis)))))]
(* (lp args arg-types)
coeff-functions)))
(with-argument-types arg-types))))