;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; gate-compression.lisp
;;
;; c) 1999-2003, Lee Spector
;; original source of much was matrices.lisp in qc-matrices
;; Oct 5 1999, made compatible with limited-oracle
;; Nov 12 2003, many updates, disentangled from LGP

;; compiler optimization settings

; for debugging
; (eval-when (compile)
;   (declaim (optimize (speed 2) (safety 1) (space 1) (debug 3))))

; for maximum reasonably safe speed
(eval-when (compile)
  (declaim (optimize (speed 3) (safety 1) (space 0) (debug 0))))

(defparameter *max-depth-for-compression* 5)
(defparameter *max-depth-for-history* 5)
(defparameter *uncompressible*
  (list 'oracle 'limited-oracle 'measure 'end 'halt))
(defvar *NUMBER-OF-QUBITS*)
(defvar *ALL-QUBITS*)


(defun matrix-matrix-multiply (matrix1 matrix2)
  (let* ((matrix-size (car (array-dimensions matrix1)))
         (result (make-array (list matrix-size matrix-size))))
    (dotimes (i matrix-size)
      (dotimes (j matrix-size)
        (setf (aref result i j)
              (let ((element 0))
                (dotimes (k matrix-size)
                  (incf element
                        (* (aref matrix1 i k)
                           (aref matrix2 k j))))
                element))))
    result))

(defun expand-matrix (gate targets)
  "Expands the operator matrix gate to a full matrix for operating
on a system of *number-of-qubits* qubits, with the operator being applied
to the qubits specified in targets. Written by Lee Spector, 1999.
Targets reversal added Sept 8, 1999."
  (let* ((targets (reverse targets))
         (m-size (expt 2 *number-of-qubits*))
         (m (make-array (list m-size m-size)))
         (non-targets ;(loop for i from 0 to (- *number-of-qubits* 1)
                      ;      unless (member i targets)
                      ;      collect i)))
	  (let ((collection nil))
	    (do ((i 0 (+ i 1)))
		((>= i *number-of-qubits*) (reverse collection))
		(unless (member i targets) (push i collection))))))
    (dotimes (i m-size)
      (dotimes (j m-size)
        (setf (aref m i j)
              (if (=in-positions non-targets i j)
                (aref gate
                      (extract@positions targets i)
                      (extract@positions targets j))
                0))))
    m))

(defun =in-positions (positions int1 int2)
  "Returns non-nil if int1 and int2 are have the same bits at all
positions in positions. Written by Lee Spector, 1999."
  (every #'(lambda (index)
             (eq (logbitp index int1)
                 (logbitp index int2)))
         positions))

(defun extract@positions (positions int)
  "Returns the number formed by extracting and concatenating 
the bits of int indexed by the positions in positions.
Written by Lee Spector, 1999."
  (let ((exponent -1))
    ;(loop for index in positions
    ;      do (incf exponent)
    ;      when (logbitp index int)
    ;      sum (expt 2 exponent))))
    (let ((sum 0))
      (dolist (index positions)
        (incf exponent)
	(when (logbitp index int)
	  (incf sum (expt 2 exponent))))
      sum)))


