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