;;; We're going to try to be *very* general here, and define things in
;;; terms of Vector and Operator.  Note: since we distinguish between
;;; vector and Vector, we are relying on Gambit-C being case sensitive.

;;; OK, I've now come to the conclusion that there is no advantage
;;; in a CLOS-like system like Meroon with separate class and generic
;;; hierarchies to defining a class without
;;; slots, unless for some error-checking value in generics.
;;; Vector doesn't pass this test, since the generics have to work
;;; on f64vectors, too.  So the Vector class is history.
;;; The only things that subclassed vectors were Operators and
;;; Finite-element-vectors; another note that they weren't really
;;; needed.

;;;(define-class Vector Object ()) ; vectors over reals

(define-class Operator Object 
  ((= domain? maybe-uninitialized:)      ; a predicate, but how do we test whether two domains are the same?
   (= range? maybe-uninitialized:)       ; a predicate, but how do we test whether two ranges are the same?
   (= mapping maybe-uninitialized:)))   ; #f if it cannot be applied easily

;;; Right now, Meroon cannot dispatch on built-in types, so in the default methods
;;; we check for f64vectors of the same length and define operations for those.
;;; Define some generic functions whose default behavior is defined only
;;; for f64vectors.

(define-generic (Vector-add (x) y)
  (if (and (f64vector? x)
	   (f64vector? y)
	   (fx= (f64vector-length x)
		(f64vector-length y)))
      (let ((n (f64vector-length x)))
	(let ((result (make-f64vector n)))
	  (do ((i (fx- n 1) (fx- i 1)))
	      ((fx< i 0) result)
	    (f64vector-set! result i (fl+ (f64vector-ref x i)
					  (f64vector-ref y i))))))
      (error "Vector-add is not defined for these objects.")))

(define-generic (Vector-subtract (x) y)
  (if (and (f64vector? x)
	   (f64vector? y)
	   (fx= (f64vector-length x)
		(f64vector-length y)))
      (let ((n (f64vector-length x)))
	(let ((result (make-f64vector n)))
	  (do ((i (fx- n 1) (fx- i 1)))
	      ((fx< i 0) result)
	    (f64vector-set! result i (fl- (f64vector-ref x i)
					  (f64vector-ref y i))))))
      (error "Vector-subtract is not defined for these objects.")))

(define-generic (Vector-scale x (y))  ; scalar multiplication
  (if (f64vector? y)
      (let ((n (f64vector-length y)))
	(let ((result (make-f64vector n)))
	  (do ((i (fx- n 1) (fx- i 1)))
	      ((fx< i 0) result)
	    (f64vector-set! result i (fl* x (f64vector-ref y i))))))
      (error "Vector-scale is not defined for this object.")))

(define-generic (Vector-*axpy a (x) y)
  (if (and (f64vector? x)
	   (f64vector? y)
	   (fx= (f64vector-length x)
		(f64vector-length y)))
      (let ((n (f64vector-length x)))
	(let ((result (make-f64vector n)))
	  (do ((i (fx- n 1) (fx- i 1)))
	      ((fx< i 0) result)
	    (f64vector-set! result i (fl+ (fl* a (f64vector-ref x i))
					  (f64vector-ref y i))))))
      (error "Vector-*axpy is not defined for these objects.")))

(define-generic (Vector-copy (x))
  (if (f64vector? x)
      (let ((n (f64vector-length x)))
	(let ((result (make-f64vector n)))
	  (do ((i (fx- n 1) (fx- i 1)))
	      ((fx< i 0) result)
	    (f64vector-set! result i (f64vector-ref x i)))))
      (error "Vector-copy is not defined for this object.")))

;;; max and min are defined for lattices.  We'll define it
;;; here for f64vectors, since it is  useful; there is no
;;; Lattice class for the same reason there is no longer a
;;; Vector class.

