;   ctensor.lisp
;   numerical tensor operation
;   (c) Kenrou, Adachi, 2007/2/8, 2007/3/6
;   version 0.3

;;; 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)))

(defun copylist (lst)
  (if (atom lst)
    lst
    (cons (car lst) (copylist (cdr lst)))))

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


;;; type of tensor ;;;

(defun tensortype (a)
  (labels
   ((rectyp (typ tmp)
	    (progn
	      (setf tmp (car tmp))
	      (if (listp tmp)
		  (progn
		    (setf typ (cons (length tmp) typ))
		    (rectyp typ tmp))
		typ))))
   (let* ((ttmp a) (ttyp (cons (length ttmp) nil)))
     (reverse (rectyp ttyp ttmp)))))


;;; rank of tensor ;;;

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


;;; make a zero tensor ;;;

(defun zerotensor (ttype)
  (let*
    ((rr (length ttype))
     (tt1 nil)
     (rtyp (reverse ttype))
     (rk (car rtyp))
     (rst (cdr rtyp)))

    (for i 1 rk
      (setf tt1 (cons 0 tt1)))
    (if (not (eq rst nil))
      (progn
        (setf rk (car rst))
        (for j 2 rr
          (let ((tt2 nil))
            (for i 1 rk
              (setf tt2 (cons tt1 tt2)))
            (setf tt1 tt2)
            (setf rst (cdr rst))
            (if (not (eq rst nil))
              (setf rk (car rst))
              nil))))
        nil)
     tt1))

;;; another zero tensor
;;; こっちの方がきれい
;;; 最後に残った for を再帰に書けば、
;;; scheme にも移植できるかも。

(defun zerotensor2 (ttype)
  (let*
      ((n (length ttype))
       (rtyp (reverse ttype))
       (rlast (car rtyp))
       (zerot 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 zerot (zerolist rlast))
	       (setf rtyp (cdr rtyp)))
	     (progn
	       (setf rlast (car rtyp))
	       (setf zerot (rec rlast zerot))
	       (setf rtyp (cdr rtyp)))))

      zerot)))

;;; tensor product ;;;

;;; Warning!,    KCL, AKCL, GCL on 4.3BSD Qusijarus0c           ;;;
;;; the mtimes below causes "Invocation history stack overflow" ;;;

