tebdol

tensor.lisp (21817B)

---

 (defpackage :tensor
   (:use :common-lisp :util :blas :array :exp)
   (:export :tensor
 	   :make-tensor
 	   :tensor-indices
 	   :tensor-sectors
 	   :segment
 	   :make-segment
 	   :segment-numbers
 	   :segment-dimension
 	   :subscript
 	   :make-subscript
 	   :subscript-numbers
 	   :subscript-subscript
 	   :sector
 	   :make-sector
 	   :sector-numbers
 	   :sector-array
 	   :list-+
 	   :list--
 	   :list-zerop
 	   :functional-tensor
 	   :conjugate-index
 	   :tensor-rank
 	   :tensor-dimensions
 	   :tensor-total-size
 	   :tensor-scalar
 	   :tensor-normalize
 	   :tensor-conjugate
 	   :tensor-contract
 	   :tensor-permute
 	   :tensor-fuse
 	   :tensor-decompose
 	   :tensor-exponential))
 
 (in-package :tensor)
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 (defstruct tensor
   (indices nil :type list)
   (sectors nil :type list))
 
 (defstruct segment
   (numbers nil :type list)
   (dimension 0 :type fixnum))
 
 (defstruct subscript
   (numbers nil :type list)
   (subscript 0 :type fixnum))
 
 (defstruct sector
   (numbers nil :type list)
   array)
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 (defun list-+ (&rest args)
   (apply #'mapcar #'+ args))
 
 ;; (defun list-- (&rest args)
 ;;   (apply #'mapcar #'- args))
 
 ;; optimized version with a single argument
 
 (defun list-- (l)
   (loop for i in l collect (- i)))
 
 (defun list-zerop (list)
   (if list
       (and (zerop (the fixnum (car list))) (list-zerop (cdr list)))
       t))
 
 (defun list-< (a b)
   (cond
     ((null a) (not (null b)))
     ((null b) nil)
     ((= (the fixnum (car a)) (the fixnum (car b))) (list-< (cdr a) (cdr b)))
     (t (< (the fixnum (car a)) (the fixnum (car b))))))
 
 (defun array-zerop (array)
   (declare (type (array blas-float) array))
   (dotimes (i (array-total-size array) t)
     (unless (zerop (row-major-aref array i)) (return nil))))
 
 (defun sort-index (index)
   (sort (copy-list index) #'list-< :key #'segment-numbers))
 
 (defun remove-zero-sectors (sectors)
   (mapcan
    #'(lambda (s)
        (unless (array-zerop (sector-array s))
 	 (list s)))
    sectors))
 
 (defun tensor-index-sequence (indices)
   (mapcar
    #'(lambda (mi) (mapcar #'nth mi indices))
    (multi-index-sequence (mapcar #'length indices))))
 
 (defun functional-tensor (indices function)
   (let ((si (mapcar #'sort-index indices)))
     (make-tensor
      :indices
      si
      :sectors
      (remove-zero-sectors
       (mapcan
        #'(lambda (ts)
 	   (let ((mn (mapcar #'segment-numbers ts))
 		 (md (mapcar #'segment-dimension ts)))
 	     (when (list-zerop (apply #'list-+ mn)) ; symmetry condition
 	       (let ((a (make-blas-array md)))
 		 (dolist (mi (multi-index-sequence md)
 			  (list (make-sector :numbers mn :array a)))
 		   (setf (apply #'aref a mi)
 			 (coerce (funcall function
 					  (mapcar
 					   #'(lambda (n i)
 					       (make-subscript
 						:numbers n
 						:subscript i))
 					   mn
 					   mi))
 				 'blas-float)))))))
        (tensor-index-sequence si))))))
 
 (defun conjugate-index (index)
   (reverse
    (mapcar #'(lambda (s)
 	       (make-segment
 		:numbers (list-- (segment-numbers s))
 		:dimension (segment-dimension s)))
 	   index)))
 
 (defun tensor-rank (tensor)
   (length (tensor-indices tensor)))
 
 (defun tensor-dimensions (tensor)
   (loop
      for i in (tensor-indices tensor)
      collect (reduce #'+ i :key #'segment-dimension)))
 
 (defun tensor-total-size (tensor)
   (reduce #'+ (tensor-sectors tensor) :key #'(lambda (s) (array-total-size (sector-array s)))))
 
 (defun tensor-scalar (tensor)
   (let ((s (tensor-sectors tensor)))
     (if s
 	(aref (sector-array (first s)))
 	+blas-0+)))
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 ;; tensor-normalize
 
 (defun tensor-norm (tensor)
   (sqrt
    (loop
       for s in (tensor-sectors tensor)
       sum (expt (array-norm (sector-array s)) 2))))
 
 (defun tensor-normalize (tensor)
   (let ((n (tensor-norm tensor)))
     (make-tensor
      :indices
      (tensor-indices tensor)
      :sectors
      (mapcar
       #'(lambda (s)
 	  (make-sector
 	   :numbers (sector-numbers s)
 	   :array (array-scalar-/ (sector-array s) n)))
       (tensor-sectors tensor)))))
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 ;; tensor-conjugate
 
 (defun tensor-conjugate (tensor)
   (make-tensor
    :indices
    (mapcar
     #'(lambda (i)
 	(reverse
 	 (mapcar #'(lambda (s)
 		     (make-segment
 		      :numbers (list-- (segment-numbers s))
 		      :dimension (segment-dimension s)))
 		 i)))
     (tensor-indices tensor))
    :sectors
    (mapcar
     #'(lambda (s)
 	(make-sector
 	 :numbers (mapcar #'list-- (sector-numbers s))
 	 :array (array-conjugate (sector-array s))))
     (tensor-sectors tensor))))
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 ;; tensor-permute
 
 (defun tensor-permute (tensor permutation)
   (make-tensor
    :indices
    (list-permute (tensor-indices tensor) permutation)
    :sectors
    (mapcar
     #'(lambda (s)
 	(make-sector
 	 :numbers (list-permute (sector-numbers s) permutation)
 	 :array (array-permute (sector-array s) permutation)))
     (tensor-sectors tensor))))
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 ;; tensor-contract
 
 (defun index-match-p (x y)
   (and (= (length x) (length y))
        (loop
 	  for i in x
 	  for j in (reverse y)
 	  unless (and (list-zerop (list-+ (segment-numbers i)
 					  (segment-numbers j)))
 		      (zerop (- (segment-dimension i)
 				(segment-dimension j))))
 	  return nil
 	  finally (return t))))
 
 (defun indices-match-p (x y)
   (and (= (length x) (length y))
        (loop
 	  for i in x
 	  for j in y
 	  unless (index-match-p i j)
 	  return nil
 	  finally (return t))))
 
 (defun contraction-hash-table (tensor indices)
   (let ((h (make-hash-table :test 'equal)))
     (dolist (x (tensor-sectors tensor) h)
       (push x
 	    (gethash (loop
 			with n = (sector-numbers x)
 			for i in indices
 			collect (nth i n))
 		     h)))))
 
 (defun sector-contract (sector1 indices1 sector2 indices2)
   (make-sector
    :numbers
    (nconc
     (loop
        with n = (sector-numbers sector1)
        for i below (length n)
        unless (member i indices1)
        collect (nth i n))
     (loop with n = (sector-numbers sector2)
        for i below (length n)
        unless (member i indices2)
        collect (nth i n)))
    :array
    (array-contract (sector-array sector1)
 		   indices1
 		   (sector-array sector2)
 		   indices2)))
 
 (defun tensor-contract (tensor1 indices1 tensor2 indices2)
   (if (numberp indices1)
       (setf indices1 (list indices1)))
   (if (numberp indices2)
       (setf indices2 (list indices2)))
   (unless (indices-match-p (loop
 			      with n = (tensor-indices tensor1)
 			      for i in indices1
 			      collect (nth i n))
 			   (loop
 			      with n = (tensor-indices tensor2)
 			      for i in indices2
 			      collect (nth i n)))
     (error "Contraction indices do not match."))
   (let ((h1 (contraction-hash-table tensor1 indices1))
 	(h2 (contraction-hash-table tensor2 indices2))
 	(h (make-hash-table :test 'equal)))
     (loop
        for k being the hash-keys in h1
        using (hash-value v1)
        do (let ((v2 (gethash (loop for i in k collect (list-- i)) h2)))
 	    (when v2
 	      (dolist (s1 v1)
 		(dolist (s2 v2)
 		  (let* ((s (sector-contract s1 indices1 s2 indices2))
 			 (n (sector-numbers s))
 			 (v (gethash n h)))
 		    (if v
 			(array-addf (sector-array v)
 				    (sector-array s))
 			(setf (gethash n h) s))))))))
     (make-tensor
      :indices
      (nconc
       (loop
 	 with n = (tensor-indices tensor1)
 	 for i below (length n)
 	 unless (member i indices1)
 	 collect (nth i n))
       (loop
 	 with n = (tensor-indices tensor2)
 	 for i below (length n)
 	 unless (member i indices2)
 	 collect (nth i n)))
      :sectors
      (remove-zero-sectors
       (loop
 	 for v being the hash-values of h
 	 collect v)))))
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 ;; tensor-fuse
 
 (defstruct joint
   segments
   subscript)
 
 (defun joint-numbers (joint)
   (mapcar #'segment-numbers (joint-segments joint)))
 
 (defun joint-number (joint)
   (apply #'list-+ (joint-numbers joint)))
 
 (defun joint-dimensions (joint)
   (mapcar #'segment-dimension (joint-segments joint)))
 
 (defun joint-dimension (joint)
   (apply #'* (joint-dimensions joint)))
 
 (defun fused-index-hash-table (indices type)
   (let ((md (mapcar #'length indices))
 	(hi (make-hash-table :test #'equal)))
     (dolist (mi (multi-index-sequence md))
       (let ((f (make-joint :segments (mapcar #'nth mi indices))))
 	(push f (gethash (joint-number f) hi))))
     (loop
        for k being the hash-keys in hi
        using (hash-value v)
        do (let ((l (sort v #'list-< :key #'(lambda (n) (apply #'append (joint-numbers n)))))
 		(hf (make-hash-table :test #'equal)))
 	    (loop
 	       with s = 0
 	       for i in (if (eq type :reverse) (reverse l) l)
 	       do
 		 (setf (joint-subscript i) s)
 		 (setf (gethash (joint-numbers i) hf) i)
 		 (incf s (joint-dimension i)))
 	    (setf (gethash k hi) hf))
        finally (return hi))))
 
 (defun fused-index-joints (number fiht)
   (loop
      for v being the hash-values in (gethash number fiht)
      collect v))
 
 (defun fused-index-dimension (numbers fiht)
   (loop
      for v being the hash-values in (gethash numbers fiht)
      sum (joint-dimension v)))
 
 (defun fused-index-subscript (numbers fiht)
   (joint-subscript (gethash numbers (gethash (apply #'list-+ numbers) fiht))))
 
 (defun fused-index (fiht)
   (sort-index
    (loop
       for k being the hash-keys in fiht
       using (hash-value v)
       collect
 	(make-segment
 	 :numbers
 	 k
 	 :dimension
 	 (loop
 	    for w being the hash-values in v
 	    sum (joint-dimension w))))))
 
 (defun partite-list (list borders)
   (loop
      with m
      with s
      with b = (pop borders)
      for i below (length list)
      for j in list
      when (and b (= i b))
      do
        (push (nreverse s) m)
        (setf s nil
 	     b (pop borders))
      do
        (push j s)
      finally
        (push (nreverse s) m)
        (return (nreverse m))))
 
 (defun array-fuse (array borders)
   (make-blas-array
    (mapcar #'(lambda (m) (apply #'* m))
 	   (partite-list (array-dimensions array) borders))
    :displaced-to array))
 
 ;; (defun array-insert-subarray (array subscripts subarray)
 ;;   (dolist (mi (multi-index-sequence (array-dimensions subarray))
 ;; 	   subarray)
 ;;     (setf (apply #'aref array (mapcar #'+ subscripts mi))
 ;; 	  (apply #'aref subarray mi))))
 
 ;; optimized version
 
 (defun make-array-insert-subarray-fn (rank)
   (let ((mi (loop for i below rank collect (gensym)))
 	(md (loop for i below rank collect (gensym)))
 	(ms (loop for i below rank collect (gensym))))
     (labels ((nest (&optional (i 0))
     	       (if (eql i rank)
     		   `(setf (aref array
 				,@(loop for j below rank
 				     collect `(+ ,(nth j ms) ,(nth j mi))))
 			  (aref subarray ,@mi))
 		   `(dotimes (,(nth i mi) ,(nth i md)
 			      ,@(when (eql i 0) '(subarray)))
 		      ,(nest (1+ i))))))
       `(lambda (array subscripts subarray)
 	 (declare (type (array blas-float) array subarray))
 	 (let (,@(loop for i below rank
 		    collect `(,(nth i md) (array-dimension subarray ,i)))
 	       ,@(loop for i below rank
 		    collect `(,(nth i ms) (nth ,i subscripts))))
 	   (declare (type fixnum ,@ms))
 	   ,(nest))))))
 
 (let ((cache (make-hash-table)))
   (defun array-insert-subarray (array subscripts subarray)
     (let ((l (length subscripts)))
       (funcall (or (gethash l cache)
 		   (setf (gethash l cache)
 			 (symbol-function
 			  (compile (intern (format nil "ARRAY-INSERT-SUBARRAY-~D" l))
 				   (make-array-insert-subarray-fn l)))))
 	       array
 	       subscripts
 	       subarray))))
 
 ;; (defun array-extract-subarray (array subscripts dimensions)
 ;;   (let ((a (make-blas-array dimensions)))
 ;;     (dolist (mi (multi-index-sequence dimensions) a)
 ;;       (setf (apply #'aref a mi)
 ;; 	    (apply #'aref array (mapcar #'+ subscripts mi))))))
 
 ;; optimized version
 
 (defun make-array-extract-subarray-fn (rank)
   (let ((mi (loop for i below rank collect (gensym)))
 	(md (loop for i below rank collect (gensym)))
 	(ms (loop for i below rank collect (gensym))))
     (labels ((nest (&optional (i 0))
     	       (if (eql i rank)
     		   `(setf (aref a ,@mi)
 			  (aref array
 				,@(loop for j below rank
 				     collect `(+ ,(nth j ms) ,(nth j mi)))))
 		   `(dotimes (,(nth i mi) ,(nth i md)
 			      ,@(when (eql i 0) '(a)))
 		      ,(nest (1+ i))))))
       `(lambda (array subscripts dimensions)
 	 (declare (type (array blas-float) array)
 		  (type list dimensions))
 	 (let ((a (make-blas-array dimensions))
 	       ,@(loop for i below rank
 		    collect `(,(nth i md) (nth ,i dimensions)))
 	       ,@(loop for i below rank
 		    collect `(,(nth i ms) (nth ,i subscripts))))
 	   (declare (type fixnum ,@md ,@ms))
 	   ,(nest))))))
 
 (let ((cache (make-hash-table)))
   (defun array-extract-subarray (array subscripts dimensions)
     (let ((l (length subscripts)))
       (funcall (or (gethash l cache)
 		   (setf (gethash l cache)
 			 (symbol-function
 			  (compile (intern (format nil "ARRAY-EXTRACT-SUBARRAY-~D" l))
 				   (make-array-extract-subarray-fn l)))))
 	       array
 	       subscripts
 	       dimensions))))
 
 (defun make-fused-array (numbers fiht)
   (make-blas-array
    (mapcar #'(lambda (n h) (fused-index-dimension n h)) numbers fiht)
    :initial-element +blas-0+))
 
 (defun fused-array-subscripts (partition fiht)
   (mapcar #'(lambda (n h) (fused-index-subscript n h)) partition fiht))
 
 (defun tensor-fuse (tensor borders types)
   (let ((f (mapcar #'fused-index-hash-table
 		   (partite-list (tensor-indices tensor) borders)
 		   types))
 	(h (make-hash-table :test 'equal)))
     (dolist (s (tensor-sectors tensor))
       (let* ((p (partite-list (sector-numbers s) borders))
 	     (n (mapcar #'(lambda (l) (apply #'list-+ l)) p))
 	     (v (gethash n h)))
 	;; create the array if it does not exist
 	(unless v
 	  (setf v (make-fused-array n f))
 	  (setf (gethash n h) v))
 	;; insert sector array into the fused array
 	(array-insert-subarray v
 			       (fused-array-subscripts p f)
 			       (array-fuse (sector-array s) borders))))
     (make-tensor
      :indices
      (mapcar #'fused-index f)
      :sectors
      (remove-zero-sectors
       (loop
 	 for k being the hash-keys in h
 	 using (hash-value v)
 	 collect
 	   (make-sector :numbers k :array v))))))
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 ;; tensor-decompose
 
 (defstruct svd
   numbers
   u
   s
   v)
 
 (defun svd-reduce-dimensions (svd maxdim normalization)
   (let* ((u (svd-u svd))
 	 (s (svd-s svd))
 	 (v (svd-v svd))
 	 (ud (array-dimension u 0))
 	 (sd (min maxdim (length s)))
 	 (vd (array-dimension v 1))
 	 (un (make-blas-array (list ud sd)))
 	 (vn (make-blas-array (list sd vd)))
 	 (sn (map 'vector
 		  (if (member normalization '(:none :left :right))
 		      #'identity
 		      #'sqrt)
 		  (subseq s 0 sd))))
     (case normalization
       (:none
        (dotimes (i ud)
 	 (dotimes (j sd)
 	   (setf (aref un i j) (aref u i j))))
        (dotimes (i sd)
 	 (dotimes (j vd)
 	   (setf (aref vn i j) (aref v i j)))))
       (:left
        (dotimes (i ud)
 	 (dotimes (j sd)
 	   (setf (aref un i j) (aref u i j))))
        (dotimes (i sd)
 	 (dotimes (j vd)
 	   (setf (aref vn i j) (* (aref v i j) (aref sn i))))))
       (:right
        (dotimes (i ud)
 	 (dotimes (j sd)
 	   (setf (aref un i j) (* (aref u i j) (aref sn j)))))
        (dotimes (i sd)
 	 (dotimes (j vd)
 	   (setf (aref vn i j) (aref v i j)))))
       (:mixed
        (dotimes (i ud)
 	 (dotimes (j sd)
 	   (setf (aref un i j) (* (aref u i j) (aref sn j)))))
        (dotimes (i sd)
 	 (dotimes (j vd)
 	   (setf (aref vn i j) (* (aref v i j) (aref sn i)))))))
     (make-svd :numbers (svd-numbers svd) :u un :v vn)))
 
 (defun diagonal-matrix (diagonal maxdim)
   (let* ((n (min maxdim (length diagonal)))
 	 (m (make-blas-array (list n n) :initial-element +blas-0+)))
     (dotimes (i n m)
       (setf (aref m i i) (coerce (aref diagonal i) 'blas-float)))))
 
 (defun array-reshape (array dimensions)
   (let ((a (make-blas-array dimensions)))
     (dotimes (i (array-total-size a) a)
       (setf (row-major-aref a i) (row-major-aref array i)))))
 
 (defun tensor-decompose (tensor indices maxdim normalization)
   (if (numberp indices)
       (setf indices (list indices)))
   ;; permute if necessary
   (let ((r (tensor-rank tensor)))
     (unless (and (eql (indices-position r indices) :left)
 		 (apply #'< indices))
       (let ((l (loop for i below (length indices) collect i)))
 	(setf tensor (tensor-permute tensor (make-indices-permutation r indices l)))
 	(setf indices l))))
   (let* ((b (list (length indices)))
 	 (p (partite-list (tensor-indices tensor) b))
 	 (f (mapcar #'fused-index-hash-table p '(:normal :reverse)))
 	 (ah (make-hash-table :test #'equal))
 	 (sh (make-hash-table :test #'eq))
 	 sl i1 i2 is1 is2 s1 s2 ss)
     ;; create fused arrays
     (loop
        for s in (tensor-sectors tensor)
        do
 	 (let* ((r (partite-list (sector-numbers s) b))
 		(n (mapcar #'(lambda (l) (apply #'list-+ l)) r))
 		(a (gethash n ah)))
 	   (unless a
 	     (setf a (make-fused-array n f))
 	     (setf (gethash n ah) a))
 	   (array-insert-subarray
 	    a
 	    (fused-array-subscripts r f)
 	    (array-fuse (sector-array s) b))))
     ;; decompose arrays
     (loop
        for n being the hash-keys in ah
        using (hash-value a)
        do
 	 (multiple-value-bind (u s v) (array-decompose a 0)
 	   (let ((x (make-svd :numbers n :u u :s s :v v)))
 	     (loop
 		for i across s
 		unless (essentially-zerop i)
 		do (push (cons i x) sl)))))
     (loop
        with r = (sort sl #'> :key #'car)
        with l = (length r)
        with e = (car maxdim)
        with m = (cadr maxdim)
        with n = (caddr maxdim)
        with u = (and e (* (- 1 e)
 			  (loop
 			     for (s . x) in r
 			     sum (* s s))))
        with c = (if u
 		    (loop
 		       for (s . x) in r
 		       count s
 		       sum (* s s) into v
 		       until (> v u))
 		    l)
        with k = (cond ((and m (< c m) (/= c l))
 		       ;; (format *error-output* "hit mindim (~D)~%" c)
 		       (min l m))
 		      ((and n (> c n))
 		       ;; (format *error-output* "hit maxdim (~D/~D)~%" c l)
 		       n)
 		      (t c))
        repeat k
        for (s . x) in r
        do (if (gethash x sh)
 	      (incf (gethash x sh))
 	      (setf (gethash x sh) 1)))
     ;; reduce the array dimensions
     (loop
        for s being the hash-keys in sh
        using (hash-value d)
        do
 	 (let* ((r (svd-reduce-dimensions s d normalization))
 		(n (svd-numbers r))
 		(n1 (first n))
 		(n2 (second n))
 		(l1 (fused-index-joints n1 (first f)))
 		(l2 (fused-index-joints n2 (second f))))
 	   (push (make-segment :numbers n2 :dimension d) i1)
 	   (push (make-segment :numbers n1 :dimension d) i2)
 	   (push (make-segment :numbers n1 :dimension d) is1)
 	   (push (make-segment :numbers n2 :dimension d) is2)
 	   ;; first sectors
 	   (dotimes (i (length l1))
 	     (let ((j (nth i l1)))
 	       (push (make-sector
 		      :numbers
 		      (append (joint-numbers j) (list n2))
 		      :array
 		      (array-reshape
 		       (array-extract-subarray
 			(svd-u r)
 			(list (joint-subscript j) 0)
 			(list (joint-dimension j) d))
 		       (append (joint-dimensions j) (list d))))
 		     s1)))
 	   ;; second sectors
 	   (dotimes (i (length l2))
 	     (let ((j (nth i l2)))
 	       (push (make-sector
 		      :numbers
 		      (append (list n1) (joint-numbers j))
 		      :array
 		      (array-reshape
 		       (array-extract-subarray
 			(svd-v r)
 			(list 0 (joint-subscript j))
 			(list d (joint-dimension j)))
 		       (append (list d) (joint-dimensions j))))
 		     s2)))
 	   ;; singular values
 	   (push (make-sector
 		  :numbers (list n1 n2)
 		  :array (diagonal-matrix (svd-s s) d))
 		 ss)
 	   ))
     ;; return decomposition
     (values
      (make-tensor :indices (append (first p) (list (sort-index i1)))
 		  :sectors s1)
      (make-tensor :indices (list (sort-index is1) (sort-index is2))
 		  :sectors ss)
      (make-tensor :indices (append (list (sort-index i2)) (second p))
 		  :sectors s2))))
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 ;; tensor-exponential
 
 (defun tensor-permuted-exponential (tensor borders parameter)
   (let* ((f (mapcar #'fused-index-hash-table
 		    (partite-list (tensor-indices tensor) borders)
 		    '(:normal :reverse)))
 	 (h (make-hash-table :test 'equal)))
     ;; add sectors; zero sectors will turn into unit sectors in the exponential
     (loop
        for k being the hash-keys in (first f)
        do (let ((n (list k (list-- k))))
 	    (setf (gethash n h) (make-fused-array n f))))
     (dolist (s (tensor-sectors tensor))
       (let* ((p (partite-list (sector-numbers s) borders))
 	     (n (mapcar #'(lambda (l) (apply #'list-+ l)) p))
 	     (v (gethash n h)))
 	;; insert sector array into the fused array
 	(array-insert-subarray v
 			       (fused-array-subscripts p f)
 			       (array-fuse (sector-array s) borders))))
     (make-tensor
      :indices
      (tensor-indices tensor)
      :sectors
      ;; todo: are zero sectors possible here?
      (remove-zero-sectors
       (loop
 	 with l
 	 for k being the hash-keys in h
 	 using (hash-value v)
 	 do
 	   (setf v (hermitian-matrix-exponential v parameter))
 	   (let ((ml (mapcar #'fused-index-joints k f)))
 	     (dolist (mi (multi-index-sequence (mapcar #'length ml)))
 	       (let* ((j (mapcar #'nth mi ml)))
 		 (push (make-sector
 			:numbers
 			(mapcan #'joint-numbers j)
 			:array
 			(array-reshape
 			 (array-extract-subarray
 			  v
 			  (mapcar #'joint-subscript j)
 			  (mapcar #'joint-dimension j))
 			 (mapcan #'joint-dimensions j)))
 		       l))))
 	 finally (return l))))))
 
 (defun tensor-exponential (tensor groups parameter)
   (unless (= (length groups) 2)
     (error "Incorrect numbers of index groups."))
   (let* ((r (tensor-rank tensor))
 	 (o (apply #'append groups))
 	 (p (loop for i below r collect i)))
     (tensor-permute
      (tensor-permuted-exponential
       (tensor-permute tensor (make-indices-permutation r o p))
       (list (length (first groups)))
       parameter)
      (make-indices-permutation r p o))))