[Openmcl-devel] Error when setting a value in a matrix unless safety 3

[Keith James]

Hi,

I'm porting some of my libraries that were developed on SBCL and
LispWorks. I'm experiencing an error that I can't explain when setting
a value (a single-float) into a type-declared matrix of single
floats. I've discovered that I can avoid the error by setting safety 3.

The error message is

value -5.0 is not of the expected type (SIMPLE-ARRAY DOUBLE-FLOAT (* *))

As the algorithm is very compute intensive (Smith-Waterman alignment),
I would like to optimize for speed and slightly reduced safety. I have
extracted the smallest self-contained example that I can, illustrating
the problem.

Most of the code has been macroexpanded, hence the odd variable
names. I've commented the exact point of the error. Any suggestions
would be much appreciated.

I'm using Version 1.3-r11936 (LinuxX8664). 

Thanks!


;; Run (test-alignment "AACGTTTAATGTAC" "CGTTTGTAA" #'simple-dna-subst
;;                     -5.0f0 -1.0f0)

(deftype array-index ()
  "Array index type."
  '(and fixnum (integer 0 *)))

(deftype path-pointer ()
  "Dynamic programming backtrace pointer."
  '(unsigned-byte 2))

(defconstant +insertx+ 1
  "Dynamic programming backtrace delete pointer.")
(defconstant +inserty+ 2
  "Dynamic programming backtrace insert pointer.")
(defconstant +match+ 3
  "Dynamic programming backtrace match pointer.")
(defconstant +none+ 0
  "Dynamic programming backtrace null pointer.")

(defvar *simple-dna-matrix*
  (make-array '(5 5)
              :element-type 'single-float
              :initial-contents '(( 5.0f0 -4.0f0 -4.0f0 -4.0f0 -1.0f0)
                                  (-4.0f0  5.0f0 -4.0f0 -4.0f0 -1.0f0)
                                  (-4.0f0 -4.0f0  5.0f0 -4.0f0 -1.0f0)
                                  (-4.0f0 -4.0f0 -4.0f0  5.0f0 -1.0f0)
                                  (-1.0f0 -1.0f0 -1.0f0 -1.0f0 -1.0f0)))
  "A DNA substitution matrix for sequences containing the bases A, C,
G, T and N.")

(defun simple-dna-index (value)
         (declare (optimize (speed 3) (safety 1)))
         (declare (type base-char value))
         (cond ((char= value #\A) 0)
               ((char= value #\C) 1)
               ((char= value #\G) 2)
               ((char= value #\T) 3)
               ((char= value #\N) 4)
               (t
                (error "invalid base ~a." value))))

(defun simple-dna-subst (x y)
  (aref *simple-dna-matrix* (simple-dna-index x) (simple-dna-index y)))

(defun test-alignment (vecm vecn subst-fn gap-open gap-extend)
  ;; (declare (optimize (safety 3))) ; Success
  (declare (optimize (safety 2))) ; Fail
  ;; local fn to avoid boxing of returned floats
  (flet ((subn (x y) (funcall subst-fn x y)))
    (let ((m (length vecm))
          (n (length vecn)))
      (let ((sc (make-array (list (1+ m) (1+ n))
                            :element-type 'single-float
                            :initial-element 0.0))
            (ix (make-array (list (1+ m) (1+ n))
                            :element-type 'single-float
                            :initial-element 0.0))
            (iy (make-array (list (1+ m) (1+ n))
                            :element-type 'single-float
                            :initial-element 0.0))
            (bt (make-array (list (1+ m) (1+ n))
                            :element-type 'path-pointer
                            :initial-element 0)))
        (declare (type (simple-array single-float (* *)) sc ix iy)
                 (type (simple-array path-pointer (* *)) bt))
        (let ((var-g0 (array-dimension sc 0))
              (var-g1 (array-dimension sc 1))
              (max-score 0.0)
              (max-row 0)
              (max-col 0))
          (loop
             for row of-type array-index from 1 below var-g0
             for prev-row of-type array-index = (1- row)
             do (loop
                   for col of-type array-index from 1 below var-g1
                   for prev-col of-type array-index = (1- col)
                   do (let ((var-g4
                             (the single-float
                               (subn (aref vecm prev-row)
                                     (aref vecn prev-col)))))
                        (let ((var-g5 (+ (aref sc row prev-col) gap-open))
                              (var-g6 (+ (aref ix row prev-col) gap-extend))
                              (var-g7 (+ (aref iy row prev-col) gap-open))
                              (var-g8 (+ (aref sc prev-row col) gap-open))
                              (var-g9 (+ (aref iy prev-row col) gap-extend))
                              (var-g10 (+ (aref ix prev-row col) gap-open))
                              (var-g11 (+ (aref sc prev-row prev-col) var-g4))
                              (var-g12 (+ (aref ix prev-row prev-col) var-g4))
                              (var-g13 (+ (aref iy prev-row prev-col) var-g4)))
                          (let ((var-g2 (if (= 1 col)
                                            var-g5
                                          (max var-g5 var-g6 var-g7)))
                                (var-g3 (if (= 1 row)
                                            var-g8
                                          (max var-g8 var-g9 var-g10)))
                                (cell-score (max 0.0 var-g11 var-g12 var-g13)))
                            ;; vvv error here vvv
                            (setf (aref ix row col) var-g2
                                  (aref iy row col) var-g3
                                  (aref sc row col) cell-score)
                            ;; ^^^            ^^^
                            (let ((var-g14 (max cell-score var-g2 var-g3)))
                              (setf (aref bt row col)
                                    (cond ((= var-g14 var-g2) +insertx+)
                                          ((= var-g14 var-g3) +inserty+)
                                          (t +match+))))
                            (when (> cell-score max-score)
                              (setf max-score cell-score
                                    max-row row
                                    max-col col)))))))
          (values max-score max-row max-col))))))


-- 

- Keith James - Wellcome Trust Sanger Institute, UK -


-- 
 The Wellcome Trust Sanger Institute is operated by Genome Research 
 Limited, a charity registered in England with number 1021457 and a 
 company registered in England with number 2742969, whose registered 
 office is 215 Euston Road, London, NW1 2BE.