135 lines
2.2 KiB
CoffeeScript
135 lines
2.2 KiB
CoffeeScript
### contract =====================================================================
|
|
|
|
Tags
|
|
----
|
|
scripting, JS, internal, treenode, general concept
|
|
|
|
Parameters
|
|
----------
|
|
a,i,j
|
|
|
|
General description
|
|
-------------------
|
|
Contract across tensor indices i.e. returns "a" summed over indices i and j.
|
|
If i and j are omitted then 1 and 2 are used.
|
|
contract(m) is equivalent to the trace of matrix m.
|
|
|
|
###
|
|
|
|
|
|
|
|
Eval_contract = ->
|
|
push(cadr(p1))
|
|
Eval()
|
|
if (cddr(p1) == symbol(NIL))
|
|
push_integer(1)
|
|
push_integer(2)
|
|
else
|
|
push(caddr(p1))
|
|
Eval()
|
|
push(cadddr(p1))
|
|
Eval()
|
|
contract()
|
|
|
|
contract = ->
|
|
save()
|
|
yycontract()
|
|
restore()
|
|
|
|
yycontract = ->
|
|
h = 0
|
|
i = 0
|
|
j = 0
|
|
k = 0
|
|
l = 0
|
|
m = 0
|
|
n = 0
|
|
ndim = 0
|
|
nelem = 0
|
|
ai = []
|
|
an = []
|
|
|
|
p3 = pop()
|
|
p2 = pop()
|
|
p1 = pop()
|
|
|
|
if (!istensor(p1))
|
|
if (!isZeroAtomOrTensor(p1))
|
|
stop("contract: tensor expected, 1st arg is not a tensor")
|
|
push(zero)
|
|
return
|
|
|
|
push(p2)
|
|
l = pop_integer()
|
|
|
|
push(p3)
|
|
m = pop_integer()
|
|
|
|
ndim = p1.tensor.ndim
|
|
|
|
if (l < 1 || l > ndim || m < 1 || m > ndim || l == m \
|
|
|| p1.tensor.dim[l - 1] != p1.tensor.dim[m - 1])
|
|
stop("contract: index out of range")
|
|
|
|
l--
|
|
m--
|
|
|
|
n = p1.tensor.dim[l]
|
|
|
|
# nelem is the number of elements in "b"
|
|
|
|
nelem = 1
|
|
for i in [0...ndim]
|
|
if (i != l && i != m)
|
|
nelem *= p1.tensor.dim[i]
|
|
|
|
#console.log "nelem:" + nelem
|
|
p2 = alloc_tensor(nelem)
|
|
#console.log "p2:" + p2
|
|
|
|
p2.tensor.ndim = ndim - 2
|
|
|
|
j = 0
|
|
for i in [0...ndim]
|
|
if (i != l && i != m)
|
|
p2.tensor.dim[j++] = p1.tensor.dim[i]
|
|
|
|
|
|
a = p1.tensor.elem
|
|
b = p2.tensor.elem
|
|
|
|
#console.log "a: " + a
|
|
#console.log "b: " + b
|
|
|
|
for i in [0...ndim]
|
|
ai[i] = 0
|
|
an[i] = p1.tensor.dim[i]
|
|
|
|
for i in [0...nelem]
|
|
push(zero)
|
|
for j in [0...n]
|
|
ai[l] = j
|
|
ai[m] = j
|
|
h = 0
|
|
for k in [0...ndim]
|
|
h = (h * an[k]) + ai[k]
|
|
push(a[h])
|
|
#console.log "a[h]: " + a[h]
|
|
add()
|
|
#console.log "tos: " + stack[tos-1]
|
|
b[i] = pop()
|
|
#console.log "b[i]: " + b[i]
|
|
for j in [(ndim - 1)..0]
|
|
if (j == l || j == m)
|
|
continue
|
|
if (++ai[j] < an[j])
|
|
break
|
|
ai[j] = 0
|
|
|
|
if (nelem == 1)
|
|
push(b[0])
|
|
else
|
|
push(p2)
|
|
|
|
#console.log "returning: " + stack[tos-1]
|