home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 3 / TheARMClub_PDCD3.iso / hensa / maths / b116_1 / jacal / tensor < prev    next >
Text File  |  1993-12-21  |  7KB  |  269 lines
  1. ;;; TENSOR.SCM -- Tensor-like support functions for JACAL
  2. ;;; Copyright (C) 1993 Jerry D. Hedden
  3. ;;; See the file `COPYING' for terms applying to this program.
  4.  
  5. ; Does not support the notion of contra-/covariant indices.
  6. ;  Users must keep track of this information themselves.
  7.  
  8. ; Assumes that all matrices are "proper" (i.e., that all "dimensions"
  9. ;  of the matrix are the same length, e.g., 4x4x4) and "compatible"
  10. ;  (e.g., a 3x3 matrix is not compatible with a 4x4x4 matrix).
  11.  
  12.  
  13. (definfo 'indexshift
  14.   "Shifts an index within a tensor")
  15.  
  16. (definfo 'indexswap
  17.   "Swaps two indices within a tensor")
  18.  
  19. (definfo 'contract
  20.   "Tensor contraction")
  21.  
  22. (definfo 'tmult
  23.   "Tensor multiplication")
  24.  
  25.  
  26. ; helper function to determine the "rank" of a tensor
  27. (define (tnsr:rank m)
  28.   (let rnk ((rank  0)
  29.         (mm    m))
  30.     (if (bunch? mm)
  31.       (rnk (+ rank 1) (car mm))
  32.       rank)))
  33.  
  34.  
  35. (defbltn 'indexshift
  36.   (lambda (m . args)
  37.     (let ((rank  (tnsr:rank m)))
  38.       (cond
  39.     ((= rank 0)
  40.        ; scalar -- no-op
  41.        m)
  42.     ((= rank 1)
  43.        ; vector -- returns a pseudo-tensor
  44.        (map list m))
  45.     ((or (= rank 2) (null? args))
  46.        ; matrix -- transpose
  47.        (apply map list m))
  48.     (else
  49.        ; tensor
  50.        ; set and constrain the "from" and "to" positions
  51.        (let* ((a  (car args))
  52.           (b  (if (null? (cdr args))
  53.             (+ a 1)
  54.             (cadr args))))
  55.          (if (< a 1)
  56.            (set! a 1)
  57.            (if (> a rank)
  58.          (set! a rank)))
  59.          (if (< b 1)
  60.            (set! b 1)
  61.            (if (> b rank)
  62.          (set! b rank)))
  63.          (if (= a b)
  64.            (if (= a rank)
  65.          (set! a (- b 1))
  66.          (set! b (+ a 1))))
  67.  
  68.          (if (< a b)
  69.            ; index shift right
  70.            (let isr1 ((ma  m)
  71.               (nn  1))
  72.          (if (< nn a)
  73.            (map (lambda (mm) (isr1 mm (+ nn 1))) ma)
  74.            (let isr2 ((mb  ma)
  75.                   (aa  (+ a 1)))
  76.              (if (= aa b)
  77.                (apply map list mb)
  78.                (map (lambda (mm) (isr2 mm (+ aa 1)))
  79.                 (apply map list mb))))))
  80.            ; index shift left
  81.            (let isl1 ((ma  m)
  82.               (nn  1))
  83.          (if (< nn b)
  84.            (map (lambda (mm) (isl1 mm (+ nn 1))) ma)
  85.            (let isl2 ((mb  ma)
  86.                   (bb  (+ b 1)))
  87.              (if (= bb a)
  88.                (apply map list mb)
  89.                (apply map list (map (lambda (mm) (isl2 mm (+ bb 1)))
  90.                         mb)))))))))))))
  91.  
  92.  
  93. (defbltn 'indexswap
  94.   (lambda (m . args)
  95.     (let ((rank  (tnsr:rank m)))
  96.       (cond
  97.     ((= rank 0)
  98.        ; scalar -- no-op
  99.        m)
  100.     ((= rank 1)
  101.        ; vector -- returns a pseudo-tensor
  102.        (map list m))
  103.     ((or (= rank 2) (null? args))
  104.        ; matrix -- transpose
  105.        (apply map list m))
  106.     (else
  107.        ; tensor
  108.        ; set and constrain the indices to be swapped
  109.        (let* ((a  (car args))
  110.           (b  (if (null? (cdr args))
  111.             (+ a 1)
  112.             (cadr args))))
  113.          (if (< b a)
  114.            (let ((c  a))
  115.          (begin (set! a b)
  116.             (set! b c))))
  117.          (if (< a 1)
  118.            (begin (set! a 1)
  119.               (if (<= b a)
  120.             (set! b 2)))
  121.            (if (> b rank)
  122.          (begin (set! b rank)
  123.             (if (<= b a)
  124.               (set! a (- b 1))))
  125.          (if (= a b)
  126.            (if (= a rank)
  127.              (set! a (- b 1))
  128.              (set! b (+ a 1))))))
  129.  
  130.          ; perform the swapping operation
  131.          (let swap1 ((ma  m)
  132.              (n   1))
  133.            (if (< n a)
  134.            (map (lambda (mm) (swap1 mm (+ n 1))) ma)
  135.            (let swap2 ((mb  ma)
  136.                    (aa  (+ a 1)))
  137.              (if (= aa b)
  138.              (apply map list mb)
  139.              (apply map list
  140.                     (map (lambda (mm) (swap2 mm (+ aa 1)))
  141.                      (apply map list mb)))))))))))))
  142.  
  143.  
  144. ; helper function for the contraction operation
  145. ;  sums the diagonal elements of a matrix
  146. (define (tnsr:contract m)
  147.   (let cxt ((mm  (map cdr (cdr m)))
  148.         (ss  (car (car m))))
  149.     (if (null? mm)
  150.       ss
  151.       (cxt (map cdr (cdr mm)) (app* $1+$2 ss (car (car mm)))))))
  152.  
  153.  
  154. (defbltn 'contract
  155.   (lambda (m . args)
  156.     (let ((rank  (tnsr:rank m)))
  157.       (cond
  158.     ((= rank 0)
  159.        ; scalar -- no-op
  160.        m)
  161.     ((= rank 1)
  162.        ; vector -- sum elements
  163.        (reduce (lambda (x y) (app* $1+$2 x y)) m))
  164.     ((= rank 2)
  165.        ; matrix -- sum diagonal elements
  166.        (tnsr:contract m))
  167.     (else
  168.        ; tensor
  169.        ; set and constrain the indices for the contraction operation
  170.        (let* ((a  (car args))
  171.           (b  (if (null? (cdr args))
  172.             (+ a 1)
  173.             (cadr args))))
  174.          (if (< b a)
  175.            (let ((c  a))
  176.          (begin (set! a b)
  177.             (set! b c))))
  178.          (if (< a 1)
  179.            (begin (set! a 1)
  180.               (if (<= b a)
  181.             (set! b 2)))
  182.            (if (> b rank)
  183.          (begin (set! b rank)
  184.             (if (<= b a)
  185.               (set! a (- b 1))))
  186.          (if (= a b)
  187.            (if (= a rank)
  188.              (set! a (- b 1))
  189.              (set! b (+ a 1))))))
  190.  
  191.        ; perform the contraction operation
  192.        (let cxt1 ((ma  m)
  193.               (nn  1))
  194.          (if (< nn a)
  195.            (map (lambda (mm) (cxt1 mm (+ nn 1))) ma)
  196.            (let cxt2 ((mb  ma)
  197.               (aa  (+ a 1)))
  198.          (if (< aa b)
  199.            (map (lambda (mm) (cxt2 mm (+ aa 1))) (apply map list mb))
  200.            (let cxt3 ((mc  mb)
  201.                   (bb  b))
  202.              (if (= bb rank)
  203.                (tnsr:contract mc)
  204.                (map (lambda (mm) (cxt3 mm (+ bb 1)))
  205.                 (apply map list (map (lambda (mx)
  206.                            (apply map list mx))
  207.                          mc)))))))))))))))
  208.  
  209.  
  210. (defbltn 'tmult
  211.   (lambda (m1 m2 . args)
  212.     (let ((r1  (tnsr:rank m1))
  213.       (r2  (tnsr:rank m2)))
  214.       (cond
  215.     ((or (= r1 0) (= r2 0))
  216.        ; scalar multiplication
  217.        (app* $1*$2 m1 m2))
  218.     ((null? args)
  219.        ; outerproduct -- scalar multiplication of the second
  220.        ;  tensor by each element of the first tensor
  221.        (let outerproduct ((ma  m1)
  222.                   (r   1))
  223.          (if (< r r1)
  224.            (map (lambda (mm) (outerproduct mm (+ r 1))) ma)
  225.            (map (lambda (x) (app* $1*$2 x m2)) ma))))
  226.     (else
  227.        ; innerproduct
  228.        ; set and contrain indices to be used
  229.        (let* ((a  (car args))
  230.           (b  (if (null? (cdr args))
  231.             a
  232.             (cadr args))))
  233.          (if (< a 1)
  234.            (set! a 1)
  235.            (if (> a r1)
  236.          (set! a r1)))
  237.          (if (< b 1)
  238.            (set! b 1)
  239.            (if (> b r2)
  240.          (set! b r2)))
  241.  
  242.          ; perform the multiplication operation
  243.          (let mult1 ((ma1  m1)
  244.              (n1   1))
  245.            (if (< n1 a)
  246.          ; find index to multiply in first tensor
  247.          (map (lambda (mm) (mult1 mm (+ n1 1))) ma1)
  248.          (let mult2 ((mb1  ma1)
  249.                  (a1   a))
  250.            (if (< a1 r1)
  251.              ; shift index to last position in first tensor
  252.              (map (lambda (mm) (mult2 mm (+ a1 1)))
  253.               (apply map list mb1))
  254.              (let mult3 ((ma2  m2)
  255.                  (n2   1))
  256.                (if (< n2 b)
  257.              ; find index to multiply in second tensor
  258.              (map (lambda (mm) (mult3 mm (+ n2 1))) ma2)
  259.              (let mult4 ((mb2  ma2)
  260.                      (a2   b))
  261.                (if (< a2 r2)
  262.                  ; shift index to last position in second tensor
  263.                  (map (lambda (mm) (mult4 mm (+ a2 1)))
  264.                   (apply map list mb2))
  265.                  ; the actual multiplication is done here
  266.                  (reduce (lambda (x y) (app* $1+$2 x y))
  267.                      (map (lambda (x y) (app* $1*$2 x y))
  268.                       mb1 mb2))))))))))))))))
  269.