(define-generic (Lattice-max (x) y)
  (if (and (f64vector? x)
	   (f64vector? y)
	   (fx= (f64vector-length x)
		(f64vector-length y)))
      (let ((n (f64vector-length x)))
	(let ((result (make-f64vector n)))
	  (do ((i (fx- n 1) (fx- i 1)))
	      ((fx< i 0) result)
	    (f64vector-set! result i (flmax (f64vector-ref x i)
					    (f64vector-ref y i))))))
      (error "Lattice-max is not defined for these objects.")))

(define-generic (Lattice-min (x) y)
  (if (and (f64vector? x)
	   (f64vector? y)
	   (fx= (f64vector-length x)
		(f64vector-length y)))
      (let ((n (f64vector-length x)))
	(let ((result (make-f64vector n)))
	  (do ((i (fx- n 1) (fx- i 1)))
	      ((fx< i 0) result)
	    (f64vector-set! result i (flmin (f64vector-ref x i)
					    (f64vector-ref y i))))))
      (error "Lattice-min is not defined for these objects.")))

(define-generic (Vector-dimension (x))
  (if (f64vector? x)   ; we'll define it for Rn
      (f64vector-length x)
      (error "Vector-dimension is not defined for this object.")))

(define-generic (Vector-zero (x))
  (if (f64vector? x)
      (make-f64vector (f64vector-length x) 0.0)
      (error "Vector-zero is not defined for this object.")))

(define-generic (Vector-dot-product (x) y)
  (declare (inlining-limit 0))  ;; no more unrolling
  (if (and (f64vector? x)
	   (f64vector? y)
	   (fx= (f64vector-length x)
		(f64vector-length y)))
      (let outer ((i (fx- (f64vector-length x) 1))
		  (sum 0.0))
	(if (fx<= 7 i)
	    (outer (fx- i 8)
		   (fl+ sum
			(fl* (f64vector-ref x i) 
			     (f64vector-ref y i))
			(fl* (f64vector-ref x (fx- i 1))
			     (f64vector-ref y (fx- i 1)))
			(fl* (f64vector-ref x (fx- i 2)) 
			     (f64vector-ref y (fx- i 2)))
			(fl* (f64vector-ref x (fx- i 3)) 
			     (f64vector-ref y (fx- i 3)))
			(fl* (f64vector-ref x (fx- i 4)) 
			     (f64vector-ref y (fx- i 4)))
			(fl* (f64vector-ref x (fx- i 5)) 
			     (f64vector-ref y (fx- i 5)))
			(fl* (f64vector-ref x (fx- i 6)) 
			     (f64vector-ref y (fx- i 6)))
			(fl* (f64vector-ref x (fx- i 7)) 
			     (f64vector-ref y (fx- i 7)))))
	    (let inner ((i i)
			(sum sum))
	      (if (fx<= 0 i)
		  (inner (fx- i 1)
			 (fl+ (fl* (f64vector-ref x i)
				   (f64vector-ref y i))
			      sum))
		  sum))))
      (error "Vector-dot-product is not defined for these objects.")))

;;; This next naive implementation may seem memory inefficient, but
;;; it actually takes about half the memory of the "optimized" loop
;;; for Linear-finite-element-vectors because Gambit does not keep
;;; flonums unboxed across jumps.

(define-generic (Vector-inner-product (x) (operator Operator) y)
  (Vector-dot-product x (Operator-apply operator y)))

(define-generic (Vector-L2-norm (x))
  (flsqrt (Vector-dot-product x x)))

(define-generic (Vector-Linf-norm (x))
  (if (f64vector? x)
      (reduce-from-left
       (lambda (x y) (flmax x (flabs y)))
       -inf.0
       (f64vector->list x))
      (error "Vector-Linf-norm: argument is not a Vector")))

(define-generic (Vector-L1-norm (x))
  (if (f64vector? x)
      (reduce-from-left
       (lambda (x y) (fl+ x (flabs y)))
       0.0
       (f64vector->list x))
      (error "Vector-L1-norm: argument is not a Vector")))

;;; We'll define a blocked vector, that simply has a vector of blocks,
;;; each of which is a Vector (which is purely conceptual).  Note: The
;;; individual blocks do *not* need to be f64vectors, they can be any
;;; type of Vectors.

(define-class Blocked-vector (Object)
  ((= blocks)))  ; a vector of blocks

