Changeset 28194 in project


Ignore:
Timestamp:
01/30/13 08:16:33 (8 years ago)
Author:
Ivan Raikov
Message:

srfi-4-utils: moving repmat and meshgrid to matrix-utils

Location:
release/4
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • release/4/matrix-utils/trunk/matrix-utils.scm

    r26824 r28194  
    44;;
    55;;
    6 ;; Copyright 2007-2012 Ivan Raikov and the Okinawa Institute of
     6;; Copyright 2007-2013 Ivan Raikov and the Okinawa Institute of
    77;; Science and Technology.
    88;;
     
    3838                 
    3939                 define-utility-matrices
    40                  with-utility-matrices)
    41 
    42    (import scheme chicken data-structures srfi-4)
     40                 with-utility-matrices
     41
     42                 f64vector-repmat f64vector-meshgrid
     43                 )
     44
     45   (import scheme chicken foreign data-structures srfi-4)
    4346   (require-extension srfi-4 srfi-42 srfi-4-comprehensions blas)
    4447
     
    428431  )
    429432
     433
     434#>
     435
     436/* repeat a block of memory rep times */
     437void cmemrep(double *dest, size_t chunk, int rep)
     438{
     439  if(rep == 1) return;
     440  memcpy(dest + chunk, dest, sizeof(double)*chunk);
     441  if(rep & 1) {
     442    dest += chunk;
     443    memcpy(dest + chunk, dest, sizeof(double)*chunk);
     444  }
     445  /* now repeat using a block twice as big */
     446  cmemrep(dest, chunk<<1, rep>>1);
     447}
     448
     449void crepmat(double *dest, const double *src, int ndim,
     450             int *destdimsize, int *dimsize, const int *dims, int *rep)
     451{
     452  int i, chunk;
     453  int d = ndim-1;
     454  /* copy the first repetition into dest */
     455  if(d == 0)
     456  {
     457    chunk = dimsize[0];
     458    memcpy(dest,src,sizeof(double)*chunk);
     459  }
     460  else
     461  {
     462    /* recursively repeat each slice of src */
     463    for(i=0; i<dims[d]; i++) {
     464      crepmat(dest + i*destdimsize[d-1], src + i*dimsize[d-1],
     465              ndim-1, destdimsize, dimsize, dims, rep);
     466    }
     467    chunk = destdimsize[d-1]*dims[d];
     468  }
     469  /* copy the result rep-1 times */
     470  cmemrep(dest,chunk,rep[d]);
     471}
     472
     473<#
     474
     475(define crepmat (foreign-safe-lambda void "crepmat"  f64vector f64vector int s32vector s32vector s32vector s32vector))
     476
     477(define cmemrep (foreign-safe-lambda void "cmemrep"  f64vector int int))
     478
     479(define (f64vector-repmat src dims reps)
     480  (let* ((ndim (length dims))
     481         (dimsize (make-s32vector ndim))
     482         (vdims (list->s32vector dims)))
     483   
     484    (assert (every positive? reps))
     485   
     486    ;; compute dimension sizes
     487    (s32vector-set! dimsize 0 (car dims))
     488    (let recur ((i 1) (dims (cdr dims)))
     489      (if (< i ndim)
     490          (begin
     491            (s32vector-set! dimsize i (* (car dims) (s32vector-ref dimsize (- i 1))))
     492            (recur (+ 1 i) (cdr dims)))
     493          ))
     494
     495    ;; determine repetition vector
     496    (let* ((nrep (length reps))
     497           (ndimdest (if (> nrep ndim) nrep ndim))
     498           (rep (make-s32vector ndimdest)))
     499      (let recur ((i 0) (reps reps))
     500        (if (< i nrep)
     501            (let ((repv (car reps)))
     502              (s32vector-set! rep i repv)
     503              (recur (+ 1 i) (cdr reps)))))
     504
     505      ;; compute output size
     506      (let ((destdims (make-s32vector ndimdest 0)))
     507        (let recur ((i 0))
     508          (if (< i ndim)
     509              (begin (s32vector-set! destdims i
     510                                    (* (s32vector-ref vdims i)
     511                                       (s32vector-ref rep i )))
     512                     (recur (+ 1 i)))))
     513
     514        (let ((extrarep
     515               (let recur ((i ndim) (extrarep 1))
     516                 (if (< i ndimdest)
     517                     (let ((ri (s32vector-ref rep i)))
     518                       (s32vector-set! destdims i ri)
     519                       (recur (+ 1 i) (* extrarep ri)))
     520                     extrarep
     521                     ))))
     522
     523          (let ((destdimsize (make-s32vector ndimdest)))
     524            (s32vector-set! destdimsize 0 (s32vector-ref destdims 0))
     525            (let recur ((i 1))
     526              (if (< i ndimdest)
     527                  (begin
     528                    (s32vector-set! destdimsize i
     529                                    (* (s32vector-ref destdimsize (- i 1))
     530                                       (s32vector-ref destdims i)))
     531                    (recur (+ 1 i)))))
     532
     533            ;; return array
     534            (let* ((destlen (s32vector-ref destdimsize (- ndimdest 1)))
     535                   (dest (make-f64vector destlen)))
     536
     537              (crepmat dest src ndim destdimsize dimsize vdims rep)
     538              (if (> ndimdest ndim)
     539                  (let ((n (s32vector-ref destdimsize (- ndim 1))))
     540                    (cmemrep dest n extrarep)))
     541              dest
     542              ))
     543          ))
     544      )))
     545
     546
     547         
     548
     549(define (f64vector-meshgrid x y z)
     550  (let ((lenx (f64vector-length x))
     551        (leny (f64vector-length y))
     552        (lenz (f64vector-length z)))
     553    (let ((dimx (list 1 lenx))
     554          (dimy (list 1 leny))
     555          (dimz (list 1 lenz)))
     556    (let ((xx (f64vector-repmat
     557               (f64vector-repmat x dimx (list leny 1))
     558               (list leny lenx)
     559               (list 1 1 lenz)))
     560          (yy (f64vector-repmat
     561               (f64vector-repmat y dimy (list 1 lenx))
     562               (list lenx leny)
     563               (list 1 1 lenz)))
     564          (zz (f64vector-repmat z dimz (list (* lenx leny) 1))))
     565      (list xx yy zz))
     566    )))
     567     
     568
    430569)
    431        
    432 
    433 
    434 
     570
     571
     572
  • release/4/srfi-4-utils/trunk/srfi-4-utils.meta

    r23287 r28194  
    1616 (category data)
    1717
    18  ; A list of eggs matrix-utils depends on.
     18 ; A list of eggs srfi-4-utils depends on.
    1919
    2020 (needs srfi-42 srfi-4-comprehensions )
     21
     22 (test-depends test)
    2123
    2224 (author "Ivan Raikov")
  • release/4/srfi-4-utils/trunk/srfi-4-utils.scm

    r28189 r28194  
    4646                 s16vector-merge-sort! u16vector-merge-sort! s8vector-merge-sort! u8vector-merge-sort!
    4747
    48                  f64vector-repmat f64vector-meshgrid
    4948                 )
    5049
    51  (import scheme chicken foreign)
     50 (import scheme chicken)
    5251 (require-library srfi-1 extras)
    5352 (import (only srfi-1 fold every)
     
    157156
    158157(define (srfi-4-vector-quick-sort! vector-ref vector-set! vector-length)
    159   (lambda (v elt< . rest)
     158  (lambda (elt< v . rest)
    160159   (let-optionals rest ((start 0) (end (vector-length v)))
    161160    (let recur ((l start) (r end))      ; Sort the range [l,r).
     
    331330(define u8vector-merge-sort!
    332331  (srfi-4-vector-merge-sort! u8vector-length u8vector-ref u8vector-set! u8vector-merge! u8vector-blit!))
    333                    
    334          
    335 
    336 
    337 
    338 #>
    339 
    340 /* repeat a block of memory rep times */
    341 void cmemrep(double *dest, size_t chunk, int rep)
    342 {
    343   if(rep == 1) return;
    344   memcpy(dest + chunk, dest, sizeof(double)*chunk);
    345   if(rep & 1) {
    346     dest += chunk;
    347     memcpy(dest + chunk, dest, sizeof(double)*chunk);
    348   }
    349   /* now repeat using a block twice as big */
    350   cmemrep(dest, chunk<<1, rep>>1);
    351 }
    352 
    353 void crepmat(double *dest, const double *src, int ndim,
    354              int *destdimsize, int *dimsize, const int *dims, int *rep)
    355 {
    356   int i, chunk;
    357   int d = ndim-1;
    358   /* copy the first repetition into dest */
    359   if(d == 0)
    360   {
    361     chunk = dimsize[0];
    362     memcpy(dest,src,sizeof(double)*chunk);
    363   }
    364   else
    365   {
    366     /* recursively repeat each slice of src */
    367     for(i=0; i<dims[d]; i++) {
    368       crepmat(dest + i*destdimsize[d-1], src + i*dimsize[d-1],
    369               ndim-1, destdimsize, dimsize, dims, rep);
    370     }
    371     chunk = destdimsize[d-1]*dims[d];
    372   }
    373   /* copy the result rep-1 times */
    374   cmemrep(dest,chunk,rep[d]);
    375 }
    376 
    377 <#
    378 
    379 (define crepmat (foreign-safe-lambda void "crepmat"  f64vector f64vector int s32vector s32vector s32vector s32vector))
    380 
    381 (define cmemrep (foreign-safe-lambda void "cmemrep"  f64vector int int))
    382 
    383 (define (f64vector-repmat src dims reps)
    384   (let* ((ndim (length dims))
    385          (dimsize (make-s32vector ndim))
    386          (vdims (list->s32vector dims)))
    387    
    388     (assert (every positive? reps))
    389    
    390     ;; compute dimension sizes
    391     (s32vector-set! dimsize 0 (car dims))
    392     (let recur ((i 1) (dims (cdr dims)))
    393       (if (< i ndim)
    394           (begin
    395             (s32vector-set! dimsize i (* (car dims) (s32vector-ref dimsize (- i 1))))
    396             (recur (+ 1 i) (cdr dims)))
    397           ))
    398 
    399     ;; determine repetition vector
    400     (let* ((nrep (length reps))
    401            (ndimdest (if (> nrep ndim) nrep ndim))
    402            (rep (make-s32vector ndimdest)))
    403       (let recur ((i 0) (reps reps))
    404         (if (< i nrep)
    405             (let ((repv (car reps)))
    406               (s32vector-set! rep i repv)
    407               (recur (+ 1 i) (cdr reps)))))
    408 
    409       ;; compute output size
    410       (let ((destdims (make-s32vector ndimdest 0)))
    411         (let recur ((i 0))
    412           (if (< i ndim)
    413               (begin (s32vector-set! destdims i
    414                                     (* (s32vector-ref vdims i)
    415                                        (s32vector-ref rep i )))
    416                      (recur (+ 1 i)))))
    417 
    418         (let ((extrarep
    419                (let recur ((i ndim) (extrarep 1))
    420                  (if (< i ndimdest)
    421                      (let ((ri (s32vector-ref rep i)))
    422                        (s32vector-set! destdims i ri)
    423                        (recur (+ 1 i) (* extrarep ri)))
    424                      extrarep
    425                      ))))
    426 
    427           (let ((destdimsize (make-s32vector ndimdest)))
    428             (s32vector-set! destdimsize 0 (s32vector-ref destdims 0))
    429             (let recur ((i 1))
    430               (if (< i ndimdest)
    431                   (begin
    432                     (s32vector-set! destdimsize i
    433                                     (* (s32vector-ref destdimsize (- i 1))
    434                                        (s32vector-ref destdims i)))
    435                     (recur (+ 1 i)))))
    436 
    437             ;; return array
    438             (let* ((destlen (s32vector-ref destdimsize (- ndimdest 1)))
    439                    (dest (make-f64vector destlen)))
    440 
    441               (crepmat dest src ndim destdimsize dimsize vdims rep)
    442               (if (> ndimdest ndim)
    443                   (let ((n (s32vector-ref destdimsize (- ndim 1))))
    444                     (cmemrep dest n extrarep)))
    445               dest
    446               ))
    447           ))
    448       )))
    449 
    450 
    451          
    452 
    453 (define (f64vector-meshgrid x y z)
    454   (let ((lenx (f64vector-length x))
    455         (leny (f64vector-length y))
    456         (lenz (f64vector-length z)))
    457     (let ((dimx (list 1 lenx))
    458           (dimy (list 1 leny))
    459           (dimz (list 1 lenz)))
    460     (let ((xx (f64vector-repmat
    461                (f64vector-repmat x dimx (list leny 1))
    462                (list leny lenx)
    463                (list 1 1 lenz)))
    464           (yy (f64vector-repmat
    465                (f64vector-repmat y dimy (list 1 lenx))
    466                (list lenx leny)
    467                (list 1 1 lenz)))
    468           (zz (f64vector-repmat z dimz (list (* lenx leny) 1))))
    469       (list xx yy zz))
    470     )))
    471      
    472 
    473332
    474333)
  • release/4/srfi-4-utils/trunk/tests/run.scm

    r28191 r28194  
    44(define elt< (lambda (i ai j bj) (< ai bj)))
    55(define elt> (lambda (i ai j bj) (> ai bj)))
    6 
    7 (print (f64vector-merge! elt< (f64vector 1 3 5 2 4) 0 3 5 (make-f64vector 5 0.) 0))
    8 
    9 (print (f64vector-merge! elt< (f64vector 1 3 5 2 4 6) 0 3 6 (make-f64vector 6 0.) 0))
    106
    117(test-group
Note: See TracChangeset for help on using the changeset viewer.