UNPKG

1.88 kBtext/coffeescriptView Raw
1### contract =====================================================================
2
3Tags
4----
5scripting, JS, internal, treenode, general concept
6
7Parameters
8----------
9a,i,j
10
11General description
12-------------------
13Contract across tensor indices i.e. returns "a" summed over indices i and j.
14If i and j are omitted then 1 and 2 are used.
15contract(m) is equivalent to the trace of matrix m.
16
17###
18
19
20
21Eval_contract = ->
22 push(cadr(p1))
23 Eval()
24 if (cddr(p1) == symbol(NIL))
25 push_integer(1)
26 push_integer(2)
27 else
28 push(caddr(p1))
29 Eval()
30 push(cadddr(p1))
31 Eval()
32 contract()
33
34contract = ->
35 save()
36 yycontract()
37 restore()
38
39yycontract = ->
40 h = 0
41 i = 0
42 j = 0
43 k = 0
44 l = 0
45 m = 0
46 n = 0
47 ndim = 0
48 nelem = 0
49 ai = []
50 an = []
51
52 p3 = pop()
53 p2 = pop()
54 p1 = pop()
55
56 if (!istensor(p1))
57 if (!iszero(p1))
58 stop("contract: tensor expected, 1st arg is not a tensor")
59 push(zero)
60 return
61
62 push(p2)
63 l = pop_integer()
64
65 push(p3)
66 m = pop_integer()
67
68 ndim = p1.tensor.ndim
69
70 if (l < 1 || l > ndim || m < 1 || m > ndim || l == m \
71 || p1.tensor.dim[l - 1] != p1.tensor.dim[m - 1])
72 stop("contract: index out of range")
73
74 l--
75 m--
76
77 n = p1.tensor.dim[l]
78
79 # nelem is the number of elements in "b"
80
81 nelem = 1
82 for i in [0...ndim]
83 if (i != l && i != m)
84 nelem *= p1.tensor.dim[i]
85
86 p2 = alloc_tensor(nelem)
87
88 p2.tensor.ndim = ndim - 2
89
90 j = 0
91 for i in [0...ndim]
92 if (i != l && i != m)
93 p2.tensor.dim[j++] = p1.tensor.dim[i]
94
95 a = p1.tensor.elem
96 b = p2.tensor.elem
97
98 for i in [0...ndim]
99 ai[i] = 0
100 an[i] = p1.tensor.dim[i]
101
102 for i in [0...nelem]
103 push(zero)
104 for j in [0...n]
105 ai[l] = j
106 ai[m] = j
107 h = 0
108 for k in [0...ndim]
109 h = (h * an[k]) + ai[k]
110 push(a[h])
111 add()
112 b[i] = pop()
113 for j in [(ndim - 1)..0]
114 if (j == l || j == m)
115 continue
116 if (++ai[j] < an[j])
117 break
118 ai[j] = 0
119
120 if (nelem == 1)
121 push(b[0])
122 else
123 push(p2)
124