(define-handy-method (Vector-add (x Blocked-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Blocked-vector)
			blocks: (vector-map Vector-add blocks (Blocked-vector-blocks y)))
      (error "Vector-add: Arguments are not the same type" x y)))

(define-handy-method (Vector-subtract (x Blocked-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Blocked-vector)
			blocks: (vector-map Vector-subtract blocks (Blocked-vector-blocks y)))
      (error "Vector-subtract: Arguments are not the same type " x y)))

(define-handy-method (Vector-scale a (x Blocked-vector))
  (instantiate-from (x Blocked-vector)
		    blocks: (vector-map (lambda (x) (Vector-scale a x)) blocks)))

(define-handy-method (Vector-*axpy a (x Blocked-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Blocked-vector)
			blocks: (vector-map (lambda (x y) (Vector-*axpy a x y))
					      blocks
					      (Blocked-vector-blocks y)))
      (error "Vector-*axpy: Arguments are not of the same type " x y)))

(define-handy-method (Vector-copy (x Blocked-vector))
  (instantiate-from (x Blocked-vector)
		    blocks: (vector-map Vector-copy blocks)))

(define-handy-method (Vector-zero (x Blocked-vector))
  (instantiate-from (x Blocked-vector)
		    blocks: (vector-map Vector-zero blocks)))

;;; Vector-dimension is not supplied for Blocked-vector because of possibly redundant
;;; representations.

;;; Similarly for Vector-dot-product

(define-handy-method (Lattice-max (x Blocked-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Blocked-vector)
			blocks: (vector-map Lattice-max blocks (Blocked-vector-blocks y)))
      (error "Lattice-max: Arguments are not the same type " x y)))

(define-handy-method (Lattice-min (x Blocked-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Blocked-vector)
			blocks: (vector-map Lattice-min blocks (Blocked-vector-blocks y)))
      (error "Lattice-min: Arguments are not the same type " x y)))

;;; applies op to v
;;; After using this code for 18 months, I discovered that
;;; Operator-apply had no methods defined for it.  This implies
;;; to me that this should not be a generic function, but a plain
;;; function.  It also implies that an Operator should be an applyable
;;; object directly without going through Operator-apply, which
;;; means that it should be a subclass of Generic, but that hurts
;;; my brain too damn much to think about.

(define (Operator-apply op v)
  (with-access
   op (Operator mapping domain?)
   (if (and (field-defined? op 'mapping)
	    mapping)
       (if (domain? v)
	   (mapping v)
	   (error "Operator-apply---vector is not in domain of operator"))
       (error "Operator-apply---application is not defined for this operator"))))

;; calculates a new operator = op1 \circle op2

(define-generic (Operator-compose (op1 Operator) (op2 Operator))
  ;; here's where we're hurting with our model---we can't test
  ;; that (Operator-range? op2) is (Operator-domain? op1)
  (instantiate Operator
	       domain?: (Operator-domain? op2)
	       range?: (Operator-range? op1)
	       mapping: (lambda (v)
			  (Operator-apply op1 (Operator-apply op2 v)))))

;; Generic identity operator

(define identity-operator
  (instantiate Operator
	       domain?: (lambda (x) #t)
	       range?: (lambda (x) #t)
	       mapping: (lambda (x) x)))

;; this is useless unless subclassed

(define-class Matrix-operator Operator
  ((= domain-dimension)
   (= range-dimension)))

(define-handy-method (initialize! (op Matrix-operator))
  (set! domain? (lambda (v) (and (f64vector? v) 
				 (fx= (f64vector-length v)
				      domain-dimension))))
  (set! range? (lambda (v) (and (f64vector? v) 
				(fx= (f64vector-length v)
				     range-dimension))))
  (call-next-method))

;; Is transpose for Matrix

(define-generic (Operator-adjoint (op Operator))
  (error "Operator-adjoint: method is not defined for this operator " op))

;; Builds an identity operator of the same class as op

(define-generic (Operator->identity (op Operator))
  (error "Operator->identity: method is not defined for this operator " op))