(defun long (n) (coerce n 'long-float))

(defun expand-gate-form (gate-form)
  (case (first gate-form)
    (qnot (expand-matrix #2A((0 1) (1 0)) (cdr gate-form)))
    (cnot (expand-matrix #2A((1 0 0 0)
                             (0 1 0 0)
                             (0 0 0 1)
                             (0 0 1 0))
                         (cdr gate-form)))
    (swap (expand-matrix #2A((1 0 0 0)
                             (0 0 1 0)
                             (0 1 0 0)
                             (0 0 0 1))
                         (cdr gate-form)))
    (hadamard (expand-matrix (make-array 
                              '(2 2)
                              :initial-contents 
                              (list (list (/ 1 (sqrt 2.0L0))  (/ 1 (sqrt 2.0L0)))
                                    (list (/ 1 (sqrt 2.0L0))  (- (/ 1 (sqrt 2.0L0))))))
                             (cdr gate-form)))
    (srn (expand-matrix (make-array 
                         '(2 2)
                         :initial-contents 
                         (list (list (/ 1 (sqrt 2.0L0))  (- (/ 1 (sqrt 2.0L0))))
                               (list (/ 1 (sqrt 2.0L0))  (/ 1 (sqrt 2.0L0)))
                               ))
                        (cdr gate-form)))
    (u-theta (expand-matrix (let ((theta (long (third gate-form))))
                              (make-array '(2 2)
                                          :initial-contents 
                                          (list (list (cos theta)  (sin theta))
                                                (list (- (sin theta))  (cos theta)))))
                            (list (second gate-form))))
    (cphase (expand-matrix (let ((alpha (long (fourth gate-form))))
                             (make-array '(4 4)
                                         :initial-contents 
                                         (list (list 1 0 0 0)
                                               (list 0 1 0 0)
                                               (list 0 0 1 0)
                                               (list 0 0 0 (exp (* (sqrt -1) alpha))))))
                           (list (second gate-form) (third gate-form))))
    (u2 (expand-matrix (let ((phi (long (third gate-form)))
                             (theta (long (fourth gate-form)))
                             (psi (long (fifth gate-form)))
                             (alpha (long (sixth gate-form)))
                             (i (sqrt -1)))
                         (make-array
                          '(2 2)
                          :initial-contents 
                          (list (list (* (exp (* i (+ (- phi) (- psi) alpha))) (cos theta))
                                      (* (exp (* i (+ (- phi) psi alpha))) (sin (- theta))))
                                (list (* (exp (* i (+ phi (- psi) alpha))) (sin theta))
                                      (* (exp (* i (+ phi psi alpha))) (cos theta))))))
                       (list (second gate-form))))
    (matrix-gate (second gate-form))))

(defun matrix-gate (qsys matrix history)
  (declare (ignore history))
  (apply-operator qsys
                  matrix
                  (reverse *All-Qubits*)))

(defun max-depth (tree)
  (if (not (listp tree))
    0
    (1+ (apply #'max (mapcar #'max-depth tree)))))
        
(defun process-for-history (gate-sequence)
  "Removes actual matrices from histories containing matrix-gate forms, 
substituting a COMPRESSED form with only the prior history. Punts and 
returns TOO-DEEP if the depth is greater than *max-depth-for-history*."
  (if (> (max-depth gate-sequence) *max-depth-for-history*)
    'too-deep
    (mapcar #'(lambda (gate-form)
                (if (eq (car gate-form) 'matrix-gate)
                  (list 'compressed (third gate-form))
                  gate-form))
            gate-sequence)))

(defun compress-compressible-gate-sequence (seq)
  "assumes all gates can be expanded"
  (cond ((<= (length seq) 1) ;; don't compress a single gate
         seq)
        ((> (max-depth seq) *max-depth-for-compression*)
         seq)
        (t (let ((composite-matrix (expand-matrix #2A((1 0)(0 1)) nil)))
             ;; start with identity
             (dolist (gate-form seq)
               (setq composite-matrix
                     (matrix-matrix-multiply (expand-gate-form gate-form)
                                             composite-matrix)))
             (if (> (check-unitarity composite-matrix) 1.0E-10)
               ;; errors too high
               seq
               ;; errors OK
               (list (list 'matrix-gate
                           composite-matrix
                           (process-for-history 
                            seq))) 
               )))))

(defun thoroughly-compress-compressible-gate-sequence (seq)
  "assumes all gates can be expanded. Just like compress-compressible-gate-sequence
except compresses even single gates into matrix-gates."
  (cond ;((<= (length seq) 1) ;; don't compress a single gate
        ; seq)
        ((> (max-depth seq) *max-depth-for-compression*)
         seq)
        (t (let ((composite-matrix (expand-matrix #2A((1 0)(0 1)) nil)))
             ;; start with identity
             (dolist (gate-form seq)
               (setq composite-matrix
                     (matrix-matrix-multiply (expand-gate-form gate-form)
                                             composite-matrix)))
             (if (> (check-unitarity composite-matrix) 1.0E-10)
               ;; errors too high
               seq
               ;; errors OK
               (list (list 'matrix-gate
                           composite-matrix
                           (process-for-history 
                            seq))) 
               )))))

(defun compress-gates (program &optional (pending nil))
  (cond ((null program)
         (if pending
           (compress-compressible-gate-sequence pending)
           nil))
        ((member (caar program) *uncompressible*)
         (if pending
           (append (compress-compressible-gate-sequence pending)
                   (cons (car program) (compress-gates (cdr program))))
           (cons (car program) (compress-gates (cdr program)))))
        (t (compress-gates (cdr program) 
                           (append pending (list (car program)))))))

(defun thoroughly-compress-gates (program &optional (pending nil))
  "like compress-gates but does it thoroughly"
  (cond ((null program)
         (if pending
           (thoroughly-compress-compressible-gate-sequence pending)
           nil))
        ((member (caar program) *uncompressible*)
         (if pending
           (append (thoroughly-compress-compressible-gate-sequence pending)
                   (cons (car program) (thoroughly-compress-gates (cdr program))))
           (cons (car program) (thoroughly-compress-gates (cdr program)))))
        (t (thoroughly-compress-gates (cdr program) 
                                      (append pending (list (car program)))))))

(defun check-unitarity (m)
  "Returns the cumulative difference of each element of (m times m*) (where
m* is the conjugate transpose of m) from the corresponding element of the
identity matrix. This will be 0 for a unitary matrix."
  (let* ((dim (car (array-dimensions m)))
         (m* (make-array (list dim dim))) ;; conjugate transpose of m
         (identity (make-array (list dim dim)))
         (product nil)
         (cumulative-error 0))
    ;; set up identity matrix
    (dotimes (i dim)
      (dotimes (j dim)
        (setf (aref identity i j) (if (= i j) 1 0))))
    ;; set up m*
    (dotimes (i dim)
      (dotimes (j dim)
        (setf (aref m* i j) 
              (conjugate (aref m j i)))))
    ;; multiply
    (setq product (matrix-matrix-multiply m m*))
    ;(print product)
    ;; sum error
    (dotimes (i dim)
      (dotimes (j dim)
        (incf cumulative-error
              (abs (- (aref product i j)
                      (aref identity i j))))))
    cumulative-error))


;; genetic operator for LGP
#|(defun gate-compression () 
  (let ((parent-program (copy-tree (second (random-individual)))))
    (list :no-computed-fitness 
          (compress-gates parent-program))))|#




#|

(compress-gates '((hadamard 0)))

(compress-gates '((hadamard 0)(hadamard 1)))

(compress-gates (append '((hadamard 0))
                        (compress-gates '((hadamard 0)(hadamard 1))))
                )


;; should get same results with and without compression
(test-quantum-program
 '((hadamard 0) (u-theta 2 1.23) (hadamard 1))  ;; without compression
 ;(compress-gates '((hadamard 0) (u-theta 2 1.23) (hadamard 1))) ;; with compression
 :num-qubits *number-of-qubits*
 :cases *fitness-cases*
 :final-measurement-qubits *measurement-qubits*
 :threshold 0.48
 :inspect nil
 :debug 0)


(thoroughly-compress-gates '((hadamard 0)))

(thoroughly-compress-gates '((hadamard 0)(hadamard 1)))

(thoroughly-compress-gates (append '((hadamard 0))
                                   (thoroughly-compress-gates '((hadamard 0)(hadamard 1))))
                           )


;; should get same results with and without compression
(test-quantum-program
 '((hadamard 0) (u-theta 2 1.23) (hadamard 1))  ;; without compression
 ;(thoroughly-compress-gates '((hadamard 0) (u-theta 2 1.23) (hadamard 1))) ;; with compression
 :num-qubits *number-of-qubits*
 :cases *fitness-cases*
 :final-measurement-qubits *measurement-qubits*
 :threshold 0.48
 :inspect nil
 :debug 0)



|#