Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /branches/vis15/src/compiler/ein/mk-operators.sml
ViewVC logotype

Annotation of /branches/vis15/src/compiler/ein/mk-operators.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3494 - (view) (download)

1 : jhr 3477 (* mk-operators.sml
2 : jhr 3476 *
3 : jhr 3477 * Functions to create the various Ein operators.
4 :     *
5 : jhr 3476 * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
6 :     *
7 :     * COPYRIGHT (c) 2015 The University of Chicago
8 :     * All rights reserved.
9 :     *)
10 :    
11 :     structure MkOperators : sig
12 :    
13 : jhr 3477 type dim = int
14 :     type shape = dim list
15 :     type ids = Ein.index_id list
16 : jhr 3476
17 : jhr 3477 val addRR : Ein.ein
18 :     val addTT : shape -> Ein.ein
19 :     val addTF : dim * shape -> Ein.ein
20 :     val addFF : dim * shape -> Ein.ein
21 :    
22 :     val subRR : Ein.ein
23 :     val subTT : shape -> Ein.ein
24 :     val subTF : dim * shape -> Ein.ein
25 : jhr 3478 val subFT : dim * shape -> Ein.ein
26 : jhr 3477 val subFF : dim * shape -> Ein.ein
27 :    
28 :     val mulRT : shape -> Ein.ein
29 :     val mulRR : Ein.ein
30 :     val mulRF : dim * shape -> Ein.ein
31 :     val mulST : dim * shape -> Ein.ein
32 :     val mulSS : dim -> Ein.ein
33 :     val mulSF : dim * shape -> Ein.ein
34 :    
35 :     val divTR : shape -> Ein.ein
36 :     val divRR : Ein.ein
37 :     val divFR : dim * shape -> Ein.ein
38 :     val divSS : dim -> Ein.ein
39 :     val divFS : dim * shape -> Ein.ein
40 :    
41 :     val negTT : shape -> Ein.ein
42 :     val negFF : dim * shape -> Ein.ein
43 :    
44 :     val cross2TT : Ein.ein
45 :     val cross3TT : Ein.ein
46 :     val cross2FF : Ein.ein
47 :     val cross3FF : Ein.ein
48 :    
49 :     (* QUESTION: why not just "outerTT : shape * shape -> Ein.ein" ? *)
50 :     val outerTT : dim * dim -> Ein.ein
51 :     val outerTM : dim * shape -> Ein.ein
52 :     val outerMT : shape * dim -> Ein.ein
53 :     (* QUESTION: should these be index_kind? *)
54 :     val outerFF : dim * Ein.index_id * Ein.index_id -> Ein.ein
55 :    
56 :     val innerTT : shape * ids -> Ein.ein
57 :     val innerFF : shape * dim * ids -> Ein.ein
58 :     val innerFT : shape * dim * ids -> Ein.ein
59 :     val innerTF : shape * dim * ids -> Ein.ein
60 :    
61 :     val colonTT : shape * ids -> Ein.ein
62 :     val colonFF : dim * shape * ids -> Ein.ein
63 :     val colonFT : dim * shape * ids -> Ein.ein
64 :     val colonTF : dim * shape * ids -> Ein.ein
65 :    
66 : cchiw 3494 val normT : shape -> Ein.ein
67 :     val normF : dim * shape -> Ein.ein
68 : jhr 3477
69 :     val normalizeTT : shape -> Ein.ein
70 :     val normalizeFF : dim * shape -> Ein.ein
71 :    
72 :     val traceT : dim -> Ein.ein
73 :     val traceF : dim * shape -> Ein.ein
74 :    
75 :     val transposeT : shape -> Ein.ein
76 :     (* QUESTION: should these be index_kind? *)
77 :     val transposeF : dim * Ein.index_id * Ein.index_id -> Ein.ein
78 :    
79 : jhr 3478 val det2T : Ein.ein
80 :     val det3T : Ein.ein
81 :     val det2F : Ein.ein
82 :     val det3F_full : Ein.ein
83 :     val det3FA : Ein.ein
84 :     val det3F : Ein.ein
85 : jhr 3477
86 :     val expF : dim -> Ein.ein
87 :     val expT : Ein.ein
88 :    
89 : jhr 3478 val powF : dim * int -> Ein.ein
90 :    
91 : jhr 3477 val sqrtF : dim -> Ein.ein
92 :     val cosF : dim -> Ein.ein
93 :     val acosF : dim -> Ein.ein
94 :     val sinF : dim -> Ein.ein
95 :     val asinF : dim -> Ein.ein
96 :     val tanF : dim -> Ein.ein
97 :     val atanF : dim -> Ein.ein
98 : cchiw 3494 val sqrtR : Ein.ein
99 :     val cosR : Ein.ein
100 :     val acosR : Ein.ein
101 :     val sinR : Ein.ein
102 :     val asinR : Ein.ein
103 :     val tanR : Ein.ein
104 :     val atanR : Ein.ein
105 : jhr 3477
106 :     val modulate : dim -> Ein.ein
107 :    
108 :     val identity : dim -> Ein.ein
109 : jhr 3478 val zeros : shape -> Ein.ein
110 : jhr 3477
111 :     val slice : shape * bool list * int list * shape -> Ein.ein
112 :    
113 :     val conv : dim * shape -> Ein.ein
114 :     val probe : shape * dim -> Ein.ein
115 :    
116 :     val curl2d : Ein.ein
117 :     val curl3d : Ein.ein
118 :     val grad : shape -> Ein.ein
119 :     val dotimes : dim * shape -> Ein.ein (* ?? *)
120 :     val divergence : dim * shape -> Ein.ein
121 :    
122 : jhr 3476 end = struct
123 :    
124 :     structure E = Ein
125 :    
126 : jhr 3477 type dim = int
127 :     type shape = dim list
128 :     type ids = int list
129 :    
130 : jhr 3476 fun specialize (alpha, inc) = List.mapi (fn (i, _) => E.V(i + inc)) alpha
131 :    
132 :     fun sumIds (n, inc, alpha) = let
133 :     val vs = List.tabulate(n, fn v => E.V(v+inc))
134 :     in
135 :     ListPair.map (fn(v, i) => (v, 0, i-1)) (vs, alpha)
136 :     end
137 :    
138 :     fun sumIds2 (n, i) = List.tabulate(n, fn v => (E.V v, 0, i))
139 :    
140 :     val subst_flag = 0 (*here*) (* ????? *)
141 :    
142 :     (******************************* Addition *****************************************)
143 :    
144 :     (* Adding tensors : < X{\alpha} + Y_{\alpha}>_{\alpha} *)
145 :     fun addTT alpha = let
146 :     val expindex = specialize(alpha, 0)
147 :     in
148 :     E.EIN{
149 :     params = [E.TEN(subst_flag, alpha), E.TEN(subst_flag, alpha)],
150 :     index = alpha,
151 :     body = E.Opn(E.Add, [E.Tensor(0, expindex), E.Tensor(1, expindex)])
152 :     }
153 :     end
154 :    
155 :     val addRR = addTT []
156 :    
157 : jhr 3477 (* Tensor and Fields *)
158 :     fun addTF (dim, shape) =let
159 :     val expindex = specialize(shape, 0)
160 : jhr 3476 in
161 :     E.EIN{
162 : jhr 3477 params = [E.TEN(subst_flag, shape), E.FLD dim],
163 : jhr 3476 index = shape,
164 : jhr 3477 body = E.Opn(E.Add, [E.Lift(E.Tensor(0, expindex)), E.Field(1, expindex)])
165 : jhr 3476 }
166 :     end
167 :    
168 : jhr 3477 (* Adding Fields : < F{\alpha} + G_{\alpha}>_{\alpha} *)
169 :     fun addFF (dim, shape) =let
170 : jhr 3476 val expindex = specialize(shape, 0)
171 :     in
172 :     E.EIN{
173 : jhr 3477 params = [E.FLD dim, E.FLD dim],
174 : jhr 3476 index = shape,
175 : jhr 3477 body = E.Opn(E.Add, [E.Field(0, expindex), E.Field(1, expindex)])
176 : jhr 3476 }
177 :     end
178 :    
179 :     (********************************* Subtraction **************************************)
180 :    
181 :     fun subTT alpha = let
182 :     val expindex = specialize(alpha, 0)
183 :     in
184 :     E.EIN{
185 : jhr 3477 params = [E.TEN(subst_flag, alpha), E.TEN(subst_flag, alpha)],
186 : jhr 3476 index = alpha,
187 : jhr 3477 body = E.Op2(E.Sub, E.Tensor(0, expindex), E.Tensor(1, expindex))
188 : jhr 3476 }
189 :     end
190 :    
191 :     val subRR = subTT []
192 :    
193 :     fun subTF (dim, shape) = let
194 : jhr 3477 val expindex = specialize(shape, 0)
195 : jhr 3476 in
196 :     E.EIN{
197 : jhr 3477 params = [E.TEN(subst_flag, shape), E.FLD dim],
198 : jhr 3476 index = shape,
199 : jhr 3477 body = E.Opn(E.Add,
200 :     [E.Lift(E.Tensor(0, expindex)), E.Op1(E.Neg, E.Field(1, expindex))])
201 : jhr 3476 }
202 :     end
203 :    
204 :     fun subFT (dim, shape) = let
205 : jhr 3477 val expindex = specialize(shape, 0)
206 : jhr 3476 in
207 :     E.EIN{
208 : jhr 3477 params = [E.TEN(subst_flag, shape), E.FLD dim],
209 : jhr 3476 index = shape,
210 :     body = E.Op2(E.Sub, E.Field(1, expindex), E.Lift(E.Tensor(0, expindex)))
211 :     }
212 :     end
213 :    
214 :     fun subFF (dim, shape) = let
215 : jhr 3477 val expindex = specialize(shape, 0)
216 : jhr 3476 in
217 :     E.EIN{
218 : jhr 3477 params = [E.FLD dim, E.FLD dim],
219 : jhr 3476 index = shape,
220 : jhr 3477 body = E.Op2(E.Sub, E.Field(0, expindex), E.Field(1, expindex))
221 : jhr 3476 }
222 :     end
223 :    
224 :     (********************************** Multiplication *************************************)
225 :    
226 :     (* scalar times tensor product: <s * T_{\alpha}>_{\alpha} *)
227 :     fun mulRT alpha = let
228 : jhr 3477 val expindex = specialize(alpha, 0)
229 : jhr 3476 in
230 :     E.EIN{
231 : jhr 3477 params = [E.TEN(subst_flag, []), E.TEN(subst_flag,alpha)],
232 : jhr 3476 index = alpha,
233 : jhr 3477 body = E.Opn(E.Prod, [E.Tensor(0, []), E.Tensor(1, expindex)])
234 : jhr 3476 }
235 :     end
236 :    
237 :     val mulRR = mulRT []
238 :    
239 : jhr 3477 fun mulRF (dim, shape) =let
240 : jhr 3476 val expindex = specialize(shape, 0)
241 :     in
242 :     E.EIN{
243 : jhr 3477 params = [E.TEN(subst_flag, []), E.FLD dim],
244 : jhr 3476 index = shape,
245 : jhr 3477 body = E.Opn(E.Prod, [E.Lift( E.Tensor(0, [])), E.Field(1,expindex)])
246 : jhr 3476 }
247 :     end
248 :    
249 :     fun mulST (dim, shape) =let
250 :     val expindex = specialize(shape, 0)
251 :     in
252 :     E.EIN{
253 : jhr 3477 params = [E.TEN(subst_flag,shape), E.FLD dim],
254 : jhr 3476 index = shape,
255 : jhr 3477 body = E.Opn(E.Prod, [E.Lift( E.Tensor(0,expindex)), E.Field(1, [])])
256 : jhr 3476 }
257 :     end
258 :    
259 :     fun mulSS dim = E.EIN{
260 : jhr 3477 params = [E.FLD dim, E.FLD dim],
261 : jhr 3476 index = [],
262 : jhr 3477 body = E.Opn(E.Prod, [E.Field(0, []), E.Field(1, [])])
263 : jhr 3476 }
264 :    
265 : jhr 3477 fun mulSF(dim, shape) =let
266 :     val expindex= specialize(shape, 0)
267 : jhr 3476 in E.EIN{
268 : jhr 3477 params = [E.FLD dim, E.FLD dim],
269 : jhr 3476 index = shape,
270 : jhr 3477 body = E.Opn(E.Prod, [E.Field(0, []), E.Field(1, expindex)])
271 : jhr 3476 }
272 :     end
273 :    
274 :     (************************************ Division ************************************)
275 :    
276 : jhr 3477 fun divTR alpha = let
277 :     val expindex = specialize(alpha, 0)
278 :     in
279 : jhr 3476 E.EIN{
280 : jhr 3477 params = [E.TEN(subst_flag, alpha), E.TEN(subst_flag, [])],
281 : jhr 3476 index = alpha,
282 : jhr 3477 body = E.Op2(E.Div, E.Tensor(0, expindex), E.Tensor(1, []))
283 :     }
284 :     end
285 : jhr 3476
286 : jhr 3477 val divRR = divTR []
287 : jhr 3476
288 : jhr 3477 fun divFR (dim, shape) = let
289 :     val expindex = specialize(shape, 0)
290 :     in
291 :     E.EIN{
292 :     params = [E.FLD dim, E.TEN(subst_flag, [])],
293 :     index = shape,
294 :     body = E.Op2(E.Div, E.Field(0, expindex), E.Lift(E.Tensor(1, [])))
295 :     }
296 :     end
297 :    
298 : jhr 3476 fun divSS dim = E.EIN{
299 : jhr 3477 params = [E.FLD dim, E.FLD dim],
300 :     index = [],
301 :     body = E.Op2(E.Div, E.Field(0, []), E.Field(1, []))
302 :     }
303 : jhr 3476
304 : jhr 3477 fun divFS(dim, shape) = let
305 :     val expindex = specialize(shape, 0)
306 :     in
307 :     E.EIN{
308 :     params = [E.FLD dim, E.FLD dim],
309 :     index = shape,
310 :     body = E.Opn(E.Prod, [E.Field(0, expindex), E.Op2(E.Div, E.B(E.Const 1), E.Field(1, []))])
311 :     }
312 :     end
313 : jhr 3476
314 : jhr 3477 (************************************* Negation **********************************)
315 :    
316 : jhr 3476 fun negTT alpha = let
317 : jhr 3477 val expindex= specialize(alpha, 0)
318 : jhr 3476 (*changed tensor lift variable here *)
319 : jhr 3477 in
320 :     E.EIN {
321 :     params = [E.TEN(subst_flag, alpha)], index = alpha,
322 :     body = E.Op1(E.Neg, E.Tensor(0, expindex))
323 :     }
324 :     end
325 : jhr 3476
326 : jhr 3477 fun negFF (dim, shape) = let
327 :     val expindex = specialize(shape, 0)
328 :     in
329 :     E.EIN{
330 :     params = [E.FLD dim], index = shape,
331 :     body = E.Op1(E.Neg, E.Field(0, expindex))
332 :     }
333 :     end
334 : jhr 3476
335 : jhr 3477 (****************************** cross product ***********************************)
336 : jhr 3476
337 : jhr 3477 (* 2-d cross product Eps_{ij}U_i V_j *)
338 : jhr 3476 val cross2TT = E.EIN{
339 : jhr 3477 params = [E.TEN(subst_flag, [2]), E.TEN(subst_flag, [2])],
340 :     index = [],
341 :     body = E.Sum([(E. V 0, 0, 1), (E.V 1, 0, 1)],
342 :     E.Opn(E.Prod, [E.G(E.Eps2(0, 1)), E.Tensor(0, [E.V 0]), E.Tensor(1, [E.V 1])]))
343 :     }
344 : jhr 3476
345 : jhr 3477 (* crossProduct is on 3D vectors ..vec3 t8=t0 × t1; *)
346 :     val cross3TT = E.EIN{
347 :     params = [E.TEN(subst_flag, [3]), E.TEN(subst_flag, [3])],
348 :     index = [3],
349 :     body = E.Sum([(E.V 1, 0, 2), (E.V 2, 0, 2)],
350 :     E.Opn(E.Prod, [
351 :     E.G(E.Epsilon(0, 1, 2)),
352 :     E.Tensor(0, [E.V 1]),
353 :     E.Tensor(1, [E.V 2])
354 :     ]))
355 :     }
356 :    
357 :     (* Field Cross Product *)
358 : jhr 3476 val cross2FF = E.EIN{
359 : jhr 3477 params = [E.FLD(2), E.FLD(2)], index = [],
360 :     body = E.Sum([(E.V 0, 0, 1), (E.V 1, 0, 1)],
361 :     E.Opn(E.Prod, [E.G(E.Eps2(0, 1)), E.Field(0, [E.V 0]), E.Field(1, [E.V 1])]))
362 :     }
363 : jhr 3476
364 : jhr 3477 (* Field Cross Product *)
365 : jhr 3476 val cross3FF = E.EIN{
366 :     params = [E.FLD(3), E.FLD(3)],index= [3],
367 : jhr 3477 body = E.Sum([(E.V 1, 0, 2), (E.V 2, 0, 2)],
368 :     E.Opn(E.Prod, [E.G(E.Epsilon(0, 1, 2)),
369 :     E.Field(0, [E.V 1]), E.Field(1, [E.V 2])]))
370 :     }
371 : jhr 3476
372 : jhr 3477 (******************** outer product ********************************)
373 : jhr 3476
374 : jhr 3477 (* QUESTION: can we generalize outer product? *)
375 :     (*Vector Examples : <T_i * T_j>_ij..t0⊗t1*)
376 :     fun outerTT (dimA, dimB) = E.EIN{
377 :     params = [E.TEN(subst_flag, [dimA]), E.TEN(subst_flag, [dimB])],
378 :     index = [dimA, dimB],
379 :     body = E.Opn(E.Prod, [E.Tensor(0, [E.V 0]), E.Tensor(1, [E.V 1])])
380 :     }
381 :    
382 :     fun outerTM (i, [j, k]) = E.EIN{
383 :     params = [E.TEN(subst_flag, [i]), E.TEN(subst_flag, [j, k])],
384 :     index = [i, j, k],
385 :     body = E.Opn(E.Prod, [E.Tensor(0, [E.V 0]), E.Tensor(1, [E.V 1, E.V 2])])
386 :     }
387 :    
388 :     fun outerMT ([i, j], k) = E.EIN{
389 :     params = [E.TEN(subst_flag, [i, j]), E.TEN(subst_flag, [k])],
390 :     index = [i, j, k],
391 :     body = E.Opn(E.Prod, [E.Tensor(0, [E.V 0, E.V 1]), E.Tensor(1, [E.V 2])])
392 :     }
393 :    
394 : jhr 3476 (*Assumes same dimension vector field *)
395 : jhr 3477 fun outerFF (dim, i, j) = E.EIN{
396 :     params = [E.FLD dim, E.FLD dim],
397 :     index = [i, j],
398 :     body = E.Opn(E.Prod, [E.Field(0, [E.V 0]), E.Field(1, [E.V 1])])
399 :     }
400 : jhr 3476
401 : jhr 3477 (*************************** inner product **********************************)
402 : jhr 3476 (* generic inner product: <T_{\alpha i} * T_{i \beta}>_{\alpha \beta} *)
403 : jhr 3477 fun innerTT (shape1, i::beta) = let
404 :     val alpha = List.take(shape1, length shape1 - 1)
405 :     val expindexA = specialize(alpha, 0)
406 :     val expindexB = specialize(beta, length alpha)
407 :     val s' = E.V(length alpha + length beta)
408 :     val s'' = [(s', 0, i-1)]
409 :     in
410 : jhr 3476 E.EIN{
411 : jhr 3477 params = [E.TEN(subst_flag, shape1), E.TEN(subst_flag, i::beta)],
412 : jhr 3476 index = alpha@beta,
413 : jhr 3477 body = E.Sum(s'', E.Opn(E.Prod, [
414 :     E.Tensor(0, expindexA@[s']), (* T_{\alpha i} *)
415 :     E.Tensor(1, [s']@expindexB ) (* T'_{i \beta} *)
416 :     ]))
417 :     }
418 :     end
419 : jhr 3476 | innerTT _ = raise Fail "Wrong shape for inner product"
420 :    
421 : jhr 3477 (* generic inner product: <T_{\alpha i} * T_{i \beta}>_{\alpha \beta} *)
422 :     fun innerFF (shape1, dim, i::beta) = let
423 :     val alpha = List.take(shape1, length shape1 - 1)
424 :     val expindexA = specialize(alpha, 0)
425 :     val expindexB = specialize(beta, length alpha)
426 :     val sid = E.V(length alpha + length beta)
427 :     in
428 :     E.EIN{
429 :     params = [E.FLD dim, E.FLD dim],
430 :     index = alpha@beta,
431 :     body = E.Sum([(sid, 0, i-1)],
432 :     E.Opn(E.Prod, [
433 :     E.Field(0, expindexA@[sid]), (* F_{\alpha i} *)
434 :     E.Field(1, [sid]@expindexB) (* F'_{i \beta} *)
435 :     ]))
436 :     }
437 :     end
438 :     | innerFF _ = raise Fail "Wrong shape for innerProductField"
439 : jhr 3476
440 : jhr 3477 fun innerFT (shape1, dim, i::beta) = let
441 :     val alpha = List.take(shape1, length shape1-1)
442 :     val expindexA = specialize(alpha, 0)
443 :     val expindexB = specialize(beta, length alpha)
444 :     val sid = E.V(length alpha + length beta)
445 :     in
446 :     E.EIN{
447 :     params = [E.FLD dim, E.TEN(subst_flag, i::beta)],
448 :     index = alpha@beta,
449 :     body = E.Sum([(sid, 0, i-1)],
450 :     E.Opn(E.Prod, [
451 :     E.Field(0, expindexA@[sid]), (* F_{\alpha i} *)
452 :     E.Lift(E.Tensor(1, [sid]@expindexB )) (* F'_{i \beta} *)
453 :     ]))
454 :     }
455 :     end
456 :     | innerFT _ = raise Fail "Wrong shape for innerProductFieldTensor"
457 : jhr 3476
458 : jhr 3477 fun innerTF (shape1, dim, i::beta) = let
459 :     val alpha = List.take(shape1, length shape1 - 1)
460 :     val expindexA = specialize(alpha, 0)
461 :     val expindexB = specialize(beta, length alpha)
462 :     val sid = E.V(length(alpha) + length beta)
463 :     in
464 :     E.EIN{
465 :     params = [E.TEN(subst_flag,shape1), E.FLD dim],index = alpha@beta,
466 :     body = E.Sum([(sid, 0, i-1)],
467 :     E.Opn(E.Prod, [
468 :     E.Lift(E.Tensor(0, expindexA@[sid])), (* F_{\alpha i} *)
469 :     E.Field(1, [sid]@expindexB) (* F'_{i \beta} *)
470 :     ]))
471 :     }
472 :     end
473 :     | innerTF _ = raise Fail "Wrong shape for innerProductTensorField"
474 : jhr 3476
475 : jhr 3477 (*************************** colon product **********************************)
476 : jhr 3476
477 : jhr 3477 (* <T_{\alpha i j} * B{i j \beta }>_\alpha \beta *)
478 :     fun colonTT (shape1, i::j::beta) = let
479 :     val lenAlpha = length shape1 - 2
480 :     val alpha = List.take(shape1, lenAlpha)
481 :     val expindexA = specialize(alpha, 0)
482 :     val expindexB = specialize(beta, lenAlpha)
483 :     val sumi = lenAlpha + length beta
484 :     val s' = [E.V sumi, E.V(sumi+1)]
485 :     val sx = [(E.V sumi, 0, i-1), (E.V(sumi+1), 0, j-1)]
486 :     in
487 : jhr 3476 E.EIN{
488 : jhr 3477 params = [E.TEN(subst_flag, shape1), E.TEN(subst_flag, i::j::beta)],
489 : jhr 3476 index = alpha@beta,
490 : jhr 3477 body = E.Sum(sx,
491 :     E.Opn(E.Prod, [E.Tensor(0, expindexA@s'), E.Tensor(1, s'@expindexB)]))
492 :     }
493 :     end
494 : jhr 3476
495 : jhr 3477 (* <F_{\alpha i j} * G_{i j \beta }>_\alpha \beta *)
496 :     fun colonFF (dim, shape1, i::j::beta) = let
497 :     val lenAlpha = length shape1 - 2
498 :     val alpha = List.take(shape1, lenAlpha)
499 :     val expindexA = specialize(alpha, 0)
500 :     val expindexB = specialize(beta, lenAlpha)
501 :     val sumi = lenAlpha + length beta
502 :     val s' = [E.V sumi, E.V(sumi+1)]
503 :     val sx = [(E.V sumi, 0, i-1), (E.V(sumi+1), 0, j-1)]
504 :     in
505 : jhr 3476 E.EIN{
506 : jhr 3477 params = [E.FLD dim, E.FLD dim],
507 : jhr 3476 index = alpha@beta,
508 : jhr 3477 body = E.Sum(sx,
509 :     E.Opn(E.Prod, [E.Field(0, expindexA@s'), E.Field(1, s'@expindexB)]))
510 :     }
511 :     end
512 : jhr 3476
513 :    
514 : jhr 3477 (* <F_{\alpha i j} * T_{i j \beta }>_\alpha \beta *)
515 :     fun colonFT (dim, shape1, i::j::beta) = let
516 :     val lenAlpha = length shape1 - 2
517 :     val alpha = List.take(shape1, lenAlpha)
518 :     val expindexA = specialize(alpha, 0)
519 :     val expindexB = specialize(beta, lenAlpha)
520 :     val sumi = lenAlpha + length beta
521 :     val s' = [E.V sumi, E.V(sumi+1)]
522 :     val sx = [(E.V sumi, 0, i-1), (E.V(sumi+1), 0, j-1)]
523 :     in
524 : jhr 3476 E.EIN{
525 : jhr 3477 params = [E.FLD dim, E.TEN(subst_flag,shape1)],
526 :     index = alpha@beta,
527 :     body = E.Sum(sx,
528 :     E.Opn(E.Prod, [E.Field(0, expindexA@s'), E.Lift(E.Tensor(1, s'@expindexB))]))
529 :     }
530 :     end
531 : jhr 3476
532 : jhr 3477 (* <T_{\alpha i j} * G{i j \beta }>_\alpha \beta *)
533 :     fun colonTF (dim, shape1, i::j::beta) = let
534 :     val lenAlpha = length shape1 - 2
535 :     val alpha = List.take(shape1, lenAlpha)
536 :     val expindexA = specialize(alpha, 0)
537 :     val expindexB = specialize(beta, lenAlpha)
538 :     val sumi = lenAlpha + length beta
539 :     val s' = [E.V sumi, E.V(sumi+1)]
540 :     val sx = [(E.V sumi, 0, i-1), (E.V(sumi+1), 0, j-1)]
541 :     in
542 : jhr 3476 E.EIN{
543 : jhr 3477 params = [E.TEN(subst_flag, i::j::beta), E.FLD dim],
544 :     index = alpha@beta,
545 :     body = E.Sum(sx,
546 :     E.Opn(E.Prod, [E.Lift(E.Tensor(0, expindexA@s')), E.Field(1, s'@expindexB)]))
547 :     }
548 :     end
549 : jhr 3476
550 : jhr 3477 (******************** Norm ********************************)
551 : jhr 3476
552 :     (*n.t.s tensor norms/normalize use high-il ops
553 :     *we can norm any size and use einapp
554 :     * normalize only implemented for vectors and use mid il-op
555 :     *)
556 : jhr 3477 (* get norm, but without the sqrt. Implemented as a summation over a modulate *)
557 : cchiw 3494 fun normTT alpha = let
558 :     val expindex= specialize(alpha,0)
559 :     val sx= sumIds(length(alpha),0,alpha)
560 :     in
561 : jhr 3476 E.EIN{
562 : cchiw 3494 params = [E.TEN(subst_flag,alpha)],
563 :     index = [],
564 :     body = E.Op1(E.Sqrt,E.Sum(sx, E.Opn(E.Prod,[E.Tensor(0, expindex), E.Tensor(0, expindex)])))
565 :     }
566 :     end
567 :    
568 :     fun normFF(dim,alpha) =let
569 :     val expindex= specialize(alpha,0)
570 :     val sx= sumIds(length(alpha),0,alpha)
571 :     in
572 :     E.EIN{
573 :     params = [E.FLD(dim)],
574 :     index = [],
575 :     body = E.Op1(E.Sqrt,E.Sum(sx, E.Opn(E.Prod,[E.Field(0, expindex),E.Field(0,expindex)])))
576 :     }
577 :     end
578 : jhr 3476
579 : cchiw 3494 fun normalizeTT alpha = let
580 : jhr 3477 val expindex = specialize(alpha, 0)
581 :     val len = length alpha
582 :     val expindexDot = specialize(alpha, len)
583 :     val sx = sumIds(len, len, alpha)
584 :     val f = E.Tensor(0, expindex)
585 :     val g = E.Tensor(1, expindexDot)
586 :     in
587 : jhr 3476 E.EIN{
588 : jhr 3477 params = [E.TEN(subst_flag, alpha), E.TEN(subst_flag,alpha)],
589 :     index = alpha,
590 :     body = E.Opn(E.Prod, [
591 :     f, E.Op2(E.Div, E.B(E.Const 1), E.Op1(E.Sqrt, E.Sum(sx, E.Opn(E.Prod, [g, g]))))
592 :     ])
593 :     }
594 :     end
595 : jhr 3476
596 : cchiw 3494 fun normalizeFF (dim, alpha as i::_) = let
597 : jhr 3477 val expindex = specialize(alpha, 0)
598 :     val len = length alpha
599 :     val expindexDot = specialize(alpha, len)
600 :     val sx = sumIds(len, len, alpha)
601 :     val f = E.Field(0, expindex)
602 :     val g = E.Field(1, expindexDot)
603 :     in
604 : jhr 3476 E.EIN{
605 : jhr 3477 params = [E.FLD dim, E.FLD dim],
606 : jhr 3476 index = [dim],
607 : jhr 3477 body = E.Opn(E.Prod, [
608 :     f, E.Op2(E.Div, E.B(E.Const 1), E.Op1(E.Sqrt, E.Sum(sx, E.Opn(E.Prod, [g, g]))))
609 :     ])
610 :     }
611 :     end
612 : jhr 3476
613 : jhr 3477 (************************* trace *************************)
614 :    
615 :     (* Trace: <M_{i, i}> This one Sx represents both i's*)
616 : jhr 3476 fun traceT dim = E.EIN{
617 : jhr 3477 params = [E.TEN(1, [dim, dim])], index = [],
618 :     body = E.Sum([(E.V 0, 0, dim-1)], E.Tensor(0, [E.V 0, E.V 0]))
619 :     }
620 : jhr 3476
621 : jhr 3477 (* Trace: <Sigma_i F_{\alpha i, i}> This one Sx represents both i's *)
622 :     fun traceF (dim, alpha) = let
623 :     val expindex = specialize(alpha, 0)
624 :     val s = E.V(length alpha)
625 :     in
626 : jhr 3476 E.EIN{
627 : jhr 3477 params = [E.FLD dim],
628 : jhr 3476 index = alpha,
629 : jhr 3477 body = E.Sum([(s, 0, dim-1)], E.Field(0, expindex@[s, s]))
630 :     }
631 :     end
632 : jhr 3476
633 : jhr 3477 (************************* tranpose *************************)
634 : jhr 3476
635 : jhr 3477 fun transposeT alpha = E.EIN{
636 :     params = [E.TEN(subst_flag, alpha)],
637 :     index = List.rev alpha,
638 :     body = E.Tensor(0, [E.V 1, E.V 0])
639 :     }
640 : jhr 3476
641 : jhr 3477 (* Transpose Field F_{ji} *)
642 :     fun transposeF (dim, i, j) = E.EIN{
643 :     params = [E.FLD dim],
644 :     index = [i, j],
645 :     body = E.Field(0, [E.V 1, E.V 0])
646 :     }
647 : jhr 3476
648 : jhr 3477 (************************* determinant *************************)
649 : jhr 3476
650 : jhr 3478 val det2T = E.EIN{
651 : jhr 3477 params = [E.TEN(0, [2, 2])],
652 :     index = [],
653 :     body = E.Op2(E.Sub,
654 :     E.Opn(E.Prod, [E.Tensor(0, [E.C 0, E.C 0]), E.Tensor(0, [E.C 1, E.C 1])]),
655 :     E.Opn(E.Prod, [E.Tensor(0, [E.C 0, E.C 1]), E.Tensor(0, [E.C 1, E.C 0])]))
656 :     }
657 : jhr 3476
658 : jhr 3478 val det3T = let
659 : jhr 3477 val a = E.Tensor(0, [E.C 0, E.C 0])
660 :     val b = E.Tensor(0, [E.C 0, E.C 1])
661 :     val c = E.Tensor(0, [E.C 0, E.C 2])
662 :     val d = E.Tensor(0, [E.C 1, E.C 0])
663 :     val e = E.Tensor(0, [E.C 1, E.C 1])
664 :     val f = E.Tensor(0, [E.C 1, E.C 2])
665 :     val g = E.Tensor(0, [E.C 2, E.C 0])
666 :     val h = E.Tensor(0, [E.C 2, E.C 1])
667 :     val i = E.Tensor(0, [E.C 2, E.C 2])
668 :     in
669 :     E.EIN{
670 :     params = [E.TEN(0, [3, 3])],
671 :     index = [],
672 :     body = E.Op2(E.Sub,
673 :     E.Opn(E.Add, [
674 :     E.Opn(E.Prod, [a, e, i]), E.Opn(E.Prod, [b, f, g]), E.Opn(E.Prod, [c, d, h])
675 :     ]),
676 :     E.Opn(E.Add, [
677 :     E.Opn(E.Prod, [c, e, g]), E.Opn(E.Prod, [b, d, i]), E.Opn(E.Prod, [a, f, h])
678 :     ]))
679 :     }
680 :     end
681 :    
682 : jhr 3478 val det2F = E.EIN{
683 : jhr 3477 params = [E.FLD 2],
684 :     index = [],
685 :     body = E.Op2(E.Sub,
686 :     E.Opn(E.Prod, [E.Field(0, [E.C 0, E.C 0]), E.Field(0, [E.C 1, E.C 1])]),
687 :     E.Opn(E.Prod, [E.Field(0, [E.C 0, E.C 1]), E.Field(0, [E.C 1, E.C 0])]))
688 :     }
689 :    
690 : jhr 3478 val det3F_full = let
691 : jhr 3477 val a = E.Field(0, [E.C 0, E.C 0])
692 :     val b = E.Field(0, [E.C 0, E.C 1])
693 :     val c = E.Field(0, [E.C 0, E.C 2])
694 :     val d = E.Field(0, [E.C 1, E.C 0])
695 :     val e = E.Field(0, [E.C 1, E.C 1])
696 :     val f = E.Field(0, [E.C 1, E.C 2])
697 :     val g = E.Field(0, [E.C 2, E.C 0])
698 :     val h = E.Field(0, [E.C 2, E.C 1])
699 :     val i = E.Field(0, [E.C 2, E.C 2])
700 :     in
701 :     E.EIN{
702 :     params = [E.FLD 3],
703 :     index = [],
704 :     body = E.Op2(E.Sub,
705 :     E.Opn(E.Add, [
706 :     E.Opn(E.Prod, [a, e, i]), E.Opn(E.Prod, [b, f, g]), E.Opn(E.Prod, [c, d, h])
707 :     ]),
708 :     E.Opn(E.Add, [
709 :     E.Opn(E.Prod, [c, e, g]), E.Opn(E.Prod, [b, d, i]), E.Opn(E.Prod, [a, f, h])
710 :     ]))
711 :     }
712 :     end
713 : jhr 3476
714 : jhr 3478 val det3FA = E.EIN{
715 : jhr 3477 params = [E.FLD 3],
716 :     index = [],
717 :     body = E.Sum([(E.V 0, 0, 2), (E.V 1, 0, 2), (E.V 2, 0, 2)],
718 :     E.Opn(E.Prod, [
719 :     E.G(E.Epsilon(0, 1, 2)),
720 :     E.Field(0, [E.C 0, E. V 0]),
721 :     E.Field(0, [E.C 1, E. V 1]),
722 :     E.Field(0, [E.C 2, E. V 2])
723 :     ]))
724 :     }
725 : jhr 3476
726 : jhr 3478 val det3F = E.EIN{
727 : jhr 3477 params = [E.FLD 3],
728 :     index = [],
729 :     body = E.Sum([(E.V 0, 0, 2)],
730 :     E.Opn(E.Prod, [
731 :     E.Field(0, [E.C 0, E. V 0]),
732 :     E.Sum([(E.V 1, 0, 2)],
733 :     E.Opn(E.Prod, [
734 :     E.Field(0, [E.C 1, E. V 1]),
735 :     E.Sum([(E.V 2, 0, 2)],
736 :     E.Opn(E.Prod, [E.G(E.Epsilon(0, 1, 2)), E.Field(0, [E.C 2, E. V 2])]))
737 :     ]))
738 :     ]))
739 :     }
740 : jhr 3476
741 : jhr 3477 (************************* Exponential **************************)
742 :     fun expF dim = E.EIN{params = [E.FLD dim], index = [], body = E.Op1(E.Exp, E.Field(0, []))}
743 :     val expT = E.EIN{params = [E.TEN(0, [])], index = [], body = E.Op1(E.Exp, E.Tensor(0, []))}
744 : jhr 3476
745 : jhr 3477 (************************* Lifted single-argument math functions *************************)
746 :     local
747 : cchiw 3494 fun liftFn op1 dim = E.EIN{
748 : jhr 3477 params = [E.FLD dim],
749 :     index = [],
750 : cchiw 3494 body = E.Op1(op1, E.Field(0, []))
751 : jhr 3477 }
752 : cchiw 3494 fun op1Fn op0= E.EIN{
753 :     params = [E.TEN(0,[])],
754 :     index = [],
755 :     body = E.Op1(op0, E.Field(0, []))
756 :     }
757 : jhr 3477 in
758 : jhr 3478 fun powF (dim, n) = E.EIN{
759 :     params = [E.FLD dim],
760 :     index = [], body = E.Op1(E.PowInt n, E.Field(0, []))
761 :     }
762 : jhr 3477 val sqrtF = liftFn E.Sqrt
763 :     val cosF = liftFn E.Cosine
764 :     val acosF = liftFn E.ArcCosine
765 :     val sinF = liftFn E.Sine
766 :     val asinF = liftFn E.ArcSine
767 :     val tanF = liftFn E.Tangent
768 :     val atanF = liftFn E.ArcTangent
769 : cchiw 3494 val sqrtR =op1Fn E.Sqrt
770 :     val acosR = op1Fn E.ArcCosine
771 :     val sinR = op1Fn E.Sine
772 :     val asinR = op1Fn E.ArcSine
773 :     val tanR = op1Fn E.Tangent
774 :     val atanR = op1Fn E.ArcTangent
775 : jhr 3477 end (* local *)
776 : jhr 3476
777 : jhr 3477 (************************* other tensor ops *************************)
778 :     fun modulate dim = E.EIN{
779 :     params = [E.TEN(subst_flag, [dim]), E.TEN(subst_flag, [dim])],
780 :     index = [dim],
781 :     body = E.Opn(E.Prod, [E.Tensor(0, [E.V 0]), E.Tensor(1, [E.V 0])])
782 :     }
783 : jhr 3476
784 : jhr 3477 fun identity dim = E.EIN{
785 :     params = [], index = [dim, dim], body = E.G(E.Delta(E.V(0), E.V(1)))
786 :     }
787 : jhr 3476
788 : jhr 3478 fun zeros shape = E.EIN{
789 :     params = [], index = shape, body = E.B(E.Const 0)
790 :     }
791 :    
792 : jhr 3477 (* QUESTION: why do we need the const list? The indices are implicit in the position of
793 :     * of the mask element! Likewise, the result type can be determined from the argTy and
794 :     * mask.
795 :     *)
796 :     fun slice (argTy, mask, const, rstTy) = let
797 :     fun iter ([], _, cnt) = []
798 :     | iter (true::es, c::cs, cnt) = E.C c :: iter(es, cs, cnt)
799 :     | iter (false::es, cs, cnt) = E.V cnt :: iter(es, cs, cnt+1)
800 :     val ix = iter(mask, const, 0)
801 :     in
802 :     E.EIN{
803 :     params = [E.TEN(subst_flag, argTy)],
804 :     index = rstTy,
805 :     body = E.Tensor(0, ix)
806 :     }
807 :     end
808 : jhr 3476
809 : jhr 3477 (******************** other field ops ********************************)
810 : jhr 3476
811 : jhr 3477 (* FLD here is bounded to image field, and dimension of h *)
812 :     fun conv (dim, shape) =let
813 :     val expindex = specialize(shape, 0)
814 :     in
815 :     E.EIN{
816 :     params = [E.IMG(dim,shape), E.KRN],
817 :     index = shape,
818 :     body = E.Conv(0,expindex, 1, [])
819 :     }
820 :     end
821 : jhr 3476
822 : jhr 3477 (* Probe: <F(x)>_{\alpha} *)
823 :     fun probe (alpha, dim) = let
824 :     val expindex = specialize(alpha, 0)
825 :     in
826 :     E.EIN{
827 :     params = [E.FLD dim, E.TEN(0, [])], index = alpha,
828 :     body = E.Probe(E.Field(0, expindex), E.Tensor(1, []))
829 :     }
830 :     end
831 : jhr 3476
832 : jhr 3477 (***************************** derivative ****************************)
833 : jhr 3476
834 : jhr 3477 (* \EinExp{\sum_{ij}\mathcal{E}_{ij} \frac{F_j}{\partial x_i} *)
835 :     val curl2d = E.EIN{
836 :     params = [E.FLD 2],
837 :     index = [],
838 :     body = E.Sum([(E.V 0, 0, 1), (E.V 1, 0, 1)],
839 :     E.Opn(E.Prod, [
840 :     E.G(E.Eps2(0, 1)),
841 :     E.Apply(E.Partial[E.V 0], E.Field(0, [E.V 1]))
842 :     ]))
843 :     }
844 : jhr 3476
845 : jhr 3477 val curl3d = E.EIN{
846 :     params = [E.TEN(1, [3])],
847 :     index = [3],
848 :     body = E.Sum([(E.V 1, 0, 2), (E.V 2, 0, 2)],
849 :     E.Opn(E.Prod, [
850 :     E.G(E.Epsilon(0, 1, 2)),
851 :     E.Apply(E.Partial[E.V 1], E.Field(0, [E.V 2]))
852 :     ]))
853 :     }
854 : jhr 3476
855 : jhr 3477 (*< d F / d_i>_i *)
856 :     fun grad (alpha as a::_) = let
857 :     val expindex = specialize(alpha, 0)
858 :     in
859 :     E.EIN{
860 :     params = [E.FLD a],
861 :     index = alpha,
862 :     body = E.Apply(E.Partial expindex, E.Field(0, []))
863 :     }
864 :     end
865 :    
866 :     (*< Sigma d F_alpha / d x_i>ALpha i CHANGE HERE *)
867 :     fun dotimes (dim, alpha) = let
868 :     val n = length alpha
869 :     (*
870 :     fun expIndex(n, inc) = List.tabulate(n, (fn(x)=>E.V (x+inc)))
871 :     val i' = expIndex(n, 0)
872 :     *)
873 :     val i' = List.tabulate (n, fn x => E.V x)
874 :     in
875 :     E.EIN{
876 :     params = [E.FLD dim], index = alpha@[dim],
877 :     body = E.Apply(E.Partial[E.V n], E.Field(0, i'))
878 :     }
879 :     end
880 :    
881 :     (* <d F_i /d_i> *)
882 :     fun divergence (dim, alpha) = let
883 :     val expindex = specialize(alpha, 0)
884 :     val sumI = length alpha
885 :     val sumIndex = E.V sumI
886 :     val sumIndexL = [sumIndex]
887 :     val S = expindex@sumIndexL
888 :     in
889 :     E.EIN{
890 :     params = [E.FLD dim],
891 :     index = alpha,
892 :     body = E.Sum([(sumIndex, 0, dim-1)], E.Apply(E.Partial sumIndexL, E.Field(0, S)))
893 :     }
894 :     end
895 :    
896 : jhr 3476 end (* mkOperators *)

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0