(defun mtimes (k b)
   (if (not (listp (car b)))
     (mapcar #'(lambda (x) (* k x)) b)
     (mapcar #'(lambda (x) (mtimes k x)) b)))


(defun tproduct (a b)
  (if (not (listp (car a)))
    (mapcar #'(lambda (x) (mtimes x b)) a)
    (mapcar #'(lambda (x) (tproduct x b)) a)))


;;; on 4.3BSD Quaijarus on simh vax simulator ;;;
;;;     must use the following definition     ;;;
;;                                             ;;
;; (defun mtimes (k b)                         ;;
;;   (labels                                   ;;
;;     ((sp (x) (* k x))                       ;;
;;      (mt (x) (mtimes k x)))                 ;;
;;                                             ;;
;;     (if (not (listp (car b)))               ;;
;;       (mapcar #'sp b)                       ;;
;;       (mapcar #'mt b))))                    ;;
;;                                             ;;
;; (defun tproduct (a b)                       ;;
;;   (labels                                   ;;
;;     ((mp (x) (mtimes x b))                  ;;
;;      (tp (x) (tproduct x b)))               ;;
;;                                             ;;
;;   (if (not (listp (car a)))                 ;;
;;     (mapcar #'mp a)                         ;;
;;     (mapcar #'tp a))))                      ;;
;;                                             ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;; scalar product ;;;

(defun sproduct (k b)
  (mtimes k b))


;;; exchange tensor type ;;;

;;;;;;;;;;; 効率悪し ;;;;;;;;;;;;;;;;;;

(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 exchttype (a i j)
   (if (= i j)
      nil
      (if (> i j)
         (exchttype a j i)
         (if (= i (1- j))
            (exchnext a i)
            (let ((b (copylist 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)))))


;;; contraction of tensor ;;;

(defun sqmatp (a)
  (if (not (numberp (car (car a))))
    nil
    (if (not (= (length a) (length (car a))))
      nil
      t)))


(defun mattrace (a)
  (if (not (sqmatp a))
    nil
    (let ((n (length a)) (tra 0))
      (for i 1 n
        (setf tra (+ tra (nth (1- i) (nth (1- i) a)))))
      tra)))


(defun cntrlst2 (a)
  (if (sqmatp a)
    (mattrace a)
    (mapcar #'(lambda (x) (cntrlst2 x)) a)))


(defun contract (a i j)
  (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 (copylist 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))))))


;;; plus function of two tensors ;;;
;; 2007/2/9 ;;

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


(defun tplus (a b)
  (if (not (equal (tensortype a) (tensortype b)))
      nil
    (tpls a b)))


(defun tminus (a b)
  (tplus a (sproduct -1 b)))


;; matrix product ;;
;; matrix = rank2 tensor ;;
;; 2007/2/12 ;;
;; 以下は参考としての実装であり、実際には非効率なので使ってはいけない。
;; (電卓程度なら使用可)

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


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


(defun mtp (a b)
  (mproduct a b))


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


(defun mpow (a n)
  (if (or (not (integerp n)) (< n 0) (not (sqmatp a)))
      nil
    (if (= n 0)
	(identitymatrix (car (tensortype a)))
      (labels
       ((mp (a n)
	    (if (= n 1)
		a
	      (mtp a (mp a (1- n))))))
       (mp a n)))))

;; 2007/3/3

(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))

;; utility ;;

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; exp(sA) ~ I + (1/1!)s^1A^1 + (1/2!)s^2A^2 + ... + (1/n!)s^nA^n  ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun exponential (a s n)
  (let ((expta (zerotensor (tensortype a))))
    (for i 0 n
	 (setf expta 
	       (tplus expta
		      (sproduct (/ 1 (factorial i)) (mpow (sproduct s a) i)))))
    expta))
 

;;;;;; 効率はこっちの方がいいはず ;;;;;;;

(defun exponential2 (a s n)
  (let*
      ((expta (zerotensor (tensortype a)))
       (m (car (tensortype a)))
       (idm (identitymatrix m)))
    (for i 0 (1- n)
	 (setf expta
	       (mproduct a
			 (tplus
			  (sproduct (/ 1 (factorial (- n i))) idm)
			  expta))))
    (setf expta (tplus idm expta))
    expta))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; exp : Mat(n) ---> GL(n) ; exponential map ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(setq *degree-of-expansion* 20)


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


(defun expmap (a)
  (exponential a 1 *degree-of-expansion*)) 



;; 2007/3/6
;; tensorp テンソルかどうか ;;
;; 以下、手抜きの実装なので、完全には動作しない。;;

(defun allnumberp (a)
  (let ((ta (copylist a))
	(ca nil)
	(n (length a))
	(bool 1))
    (for i 1 n
	 (progn
	   (setf ca (car ta))
	   (setf ta (cdr ta))
	   (setf bool
		 (if (listp ca)
		     0
		   (* bool 1)))))
    (if (= bool 1)
	t
      nil)))
	       

;;;;;; a がテンソルでないときに、少し調べてみるのに使う。;;;;;;
;                                                      ;
;(defun subttypelist (a)                               ;
;  (if (not (equal a nil))                             ;
;      (let ((stp nil) (ra nil))                       ;
;	(setf stp (cons (tensortype (car a)) stp))     ;
;	(setf ra (cdr a))                              ;
;	(if (not (equal ra nil))                       ;
;	    (setf stp (append stp (subttypelist ra)))  ;
;	  stp))                                        ;
;    nil))                                             ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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


(defun sametensortypep (a b)
  (if (or (not (listp a)) (not (listp b)))
      nil
    (let ((ta (tointeger (sproduct 0 a)))
	  (tb (tointeger (sproduct 0 b))))
      (if (equal ta tb)
	  t
	nil))))


(defun samettype (a b)
  (if (or (not (listp a)) (not (listp b)))
      0
    (let ((ta (tointeger (sproduct 0 a))) 
	  (tb (tointeger (sproduct 0 tb))))
      (if (equal ta tb)
	  1
	0))))


;; obsolete ;;

(defun old-tensorp (a)
  (if (not (listp a))
      nil
    (if (allnumberp a)
	t
      (let ((bool 1) (ma (car a)) (ra (cdr a)) (n (length a)))
	(for i 1 (1- n)
	     (progn
	       (setf bool (* bool (samettype ma (car ra))))
	       (setf ra (cdr ra))))
	(if (= bool 1)
	    t
	  nil)))))


;; new-version

(defun tensorp (a)
  (let ((za (zerotensor (tensortype a)))
	(ta (tointeger (sproduct 0 a))))
    (if (equal za ta)
	t
      nil)))


;;; end of ctensor.lisp ;;;