;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                ctensor.lisp                   ;;;
;;;          numeical tensor operation            ;;;
;;;          (c) Kenrou, Adachi, 2007/2/8-        ;;;
;;;          vwesion 1.0.1 (2007/3/9)             ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; Utility from the book ''ANSI Common Lisp'' by Paul Graham ;;;

(defmacro for (var start stop &body body)
  (let ((gstop (gensym)))
    `(do ((,var ,start (1+ ,var))
	  (,gstop ,stop))
	 ((> ,var ,gstop))
	 ,@body)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;; tensortype , tensorrank, zerotensor, tensorp ;;;

(defun pseudotype (a)
  (labels
   ((rectyp (typ tmp)
	    (let ((ftmp (car tmp)))
	      (if (listp ftmp)
		  (rectyp (cons (length ftmp) typ) ftmp)
		typ))))
   (reverse (rectyp (cons (length a) nil) a))))


(defun tointeger (a)
  (if (not (listp (car a)))
      (mapcar #'round a)
    (mapcar #'(lambda (x) (tointeger x)) a)))


(defun tensortype (a)     ;depending on zerotensor below 
  (labels
   ((sp (k c)
	(if (not (listp (car c)))
	    (mapcar #'(lambda (x) (* k x)) c)
	  (mapcar #'(lambda (x) (sp k x)) c))))
  (let*
      ((typ (pseudotype a))
       (za (zerotensor typ))
       (ta (tointeger (sp 0 a))))
    (if (equal za ta)
	typ
      nil))))
      

(defun tensorrank (f)
  (length (tensortype f)))


(defun zerotensor (ttype)
  (let*
      ((n (length ttype))
       (rtyp (reverse ttype))
       (rlst (car rtyp))
       (zero nil))
    (labels
     ((zerolist (m)
		(if (= m 0)
		    nil
		  (cons 0 (zerolist (1- m)))))
      (rec (m rz)
	   (if (= m 0)
	       nil
	     (cons rz (rec (1- m) rz)))))
     (for i 1 n
	  (if (= i 1)
	      (progn
		(setf zero (zerolist rlst))
		(setf rtyp (cdr rtyp)))
	    (progn
	      (setf rlst (car rtyp))
	      (setf zero (rec rlst zero))
	      (setf rtyp (cdr rtyp)))))
     zero)))


(defun tensorp (a)
  (labels
   ((sp (k c)
	(if (not (listp (car c)))
	    (mapcar #'(lambda (x) (* k x)) c)
	  (mapcar #'(lambda (x) (sp k x)) c))))
   (let ((za (zerotensor (pseudotype a)))
	 (ta (tointeger (sp 0 a))))
     (if (equal za ta)
	 t
       nil))))


(defun sametensortypep (a b)
  (labels
   ((sp (k c)
	(if (not (listp (car c)))
	    (mapcar #'(lambda (x) (* k x)) c)
	  (mapcar #'(lambda (x) (sp k x)) c))))
   (if (and (tensorp a) (tensorp b))
       (let ((za (tointeger (sp 0 a))) (zb (tointeger (sp 0 b))))
	 (if (equal za zb)
	     t
	   nil))
     nil)))


;;; tensorproduct, scalarproduct ;;;


(defun sclrprdct (k b)
  (if (numberp (car b))
      (mapcar #'(lambda (x) (* k x)) b)
    (mapcar #'(lambda (x) (sclrprdct k x)) b)))


(defun pseudoproduct (a b)
  (if (not (listp (car a)))
      (mapcar #'(lambda (x) (sclrprdct x b)) a)
    (mapcar #'(lambda (x) (pseudoproduct x b)) a)))


(defun tensorproduct (a b)
  (if (and (tensorp a) (tensorp b))
      (pseudoproduct a b)
    nil))


(defun scalarproduct (k b)
  (if (and (numberp k) (tensorp b))
      (sclrprdct k b)
    nil))


;;; tensorplus, tensorminus ;;;

(defun tpls (a b)
  (if (numberp (car a))
      (mapcar #'+ a b)
    (mapcar #'tpls a b)))


(defun tensorplus (a b)
  (if (sametensortypep a b)
      (tpls a b)
    nil))


(defun tensorminus (a b)
  (tensorplus a (sclrprdct -1 b)))


;;; exhangetensortype ;;;

(defun endcons (x a)
   (reverse (cons x (reverse a))))


(defun exchinext (a)
  (let
      ((ni  (length a))
       (nip (length (car a)))
       (b   nil))
    (for j 1 nip
	 (let ((c nil))
	   (for k 1 ni
		(setf c (endcons (nth (1- j) (nth (1- k) a)) c)))
	   (setf b (endcons c b))))
    b))


(defun exchnext (a i)
  (let ((n (tensorrank a)))
    (if (or (> i (1- n)) (< i 1))
	nil
      (if (= i 1)
	  (exchinext a)
        (mapcar #'(lambda (x) (exchnext x (1- i))) a)))))


(defun exchangetensortype (a i j)
  (if (not (tensorp a))
      nil
    (if (= i j)
	nil
      (if (> i j)
	  (exchangetensortype a j i)
	(if (= i (1- j))
            (exchnext a i)
	  (let ((b a))
	    (for k i (1- j)
		 (setf b (exchnext b k)))
	    (for k i (- j 2)
		 (setf b (exchnext b (- (+ (- j 2) i) k))))
	    b))))))


;;; ,squarematrixp, matrixtrace, contract ;;;

(defun squarematrixp (a)
  (let ((trk (tensorrank a)))
    (if (not (= 2 trk))
	nil
      (if (not (= (length a) (length (car a))))
	  nil
	t))))


(defun matrixtrace (a)
  (if (not (squarematrixp a))
      nil
    (let ((n (length a)) (trace 0))
      (for i 1 n
	   (setf trace (+ trace (nth (1- i) (nth (1- i) a)))))
      trace)))


(defun cntrlst2 (a)
  (if (squarematrixp a)
      (matrixtrace a)
    (mapcar #'(lambda (x) (cntrlst2 x)) a)))


(defun contract (a i j)
  (if (not (tensorp a))
      nil
    (let ((n (tensorrank a)))
      (if (or (= i j) (< i 1) (> i n) (< j 1) (> j n))
	  nil
	(if (> i j)
	    (contract a j i)
	  (if (= i (1- n))
	      (cntrlst2 a)
	    (let ((b a))
	      (if (not (= j n))
		  (for p j (1- n)
		       (setf b (exchnext b p)))
		nil)
	      (if (not (= i (1- n)))
		  (for p i (- n 2)
		       (setf b (exchnext b p)))
		nil)
	      (setf b (cntrlst2 b))
	      b)))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                matrix operation                      ;;;
;;; matrixp, matrixproduct, matrixtranspose, matrixpower ;;;
;;; identitymatrix, exponential, expmap                  ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun matrixp (a)
  (if (= (tensorrank a) 2)
      t
    nil))


(defun matrixproduct (a b)
  (if (or (not (matrixp a)) (not (matrixp b)))
      nil
    (if (not (= (second (tensortype a)) (first (tensortype b))))
	nil
      (let ((mat (tensorproduct a b)))
	(contract mat 2 3)))))


(defun matrixtranspose (a)
  (if (not (matrixp a))
      nil
    (exchinext a)))


(defun matrixpower (a n)
  (if (or (not (integerp n)) (< n 0) (not (squarematrixp a)))
      nil
    (if (= n 0)
	(identitymatrix (car (tensortype a)))
      (labels
       ((mp (a n)
	    (if (= n 1)
		a
	      (matrixproduct a (mp a (1- n))))))
       (mp a n)))))


(defun identitymatrix (n)
  (let ((id nil) (icol nil))
    (for i 0 (1- n)
	 (setf icol nil)
	 (for j 0 (1- n)
	      (if (= i j)
		  (setf icol (cons 1 icol))
		(setf icol (cons 0 icol))))
	 (setf id (cons icol id)))
    id))


(defun factorial (n)
  (if (= n 0)
      1
    (* n (factorial (1- n)))))


(setq *degree-of-expansion* 20)


(defun deg-exp (m)
  (setf *degree-of-expansion* m))


(defun exponential (a s)
  (if (numberp s)
      (let ((expta (zerotensor (tensortype a))))
	(for i 0 *degree-of-expansion*
	     (setf expta 
		   (tpls expta
			 (sclrprdct (/ 1 (factorial i))
				    (matrixpower (sclrprdct s a) i)))))
	expta)
    nil))


(defun exponential2 (a s)
  (if (numberp s)
      (let*
	  ((expta (zerotensor (tensortype a)))
	   (m (car (tensortype a)))
	   (idm (identitymatrix m)))
	(for i 0 (1- *degree-of-expansion*)
	     (setf expta
		   (matrixproduct a
			     (tpls
			      (sclrprdct
			       (/ 1 (factorial (- *degree-of-expansion* i)))
			       idm)
			      expta))))
	(setf expta (tpls idm expta))
	expta)))


(defun expmap (a)
  (exponential a 1))


;;; matrixplus, matrixminus ;;;

(defun matrixplus (a b)
  (if (and (matrixp a) (sametensortypep a b))
      (tensorplus a b)
    nil))


(defun matrixminus (a b)
  (matrixplus a (sclrprdct -1 b)))


;;; alias ;;;

(defun sprod (k a)
  (scalarproduct k a))


(defun tprod (a b)
  (tensorproduct a b))


(defun tplus (a b)
  (tensorplus a b))


(defun tminus (a b)
  (tensorminus a b))


(defun matprod (a b)
  (matrixproduct a b))


(defun matplus (a b)
  (matrixplus a b))


(defun matminus (a b)
  (matrixminus a b))


(defun mattrace (a)
  (matrixtrace a))


(defun matexp (a)
  (expmap a))


(defun matid (n)
  (identitymatrix n))


(defun mpow (a n)
  (matrixpower a n))

;;; for compatibility to old version ;;;

(defun sproduct (k a)
  (scalarproduct k a))


(defun tproduct (a b)
  (tensorproduct a b))


(defun mproduct (a b)
  (matrixproduct a b))


(defun exchttype (a i j)
  (exchangetensortype a i j))
 

;;; end of ctensor.lisp ;;;