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

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