compound.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
   3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
   4 ;;; Since 1991, ANSI was finally finished.  Edited for ANSI Common Lisp.
   5
   6 ;;; compound -- Compound data and element-wise mapping functions
   7 ;;;
   8 ;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
   9 ;;; unrestricted use.
  10 ;;;
  11
  12 ;;;
  13 ;;; Package Setup
  14 ;;;
  15
  16 (in-package :cl-user)
  17
  18 (defpackage :lisp-stat-compound-data
  19   (:use :common-lisp
  20         :lisp-stat-object-system
  21         :lisp-stat-types)
  22   (:import-from :lisp-stat-fastmap fastmap)
  23   (:shadowing-import-from :lisp-stat-object-system
  24                           slot-value
  25                           call-next-method call-method)
  26   (:export compound-data-p *compound-data-proto*
  27            compound-object-p
  28            compound-data-seq compound-data-length
  29
  30            element-list element-seq
  31
  32            sort-data order rank
  33
  34            recursive-map-elements map-elements
  35
  36            repeat
  37            ;; export matrix-related functionality (not sure??)
  38
  39            check-sequence
  40            get-next-element make-next-element set-next-element
  41            sequencep iseq
  42
  43            ;; maybe?
  44            ordered-nneg-seq
  45            select
  46
  47            which
  48            ;; vector differences
  49            difference rseq
  50            ))
  51
  52 (in-package :lisp-stat-compound-data)
  53
  54 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  55 ;;;
  56 ;;;                    Internal Support Functions
  57 ;;;
  58 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  59
  60 (defun cmpndp (x)
  61   "Predicate to determine if argument is compound. Most common
  62 non-compound types are checked first."
  63   (declare (inline numberp symbolp stringp consp arrayp array-total-size))
  64   (cond ((or (numberp x) (symbolp x) (stringp x)) nil)
  65          ((or (consp x) (and (arrayp x) (< 0 (array-total-size x)))) t)
  66          (t (compound-object-p x))))
  67
  68 (defun find-compound-data (list)
  69   "Returns first compound data item in LIST or NIL if there is none."
  70   (dolist (x list) (if (cmpndp x) (return x))))
  71
  72 (defun any-compound-elements (seq)
  73   "Checks for a compound element."
  74   (cond ((consp seq) (dolist (x seq) (if (cmpndp x) (return x))))
  75         ((vectorp seq)
  76          (let ((n (length seq)))
  77            (declare (fixnum n))
  78            (dotimes (i n)
  79              (declare (fixnum i))
  80              (let ((x (aref seq i)))
  81                (if (cmpndp x) (return x))))))
  82         (t (error "argument must be a list or vector"))))
  83
  84 (defun compound-data-sequence (x)
  85   "Returns sequence of data values for X."
  86   (declare (inline consp vectorp arrayp make-array array-total-size))
  87   (cond
  88    ((or (consp x) (vectorp x)) x)
  89    ((arrayp x) (make-array (array-total-size x) :displaced-to x))
  90    (t (send x :data-seq))))
  91
  92 (defmacro sequence-type (x) `(if (consp ,x) 'list 'vector))
  93
  94 (defun make-compound-data (shape sequence)
  95 "Construct a compound data item to match the shape of the first
  96 argument."
  97   (let ((n (length (compound-data-sequence shape))))
  98     (if (/= n (length sequence)) (error "compound data not the same shape"))
  99     (cond
 100      ((consp shape) (if (consp sequence) sequence (coerce sequence 'list)))
 101      ((vectorp shape)
 102       (if (vectorp sequence) sequence (coerce sequence 'vector)))
 103      ((arrayp shape)
 104       (make-array (array-dimensions shape)
 105                   :displaced-to (coerce sequence 'vector)))
 106      (t (send shape :make-data sequence)))))
 107
 108 (defun make-circle (x)
 109   "Make a circular list of one element."
 110   (declare (inline cons rplacd))
 111   (let ((x (cons x nil)))
 112     (rplacd x x)
 113     x))
 114
 115 (defun check-compound (x)
 116   "Signals an error if X is not compound."
 117   (if (not (cmpndp x)) (error "not a compound data item - ~a" x)))
 118
 119 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 120 ;;;
 121 ;;; MAP-ELEMENTS function
 122 ;;; Applies a function to arguments. If all arguments are simple (i. e.
 123 ;;; not compound) then MAP-ELEMENTS acts like funcall. Otherwise all
 124 ;;; compound arguments must be of the same shape and simple arguments
 125 ;;; are treated as if they were compound arguments of the appropriate
 126 ;;; shape. This is implemented by replacin all simple arguments by
 127 ;;; circular lists of one element.
 128 ;;;
 129 ;;; This implementation uses FASTMAP, a version of MAP that is assumed
 130 ;;; to
 131 ;;;
 132 ;;;     a) work reasonable fast on any combination of lists and vectors
 133 ;;;        as its arguments
 134 ;;;
 135 ;;;     b) not hang if at least one of its arguments is not a circular
 136 ;;;        list.
 137 ;;;
 138 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 139
 140 (defun fixup-map-elements-arglist (args)
 141   (do* ((args args (rest args))
 142         (x (car args) (car args)))
 143        ((null args))
 144     (declare (inline car))
 145     (setf (car args)
 146           (if (cmpndp x) (compound-data-sequence x) (make-circle x)))))
 147
 148 (defun map-elements (fcn &rest args)
 149 "Args: (fcn &rest args)
 150 Applies FCN elementwise. If no arguments are compound MAP-ELEMENTS
 151 acts like FUNCALL. Compound arguments must all be the same shape. Non
 152 compound arguments, in the presence of compound ones, are treated as
 153 if they were of the same shape as the compound items with constant data
 154 values."
 155   (let ((first-compound (find-compound-data args)))
 156     (cond ((null first-compound) (apply fcn args))
 157           (t (fixup-map-elements-arglist args)
 158              (let* ((seq (compound-data-sequence first-compound))
 159                     (type (sequence-type seq)))
 160                (make-compound-data first-compound
 161                                    (apply #'fastmap type fcn args)))))))
 162
 163 (defun recursive-map-elements (base-fcn fcn &rest args)
 164 "Args: (base-fcn fcn &rest args)
 165 The same idea as MAP-ELEMENTS, except arguments are in a list and the
 166 base and recursive cases can use different functions. Modified to check
 167 for second level of compounding and use base-fcn if there is none."
 168   (let ((first-compound (find-compound-data args)))
 169     (cond ((null first-compound) (apply base-fcn args))
 170           (t (fixup-map-elements-arglist args)
 171              (let* ((seq (compound-data-sequence first-compound))
 172                     (type (sequence-type seq))
 173                     (f (if (any-compound-elements seq) fcn base-fcn)))
 174                (make-compound-data first-compound
 175                                    (apply #'fastmap type f args)))))))
 176
 177
 178 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 179 ;;;;
 180 ;;;;             Public Predicate and Accessor Functions
 181 ;;;;
 182 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 183
 184 (defun compound-data-p (x)
 185 "Args: (x)
 186 Returns T if X is a compound data item, NIL otherwise."
 187   (cmpndp x))
 188
 189 (defun compound-data-seq (x)
 190 "Args (x)
 191 Returns data sequence in X."
 192   (check-compound x)
 193   (compound-data-sequence x))
 194
 195 (defun compound-data-length (x)
 196 "Args (x)
 197 Returns length of data sequence in X."
 198   (check-compound x)
 199   (length (compound-data-sequence x)))
 200
 201 (defun element-list (x)
 202   (cond
 203    ((compound-data-p x)
 204     (let ((x (concatenate 'list (compound-data-seq x)))) ; copies sequence
 205       (cond
 206        ((any-compound-elements x)
 207         (do ((next x (rest next)))
 208             ((not (consp next)))
 209           (setf (first next) (element-list (first next))))
 210         (do ((result (first x))
 211              (last (last (first x)))
 212              (next (rest x) (rest next)))
 213             ((not (consp next)) result)
 214           (setf (rest last) (first next))
 215           (setf last (last (first next)))))
 216        (t x))))
 217    (t (list x))))
 218
 219 (defun element-seq (x)
 220 "Args: (x)
 221 Returns sequence of the elements of compound item X."
 222   (check-compound x)
 223   (let ((seq (compound-data-seq x)))
 224     (if (any-compound-elements seq) (element-list seq) seq)))
 225
 226 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 227 ;;;;
 228 ;;;;                     Compound Data Objects
 229 ;;;;
 230 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 231
 232 (defproto *compound-data-proto*)
 233
 234 (defmeth *compound-data-proto* :data-length (&rest args) nil)
 235 (defmeth *compound-data-proto* :data-seq (&rest args) nil)
 236 (defmeth *compound-data-proto* :make-data (&rest args) nil)
 237 (defmeth *compound-data-proto* :select-data (&rest args) nil)
 238
 239 (defun compound-object-p (x) (kind-of-p x *compound-data-proto*))
 240
 241
 242
 243 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 244 ;;;;
 245 ;;;;                         Sorting Functions
 246 ;;;;
 247 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 248
 249 (defun sort-data (x)
 250 "Args: (sequence)
 251 Returns a sequence with the numbers or strings in the sequence X in order."
 252   (flet ((less (x y) (if (numberp x) (< x y) (string-lessp x y))))
 253     (stable-sort (copy-seq (compound-data-seq x)) #'less)))
 254
 255 (defun order (x)
 256 "Args (x)
 257 Returns a sequence of the indices of elements in the sequence of numbers
 258 or strings X in order."
 259   (let* ((seq (compound-data-seq x))
 260          (type (if (consp seq) 'list 'vector))
 261          (i -1))
 262     (flet ((entry (x) (setf i (+ i 1)) (list x i))
 263            (less (a b)
 264                  (let ((x (first a))
 265                        (y (first b)))
 266                    (if (numberp x) (< x y) (string-lessp x y)))))
 267       (let ((sorted-seq (stable-sort (map type #'entry seq) #'less)))
 268         (map type #'second sorted-seq)))))
 269
 270 ;; this isn't destructive -- do we document destructive only, or any
 271 ;; variant?
 272 (defun rank (x)
 273 "Args (x)
 274 Returns a sequence with the elements of the list or array of numbers or
 275 strings X replaced by their ranks."
 276   (let ((ranked-seq (order (order x))))
 277     (make-compound-data
 278      ;; compound-data-shape is undefined?
 279      (compound-data-shape x) ranked-seq)))
 280
 281
 282
 283 ;;;
 284 ;;; REPEAT function
 285 ;;;
 286
 287 (defun repeat (a b)
 288 "Args: (vals times)
 289 Repeats VALS. If TIMES is a number and VALS is a non-null, non-array atom,
 290 a list of length TIMES with all elements eq to VALS is returned. If VALS
 291 is a list and TIMES is a number then VALS is appended TIMES times. If
 292 TIMES is a list of numbers then VALS must be a list of equal length and
 293 the simpler version of repeat is mapped down the two lists.
 294 Examples: (repeat 2 5)                 returns (2 2 2 2 2)
 295           (repeat '(1 2) 3)            returns (1 2 1 2 1 2)
 296           (repeat '(4 5 6) '(1 2 3))   returns (4 5 5 6 6 6)
 297           (repeat '((4) (5 6)) '(2 3)) returns (4 4 5 6 5 6 5 6)"
 298   (cond ((compound-data-p b)
 299          (let* ((reps (coerce (compound-data-seq (map-elements #'repeat a b))
 300                               'list))
 301                 (result (first reps))
 302                 (tail (last (first reps))))
 303            (dolist (next (rest reps) result)
 304                    (when next
 305                          (setf (rest tail) next)
 306                          (setf tail (last next))))))
 307         (t (let* ((a (if (compound-data-p a)
 308                          (coerce (compound-data-seq a) 'list)
 309                          (list a)))
 310                   (result nil))
 311              (dotimes (i b result)
 312                       (let ((next (copy-list a)))
 313                         (if result (setf (rest (last next)) result))
 314                         (setf result next)))))))
 315 ;;;
 316 ;;; WHICH function
 317 ;;;
 318
 319 (defun which (x)
 320 "Args: (x)
 321 Returns a list of the indices where elements of sequence X are not NIL."
 322   (let ((x (list (compound-data-seq x)))
 323         (result nil)
 324         (tail nil))
 325     (flet ((add-result (x)
 326              (if result (setf (rest tail) (list x)) (setf result (list x)))
 327              (setf tail (if tail (rest tail) result)))
 328            (get-next-element (seq-list i)
 329              (cond ((consp (first seq-list))
 330                     (let ((elem (first (first seq-list))))
 331                       (setf (first seq-list) (rest (first seq-list)))
 332                       elem))
 333                    (t (aref (first seq-list) i)))))
 334           (let ((n (length (first x))))
 335             (dotimes (i n result)
 336                      (if (get-next-element x i) (add-result i)))))))
 337
 338
 339 ;;; Sequences are part of ANSI CL, being a supertype of vector and
 340 ;;; list (ordered set of things).
 341 ;;;
 342 ;;; Need to use the interenal structure when possible -- silly to be
 343 ;;; redundant!  However, this means we need to understand what
 344 ;;; sequences were intending to do, which I'm not clear on yet.
 345
 346 ;;; The original ordering, object-wise, was to have compound
 347 ;;; functionality passed into sequences, into other data sources.
 348 ;;; However, at this point, we will see about inverting this and
 349 ;;; having basic data types pushed through compound, to simplify
 350 ;;; packaging.
 351
 352 ;;;                      Type Checking Functions
 353
 354 (defun check-sequence (a)
 355   ;; FIXME:AJR: does this handle consp as well?  (Luke had an "or"
 356   ;; with consp).
 357   (if (not (typep a 'sequence))
 358       (error "not a sequence - ~s" a)))
 359
 360 ;;;                       Sequence Element Access
 361
 362
 363 ;;; (elt x i) -- NOT.  This is more like "pop".
 364 (defun get-next-element (x i)
 365   "Get element i from seq x.  FIXME: not really??"
 366   (let ((myseq (first x)))
 367     (if (consp myseq)
 368         (let ((elem (first myseq)))
 369           (setf (first x) (rest myseq))
 370           elem)
 371       (aref myseq i))))
 372
 373 ;;; (setf (elt x i) v)
 374 (defun set-next-element (x i v)
 375   (let ((seq (first x)))
 376     (cond ((consp seq)
 377            (setf (first seq) v)
 378            (setf (first x) (rest seq)))
 379           (t (setf (aref seq i) v)))))
 380
 381 (defun make-next-element (x) (list x))
 382
 383
 384 ;;;                         Sequence Functions
 385
 386
 387 ;; to prevent breakage.
 388 (defmacro sequencep (x)
 389   (typep x 'sequence))
 390
 391 (defun iseq (a &optional b)
 392 "Args: (n &optional m)
 393 Generate a sequence of consecutive integers from a to b.
 394 With one argumant returns a list of consecutive integers from 0 to N - 1.
 395 With two returns a list of consecutive integers from N to M.
 396 Examples: (iseq 4) returns (0 1 2 3)
 397           (iseq 3 7)  returns (3 4 5 6 7)
 398           (iseq 3 -3) returns (3 2 1 0 -1 -2 -3)"
 399   (if b
 400       (let ((n (+ 1 (abs (- b a))))
 401             (x nil))
 402         (dotimes (i n x)
 403                  (setq x (cons (if (< a b) (- b i) (+ b i)) x))))
 404       (cond
 405        ((= 0 a) nil)
 406        ((< a 0) (iseq (+ a 1) 0))
 407        ((< 0 a) (iseq 0 (- a 1))))))
 408
 409 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 410 ;;;;
 411 ;;;;               Subset Selection and Mutation Functions
 412 ;;;;
 413 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 414
 415
 416 (defun old-rowmajor-index (index indices dim olddim)
 417   "translate row major index in resulting subarray to row major index
 418 in the original array."
 419   (declare (fixnum index))
 420   (let ((rank (length dim))
 421         (face 1)
 422         (oldface 1)
 423         (oldindex 0))
 424     (declare (fixnum rank face oldface))
 425
 426     (dotimes (i rank)
 427       (declare (fixnum i))
 428       (setf face (* face (aref dim i)))
 429       (setf oldface (* oldface (aref olddim i))))
 430
 431     (dotimes (i rank)
 432       (declare (fixnum i))
 433       (setf face (/ face (aref dim i)))
 434       (setf oldface (/ oldface (aref olddim i)))
 435       (incf oldindex
 436             (* oldface (aref (aref indices i) (floor (/ index face))))) ;;*** is this floor really needed???
 437       (setf index (rem index face)))
 438     oldindex))
 439
 440 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 441 ;;;;
 442 ;;;;               Subset Selection and Mutation Functions
 443 ;;;;
 444 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 445
 446 (defun subarray-select (a indexlist &optional (values nil set_values))
 447   "extract or set subarray for the indices from a displaced array."
 448   (let ((indices nil)
 449         (index)
 450         (dim)
 451         (vdim)
 452         (data)
 453         (result_data)
 454         (olddim)
 455         (result)
 456         (rank 0)
 457         (n 0)
 458         (k 0))
 459     (declare (fixnum rank n))
 460
 461     (if (or (sequencep a) (not (arrayp a))) (error "not an array - ~a" a))
 462     (if (not (listp indexlist))  (error "bad index list - ~a" indices))
 463     (if (/= (length indexlist)  (array-rank a))
 464         (error "wrong number of indices"))
 465
 466     (setf indices (coerce indexlist 'vector))
 467
 468     (setf olddim (coerce (array-dimensions a) 'vector))
 469
 470     ;; compute the result dimension vector and fix up the indices
 471     (setf rank (array-rank a))
 472     (setf dim (make-array rank))
 473     (dotimes (i rank)
 474       (declare (fixnum i))
 475       (setf index (aref indices i))
 476       (setf n (aref olddim i))
 477       (setf index (if (fixnump index) (vector index) (coerce index 'vector)))
 478       (setf k (length index))
 479       (dotimes (j k)
 480         (declare (fixnum j))
 481         (if (<= n (check-nonneg-fixnum (aref index j)))
 482           (error "index out of bounds - ~a" (aref index j)))
 483         (setf (aref indices i) index))
 484       (setf (aref dim i) (length index)))
 485
 486     ;; set up the result or check the values
 487     (let ((dim-list (coerce dim 'list)))
 488       (cond
 489        (set_values
 490         (cond
 491          ((compound-data-p values)
 492           (if (or (not (arrayp values)) (/= rank (array-rank values)))
 493             (error "bad values array - ~a" values))
 494           (setf vdim (coerce (array-dimensions values) 'vector))
 495           (dotimes (i rank)
 496             (declare (fixnum i))
 497             (if (/= (aref vdim i) (aref dim i))
 498               (error "bad value array dimensions - ~a" values)))
 499           (setf result values))
 500          (t (setf result (make-array dim-list :initial-element values)))))
 501        (t (setf result (make-array dim-list)))))
 502
 503     ;; compute the result or set the values
 504     (setf data (compound-data-seq a))
 505     (setf result_data (compound-data-seq result))
 506     (setf n (length result_data))
 507     (dotimes (i n)
 508       (declare (fixnum i))
 509       (setf k (old-rowmajor-index i indices dim olddim))
 510       (if (or (> 0 k) (>= k (length data))) (error "index out of range"))
 511       (if set_values
 512         (setf (aref data k) (aref result_data i))
 513         (setf (aref result_data i) (aref data k))))
 514
 515     result))
 516
 517
 518 ;;;; is x an ordered sequence of nonnegative positive integers?
 519 (defun ordered-nneg-seq(x)
 520   ;; FIXME -- sbcl warning about unreachable code, might be a logic error here.
 521   (if (sequencep x)
 522       (let ((n (length x))
 523             (cx (make-next-element x))
 524             (m 0))
 525         (dotimes (i n t)
 526           (let ((elem (check-nonneg-fixnum (get-next-element cx i))))
 527             (if (> m elem) (return nil) (setf m elem)))))))
 528
 529 ;;;; select or set the subsequence corresponding to the specified indices
 530 (defun sequence-select(x indices &optional (values nil set-values))
 531   ;; FIXME -- sbcl warning about unreachable code, might be a logic error here.
 532   (let ((rlen 0)
 533         (dlen 0)
 534         (vlen 0)
 535         (data nil)
 536         (result nil))
 537     (declare (fixnum rlen dlen vlen))
 538
 539     ;; Check the input data
 540     (check-sequence x)
 541     (check-sequence indices)
 542     (if set-values (check-sequence values))
 543
 544     ;; Find the data sizes
 545     (setf data (if (ordered-nneg-seq indices) x (coerce x 'vector)))
 546     (setf dlen (length data))
 547     (setf rlen (length indices))
 548     (when set-values
 549       (setf vlen (length values))
 550       (if (/= vlen rlen) (error "value and index sequences do not match")))
 551
 552     ;; set up the result/value sequence
 553     (setf result
 554           (if set-values
 555               values
 556               (make-sequence (if (listp x) 'list 'vector) rlen)))
 557
 558     ;; get or set the sequence elements
 559     (if set-values
 560       (do ((nextx x)
 561            (cr (make-next-element result))
 562            (ci (make-next-element indices))
 563            (i 0 (+ i 1))
 564            (j 0)
 565            (index 0))
 566           ((>= i rlen))
 567         (declare (fixnum i j index))
 568         (setf index (get-next-element ci i))
 569         (if (<= dlen index) (error "index out of range - ~a" index))
 570         (let ((elem (get-next-element cr i)))
 571           (cond
 572            ((listp x)
 573             (when (> j index)
 574               (setf j 0)
 575               (setf nextx x))
 576             (do ()
 577                 ((not (and (< j index) (consp nextx))))
 578               (incf j 1)
 579               (setf nextx (rest nextx)))
 580             (setf (first nextx) elem))
 581            (t (setf (aref x index) elem)))))
 582       (do ((nextx data)
 583            (cr (make-next-element result))
 584            (ci (make-next-element indices))
 585            (i 0 (+ i 1))
 586            (j 0)
 587            (index 0)
 588            (elem nil))
 589           ((>= i rlen))
 590         (declare (fixnum i j index))
 591         (setf index (get-next-element ci i))
 592         (if (<= dlen index) (error "index out of range - ~a" index))
 593         (cond
 594          ((listp data) ;; indices must be ordered
 595           (do ()
 596               ((not (and (< j index) (consp nextx))))
 597             (incf j 1)
 598             (setf nextx (rest nextx)))
 599           (setf elem (first nextx)))
 600          (t (setf elem (aref data index))))
 601         (set-next-element cr i elem)))
 602
 603     result))
 604
 605 ;;;
 606 ;;; SELECT function
 607 ;;;
 608
 609 (defun select (x &rest args)
 610 "Args: (a &rest indices)
 611 A can be a list or an array. If A is a list and INDICES is a single number
 612 then the appropriate element of A is returned. If  is a list and INDICES is
 613 a list of numbers then the sublist of the corresponding elements is returned.
 614 If A in an array then the number of INDICES must match the ARRAY-RANK of A.
 615 If each index is a number then the appropriate array element is returned.
 616 Otherwise the INDICES must all be lists of numbers and the corresponding
 617 submatrix of A is returned. SELECT can be used in setf."
 618   (cond
 619    ((every #'fixnump args)
 620     (if (listp x) (nth (first args) x) (apply #'aref x args)))
 621    ((sequencep x) (sequence-select x (first args)))
 622    (t (subarray-select x args))))
 623
 624
 625 ;; Built in SET-SELECT (SETF method for SELECT)
 626 (defun set-select (x &rest args)
 627   (let ((indices (butlast args))
 628         (values (first (last args))))
 629     (cond
 630      ((sequencep x)
 631       (if (not (consp indices)) (error "bad indices - ~a" indices))
 632       (let* ((indices (first indices))
 633              (i-list (if (fixnump indices) (list indices) indices))
 634              (v-list (if (fixnump indices) (list values) values)))
 635         (sequence-select x i-list v-list)))
 636      ((arrayp x)
 637       (subarray-select x indices values))
 638      (t (error "bad argument type - ~a" x)))
 639     values))
 640
 641 (defsetf select set-select)
 642
 643 ;;;;
 644 ;;;; Basic Sequence Operations
 645 ;;;;
 646
 647 (defun difference (x)
 648 "Args: (x)
 649 Returns differences for a sequence X."
 650   (let ((n (length x)))
 651     (- (select x (iseq 1 (1- n))) (select x (iseq 0 (- n 2))))))
 652
 653 (defun rseq (a b num)
 654 "Args: (a b num)
 655 Returns a list of NUM equally spaced points starting at A and ending at B."
 656   (+ a (* (values-list (iseq 0 (1- num))) (/ (float (- b a)) (1- num)))))
 657
 658