source: project/release/3/mat5-lib/trunk/mat5-lib.scm @ 11550

Last change on this file since 11550 was 11550, checked in by Ivan Raikov, 12 years ago

Updates to license text; eliminated use of syntax-case.

File size: 84.9 KB
Line 
1;;
2;;
3;; Definitions and read/write routines for MAT 5.0 binary format.
4;;
5;; Copyright 2005-2008 Ivan Raikov and the Georgia Institute of Technology
6;;
7;; This program is free software: you can redistribute it and/or
8;; modify it under the terms of the GNU General Public License as
9;; published by the Free Software Foundation, either version 3 of the
10;; License, or (at your option) any later version.
11;;
12;; This program is distributed in the hope that it will be useful, but
13;; WITHOUT ANY WARRANTY; without even the implied warranty of
14;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15;; General Public License for more details.
16;;
17;; A full copy of the GPL license can be found at
18;; <http://www.gnu.org/licenses/>.
19;;
20
21(require-extension iset)
22(require-extension srfi-1)
23(require-extension srfi-4)
24(require-extension srfi-13)
25(require-extension srfi-14)
26(require-extension srfi-40)
27(require-extension srfi-47)
28(require-extension z3)
29(require-extension datatype)
30(require-extension endian-port)
31
32(require-extension lolevel)
33(require-extension extras)
34
35(define-extension mat5-lib)
36
37
38(declare (export MAT5:data-type?                 
39                 miINT8 miUINT8 miINT16 miUINT16 miINT32 miUINT32 
40                 miSINGLE miDOUBLE miINT64 miUINT64 miMATRIX 
41                 miCOMPRESSED miUTF8 miUTF16 miUTF32
42                 make-MAT5:header MAT5:header?
43                 MAT5:header-magic MAT5:header-text MAT5:header-subsys
44                 MAT5:header-version MAT5:header-eport
45                 make-MAT5:data-element MAT5:data-element?
46                 MAT5:data-element-type MAT5:data-element-bytes MAT5:data-element-data 
47                 MAT5:array? MAT5:object MAT5:structure MAT5:cell-array MAT5:sparse-array MAT5:num-array
48                 MAT5:dimensions?
49                 MAT5:numeric-type? 
50                 MAT5:array-type? 
51                 init-MAT5:cell MAT5:cell? MAT5:cell-data MAT5:cell-dims
52                 vector->MAT5:cell
53                 MAT5:cell-dims 
54                 MAT5:cell-ref 
55                 MAT5:cell-set!
56                 MAT5:array-foldi
57                 MAT5:read-header               
58                 MAT5:read-data-element
59                 MAT5:write-header
60                 MAT5:write-data-element
61                 MAT5:read-stream               
62                 MAT5:stream-read-data-element           
63                 MAT5:write-stream               
64                 MAT5:stream-write-data-element)
65         (inline)
66         (lambda-lift))
67                 
68
69;--------------------
70;  Message routines
71;
72;
73
74(define (MAT5:warning x . rest)
75  (let loop ((port (open-output-string)) (objs (cons x rest)))
76    (if (null? objs)
77        (begin
78          (newline port)
79          (print-error-message (get-output-string port) 
80                               (current-error-port) "MAT5 warning:"))
81        (begin (display (car objs) port)
82               (display " " port)
83               (loop port (cdr objs))))))
84
85(define (MAT5:error x . rest)
86  (let ((port (open-output-string)))
87    (if (endian-port? x)
88        (begin
89          (display "[" port)
90          (display (endian-port-filename x) port)
91          (display ":" port)
92          (display (endian-port-pos x) port)
93          (display "] " port)))
94    (let loop ((objs (if (endian-port? x) rest (cons x rest))))
95      (if (null? objs)
96          (begin
97            (newline port)
98            (error 'MAT5 (get-output-string port)))
99          (begin (display (car objs) port)
100                 (display " " port)
101                 (loop (cdr objs)))))))
102
103
104;----------------------------------
105;  General-purpose utility routines
106;
107
108(define (split2 x) (values (first x) (second x)))
109(define (split3 x) (values (first x) (second x) (third x)))
110
111(define (byte-vector->zstring x)
112  (let ((ws  (char-set-union char-set:whitespace (char-set (integer->char 0))))
113        (s   (byte-vector->string x)))
114    (string-trim-right s ws)))
115
116(define (natural? x) (and (integer? x) (positive? x)))
117
118(define (natural-list? x)
119  (and (list? x) (every natural? x)))
120
121(define (positive-or-zero? x) (and (integer? x) (or (zero? x) (positive? x))))
122
123(define (positive-or-zero-list? x)
124  (and (list? x) (every positive-or-zero? x)))
125
126
127(define (string-list? x)
128  (and (list? x) (every string? x)))
129
130
131; Procedure:
132; stream-fold:: FUNC AX STRM -> AX
133;
134; Taken from stream-ext library by Alejandro Forero Cuervo
135;
136(define (stream-fold func nil str)
137  (if (stream-null? str) nil
138    (stream-fold func (func (stream-car str) nil) (stream-cdr str))))
139
140
141; Procedure:
142; vector->array:: VECTOR-OPS VECTOR ARRAY-PROTOTYPE DIMENSIONS -> ARRAY
143;
144; Based on vector->array from SRFI-63 reference implementation
145; Copyright (C) 2001, 2003, 2005 Aubrey Jaffer
146;
147(define (vector->array vops vect prototype dimensions . rest)
148  (let-optionals  rest ((order 'row-major))
149   (let* ((vector-length (alist-ref 'vector-length vops))
150          (vector-ref (alist-ref 'vector-ref vops))
151          (vdx (vector-length vect)))
152     (if (not (eqv? vdx (apply * dimensions))) (MAT5:error "incompatible dimensions"))
153     (letrec ((ra (apply make-array prototype dimensions))
154              (v2ra  (lambda (dims idxs)
155                       (cond ((null? dims)
156                              (set! vdx (+ -1 vdx))
157                              (apply array-set! ra (vector-ref vect vdx) (reverse idxs)))
158                             (else
159                              (do ((idx (+ -1 (car dims)) (+ -1 idx)))
160                                  ((negative? idx) vect)
161                                (v2ra (cdr dims) (cons idx idxs)))))))
162              (v2ca   (lambda (dims idxs)
163                        (cond ((null? dims)
164                               (set! vdx (+ -1 vdx))
165                               (apply array-set! ra (vector-ref vect vdx) idxs))
166                              (else
167                               (do ((idx (+ -1 (car dims)) (+ -1 idx)))
168                                   ((negative? idx) vect)
169                                 (v2ca (cdr dims) (cons idx idxs))))))))
170       (if (eq? order 'row-major)
171           (v2ra dimensions '())
172           (v2ca (reverse dimensions) '()))
173         ra))))
174
175
176; Procedure:
177; align-eport-pos:: ENDIAN-PORT -> UNDEFINED
178;
179; Set file position to 64-bit boundary. Returns the # bytes skipped
180;
181;   * eport    an endian port
182;
183(define (align-eport-pos eport)
184  (let* ((pos  (endian-port-pos eport))
185         (pos1  (let loop ((pos1 pos))
186                   (if (not (zero? (modulo pos1 8))) 
187                       (loop (+ pos1 1)) pos1))))
188    (endian-port-setpos eport pos1)
189    (- pos1 pos)))
190
191
192
193;-------------------------------------
194;  MAT-File constants
195;
196
197(define MAT5:header-text-size    116)
198(define MAT5:header-subsys-size    8)
199
200
201; Version and magic numbers
202(define MAT5:lsb-version  #x0001)
203(define MAT5:lsb-magic    #x494D)
204(define MAT5:msb-version  #x0100)
205(define MAT5:msb-magic    #x4D49)
206
207;
208; MAT-File array flags (represented as bit positions)
209;
210(define mxLOGICAL_FLAG  1)
211(define mxGLOBAL_FLAG   2)
212(define mxCOMPLEX_FLAG  3)
213
214(define MAT5:field-name-length 32)
215
216;-------------------------------------
217;  Structure and datatype definitions
218;
219
220; Datatype: MAT5:data-type
221;
222; MAT-file data types
223(define-datatype MAT5:data-type MAT5:data-type? 
224  (miINT8)
225  (miUINT8)
226  (miINT16)     
227  (miUINT16)     
228  (miINT32)     
229  (miUINT32)     
230  (miSINGLE)     
231  (miDOUBLE)     
232  (miINT64)     
233  (miUINT64)     
234  (miMATRIX)     
235  (miCOMPRESSED)
236  (miUTF8)       
237  (miUTF16)     
238  (miUTF32))     
239
240;
241; MAT-File array classes
242;
243(define-datatype MAT5:array-class MAT5:array-class? 
244  (mxCELL_CLASS)
245  (mxSTRUCT_CLASS)
246  (mxOBJECT_CLASS)
247  (mxCHAR_CLASS)
248  (mxSPARSE_CLASS)
249  (mxDOUBLE_CLASS)
250  (mxSINGLE_CLASS)
251  (mxINT8_CLASS)
252  (mxUINT8_CLASS)
253  (mxINT16_CLASS)
254  (mxUINT16_CLASS)
255  (mxINT32_CLASS)
256  (mxUINT32_CLASS))
257
258
259;  Structure: MAT5:header
260;
261;  MAT-file Level 5 Header Format
262
263;  * text     is a byte vector of size 116
264;  * subsys   is a byte vector of size 8
265;  * version  is unsigned 16-bit version number
266;  * eport    is an endian port (see module endian-port)
267;
268(define-record MAT5:header  magic text subsys version eport)
269
270
271;  Structure: MAT5:data-element
272;
273;  MAT-file Level 5 data element 
274
275;  * type     is of type MAT5:data-type
276;  * bytes    is the size of the data element in bytes
277;  * data     is the data element object -- either a fixnum/flonum or an array
278;
279(define-record MAT5:data-element  type bytes data)
280
281
282;  Structure: MAT5:array-flags
283;
284;  MAT-file Level 5 array flags
285;
286;  * flags    a bit vector that contains the mxLOGICAL_FLAG, mxGLOBAL_FLAG,
287;             mxCOMPLEX_FLAG flags, if set.
288;  * class    array class; see the definition of MAT5:array-class datatype
289;  * nzmax    maximum number of non-zero elements in a sparse matrix
290;
291(define-record MAT5:array-flags  flags  class  nzmax)
292
293
294; Datatype: MAT5:array
295;
296; A representation of the different types of MAT-file arrays
297;
298(define-datatype MAT5:array MAT5:array? 
299
300
301  ;
302  ; Object
303  ;
304  ;  * dims        dimensions (a list of positive integers) -- these are the
305  ;                dimensions of the cells contained in the structure, if any
306  ;  * class       class name (string)
307  ;  * fields      a cell in which every element is an alist of  fields
308  ;
309  (MAT5:object      (name         string?)
310                    (dims         natural-list?)
311                    (class-name   string?)
312                    (field-names  string-list?)
313                    (fields       MAT5:cell?))
314
315  ; Structure
316  ;
317  ;  * dims        dimensions (a list of positive integers) -- these are the
318  ;                dimensions of the cells contained in the structure, if any
319  ;  * fields      a cell in which every element is an alist of  fields
320  ;
321  (MAT5:structure   (name         string?)
322                    (dims         natural-list?)
323                    (field-names  string-list?)
324                    (fields       MAT5:cell?))
325
326  ;
327  ; Cell array
328  ;
329  ;  * dims        dimensions (a list of positive integers)
330  ;  * cell        a MAT5:cell object (nested R5RS vectors)
331  ;
332  (MAT5:cell-array    (name       string?)
333                      (dims       natural-list?)
334                      (cell       MAT5:cell?))
335
336  ;
337  ; Homogeneous sparse array
338  ;
339  ;  * data-type   array element MAT5 type (can only be numeric type)
340  ;  * dims        dimensions (a list of positive integers)
341  ;  * row-index   row indices of non-zero elements
342  ;  * col-index   column indices of non-zero elements
343  ;  * real        real part (an SRFI-47 array or a stream of column vectors)
344  ;  * imag        imaginary part (array or #f)
345  ;
346  (MAT5:sparse-array    (name       string?)
347                        (data-type  MAT5:numeric-type?)
348                        (dims       MAT5:dimensions?)
349                        (row-index  natural-list?)
350                        (col-index  natural-list?)
351                        (real       (lambda (x) (or (array? x) (stream? x)))) 
352                        (imag       (lambda (x) (or (array? x) (stream? x) (not x)))))
353
354  ;
355  ; Homogeneous numeric array
356  ;
357  ;  * data-type  array element MAT5 type (can only be numeric type)
358  ;  * dims       dimensions (a list of positive integers)
359  ;  * real       real part (an SRFI-47 array or a stream of column vectors)
360  ;  * imag       imaginary part (array, stream, or #f)
361  ;
362  (MAT5:num-array  (name       string?)
363                   (data-type  MAT5:numeric-type?)
364                   (dims       MAT5:dimensions?) 
365                   (real       (lambda (x) (or (array? x) (stream? x))))
366                   (imag       (lambda (x) (or (array? x) (stream? x) (not x))))))
367
368;
369; The dimensions of a MAT5 object/struct/cell/array are defined as a
370; list of positive integers, OR a list of positive integers where the
371; last element is the symbol ??  The second format is used when the
372; array data is represented by a stream of column vectors, and it is
373; not known how many elements are in the stream ahead of time.
374;
375(define (MAT5:dimensions? x)
376  (let* ((rx   (reverse x))
377         (last (car rx)))
378    (cond ((natural? last) (natural-list? (cdr rx)))
379          ((eq? '?? last)  (natural-list? (cdr rx)))
380          (else #f))))
381
382
383
384;-----------------------
385; Record pretty printers
386;
387
388(define-record-printer (array x out)
389  (define m 5)
390  (let* ((dims  (array-dimensions x))
391         (vdx  (last dims)))
392    (letrec ((arpp  (lambda (dims idxs lst)
393                      (cond ((null? (cdr dims))
394                             (let ((v  (list-tabulate (if (> m vdx) vdx m) 
395                                                      (lambda (i) 
396                                                        (apply array-ref x (reverse (cons i idxs)))))))
397                             (cons (if (< m vdx) (append v '(...)) v) lst)))
398                            (else
399                             (let ((dim  (car dims))
400                                   (lst1 '()))
401                               (do ((idx 0 (+ 1 idx)))
402                                   ((or (>= idx dim) (>= idx m)) 
403                                    (if (>= idx m) 
404                                        (cons (reverse (cons '... lst1)) lst)
405                                        (cons (reverse lst1) lst)))
406                                 (set! lst1 (arpp (cdr dims) (cons idx idxs) lst1)))))))))
407
408      (fprintf out "#(array ~S)"
409               (reverse (arpp dims '() '()))))))
410   
411 
412
413(define-record-printer (MAT5:header x out)
414  (fprintf out "#(MAT5:header magic=0x~X text=~S subsys=~S version=~S)"
415           (MAT5:header-magic x)
416           (MAT5:header-text x) 
417           (MAT5:header-subsys x) 
418           (MAT5:header-version x)) )
419
420
421(define-record-printer (MAT5:data-element x out)
422  (fprintf out "#(MAT5:data-element type=~S bytes=~S data=~S)"
423           (MAT5:data-element-type x)
424           (MAT5:data-element-bytes x)
425           (MAT5:data-element-data x)))
426
427(define-record-printer (MAT5:data-type x out)
428  (fprintf out "~A" (MAT5:data-type->string x)))
429
430(define-record-printer (MAT5:array-class x out)
431  (fprintf out "~A" (MAT5:array-class->string x)))
432
433(define-record-printer (MAT5:array x out)
434  (cases MAT5:array x
435         (MAT5:object (name dims class field-names fields)
436                      (fprintf out "#(MAT5:object name=~S dims=~S class=~S field-names=~A fields=~S)"
437                               name dims class field-names fields))
438
439         (MAT5:structure (name dims field-names fields)
440                         (fprintf out "#(MAT5:structure name=~S dims=~S field-names=~A fields=~S)"
441                                  name dims field-names fields))
442         
443         (MAT5:cell-array (name dims cell)
444                          (fprintf out "#(MAT5:cell name=~S dims=~S cell=~S)"
445                                   name dims cell))
446         
447         (MAT5:sparse-array (name data-type dims row-index col-index real imag) 
448          (fprintf out "#(MAT5:sparse-array name=~S data-type=~A dims=~S row-index=~S col-index=~S real=~S imag=~S)"
449                   name data-type dims row-index col-index real imag))
450         
451         (MAT5:num-array (name data-type dims  real imag)
452                         (fprintf out "#(MAT5:num-array name=~S data-type=~A dims=~S real=~S imag=~S)"
453                                  name data-type dims  real imag))
454         
455         (else           "#<MAT5:array>")))
456
457(define-record-printer (MAT5:array-flags x out)
458  (let ((flags (MAT5:array-flags-flags x))
459        (class (MAT5:array-flags-class x))
460        (nzmax (MAT5:array-flags-nzmax x)))
461    (let ((sflags (filter (lambda(x) x)
462                          (list (if (bit-vector-ref flags mxLOGICAL_FLAG) "mxLOGICAL_FLAG" #f)
463                                (if (bit-vector-ref flags mxGLOBAL_FLAG)  "mxGLOBAL_FLAG" #f)
464                                (if (bit-vector-ref flags mxCOMPLEX_FLAG) "mxCOMPLEX_FLAG" #f)))))
465      (fprintf out "#(MAT5:array-flags flags=~A  class=~S  nzmax=~S)"
466               (string-intersperse sflags "|") class nzmax))))
467
468(define-record-printer (MAT5:cell x out)
469  (let ((dims (MAT5:cell-dims x))
470        (data (MAT5:cell-data x)))
471    (fprintf out "#(MAT5:cell dims=~A  data=~S)" dims data)))
472
473;------------------------------------------------------------
474; Predicates and routines for handling MAT-file datatypes and
475; structures
476;
477;
478 
479; Procedure:
480; MAT5:complex-array?:: MAT5:ARRAY-FLAGS -> BOOLEAN
481;
482; Returns true if the complex flag is set in the given array flags
483;
484(define (MAT5:complex-array? flags)
485  (bit-vector-ref (MAT5:array-flags-flags flags)  mxCOMPLEX_FLAG))
486
487
488; Procedure:
489; MAT5:sparse-class?:: MAT5:ARRAY-FLAGS -> BOOLEAN
490;
491; Returns true if the class field in the given array flags is sparse
492;
493(define (MAT5:sparse-class? flags)
494  (cases MAT5:array-class (MAT5:array-flags-class flags)
495         (mxSPARSE_CLASS () #t)
496         (else              #f)))
497
498
499; Procedure:
500; MAT5:cell-class?:: MAT5:ARRAY-FLAGS -> BOOLEAN
501;
502; Returns true if the class field in the given array flags is cell
503;
504(define (MAT5:cell-class? flags)
505  (cases MAT5:array-class (MAT5:array-flags-class flags)
506         (mxCELL_CLASS () #t)
507         (else              #f)))
508
509
510; Procedure:
511; MAT5:structure-class?:: MAT5:ARRAY-FLAGS -> BOOLEAN
512;
513; Returns true if the class field in the given array flags is structure
514;
515(define (MAT5:structure-class? flags)
516  (cases MAT5:array-class (MAT5:array-flags-class flags)
517         (mxSTRUCT_CLASS () #t)
518         (else              #f)))
519
520
521; Procedure:
522; MAT5:object-class?:: MAT5:ARRAY-FLAGS -> BOOLEAN
523;
524; Returns true if the class field in the given array flags is object
525;
526(define (MAT5:object-class? flags)
527  (cases MAT5:array-class (MAT5:array-flags-class flags)
528         (mxOBJECT_CLASS () #t)
529         (else              #f)))
530
531
532; Procedure:
533; MAT5:numeric-type?:: MAT5:DATA-TYPE -> BOOLEAN
534;
535; Return true if type is an atomic numeric type
536;
537(define (MAT5:numeric-type? t)
538  (cases MAT5:data-type t
539         (miINT8   ()     #t)
540         (miUINT8  ()     #t)
541         (miINT16  ()     #t)
542         (miUINT16 ()     #t)
543         (miINT32  ()     #t)
544         (miUINT32 ()     #t)
545         (miSINGLE ()     #t)
546         (miDOUBLE ()     #t)
547         (miINT64  ()     #t)
548         (miUINT64 ()     #t)
549         (miUTF8   ()     #t)
550         (miUTF16  ()     #t)
551         (miUTF32  ()     #t)
552         (else            #f)))
553
554
555; Procedure:
556; MAT5:array-type?:: MAT5:DATA-TYPE -> BOOLEAN
557;
558; Return true if type is array type (miMATRIX)
559;
560(define (MAT5:array-type? t)
561  (cases MAT5:data-type t
562         (miMATRIX  ()     #t)
563         (else            #f)))
564
565
566; Procedure:
567; MAT5:sizeof:: MAT5:ARRAY-FLAGS -> UINTEGER
568;
569; Returns the size in bytes of the atomic datatype t
570;
571; t must not be miMATRIX or miCOMPRESSED
572;
573(define (MAT5:sizeof t)
574  (cases MAT5:data-type t
575         (miINT8   ()     1)
576         (miUINT8  ()     1)
577         (miINT16  ()     2)
578         (miUINT16 ()     2)
579         (miINT32  ()     4)
580         (miUINT32 ()     4)
581         (miSINGLE ()     4)
582         (miDOUBLE ()     8)
583         (miINT64  ()     8)
584         (miUINT64 ()     8)
585         (miUTF8   ()     1)
586         (miUTF16  ()     2)
587         (miUTF32  ()     4)
588         (else            #f)))
589
590
591; Procedure:
592; MAT5:data-type->string:: MAT5:DATA-TYPE -> STRING
593;
594; Returns the string identifier for the given datatype
595;
596(define (MAT5:data-type->string t)
597  (cases MAT5:data-type t
598         (miINT8   ()     "miINT8")
599         (miUINT8  ()     "miUINT8")
600         (miINT16  ()     "miINT16")
601         (miUINT16 ()     "miUINT16")
602         (miINT32  ()     "miINT32")
603         (miUINT32 ()     "miUINT32")
604         (miSINGLE ()     "miSINGLE")
605         ;;      RESERVED     #x00000008
606         (miDOUBLE ()     "miDOUBLE")
607         ;;      RESERVED     #x0000000A
608         ;;      RESERVED     #x0000000B
609         (miINT64  ()     "miINT64")
610         (miUINT64 ()     "miUINT64")
611         (miMATRIX ()     "miMATRIX")
612         (miCOMPRESSED () "miCOMPRESSED")
613         (miUTF8   ()     "miUTF8")
614         (miUTF16  ()     "miUTF16")
615         (miUTF32  ()     "miUTF32")
616         (else            #f)))
617
618
619; Procedure:
620; MAT5:data-type->array-class:: MAT5:DATA-TYPE -> MAT5:ARRAY-FLAGS
621;
622; Return the numeric array class corresponding to the given datatype
623;
624(define (MAT5:data-type->array-class t)
625  (cases MAT5:data-type t
626         (miINT8   ()     (mxINT8_CLASS))
627         (miUINT8  ()     (mxUINT8_CLASS))
628         (miINT16  ()     (mxINT16_CLASS))
629         (miUINT16 ()     (mxUINT16_CLASS))
630         (miINT32  ()     (mxINT32_CLASS))
631         (miUINT32 ()     (mxUINT32_CLASS))
632         (miSINGLE ()     (mxSINGLE_CLASS))
633         ;;      RESERVED     #x00000008
634         (miDOUBLE ()     (mxDOUBLE_CLASS))
635         ;;      RESERVED     #x0000000A
636         ;;      RESERVED     #x0000000B
637         (else            #f)))
638
639; Procedure:
640; MAT5:data-type->word:: MAT5:DATA-TYPE -> UINTEGER
641;
642; Return the numeric code for the given datatype
643;
644(define (MAT5:data-type->word t)
645  (cases MAT5:data-type t
646         (miINT8   ()     #x00000001)
647         (miUINT8  ()     #x00000002)
648         (miINT16  ()     #x00000003)
649         (miUINT16 ()     #x00000004)
650         (miINT32  ()     #x00000005)
651         (miUINT32 ()     #x00000006)
652         (miSINGLE ()     #x00000007)
653         ;;      RESERVED     #x00000008
654         (miDOUBLE ()     #x00000009)
655         ;;      RESERVED     #x0000000A
656         ;;      RESERVED     #x0000000B
657         (miINT64  ()     #x0000000C)
658         (miUINT64 ()     #x0000000D)
659         (miMATRIX ()     #x0000000E)
660         (miCOMPRESSED () #x0000000F)
661         (miUTF8   ()     #x00000010)
662         (miUTF16  ()     #x00000011)
663         (miUTF32  ()     #x00000012)
664         (else            #f)))
665
666
667; Procedure:
668; MAT5:word->data-type:: UINTEGER -> MAT5:DATA-TYPE
669;
670; Return the datatype corresponding to the given numeric value, or #f
671;
672(define (MAT5:word->data-type x)
673  (cond ((= x #x00000001)   (miINT8))
674        ((= x #x00000002)   (miUINT8))
675        ((= x #x00000003)   (miINT16))
676        ((= x #x00000004)   (miUINT16))
677        ((= x #x00000005)   (miINT32))
678        ((= x #x00000006)   (miUINT32))
679        ((= x #x00000007)   (miSINGLE))
680        ((= x #x00000009)   (miDOUBLE))
681        ((= x #x0000000C)   (miINT64))
682        ((= x #x0000000D)   (miUINT64))
683        ((= x #x0000000E)   (miMATRIX))
684        ((= x #x0000000F)   (miCOMPRESSED))
685        ((= x #x00000010)   (miUTF8))
686        ((= x #x00000011)   (miUTF16))
687        ((= x #x00000012)   (miUTF32))
688        (else               #f)))
689
690         
691; Procedure:
692; MAT5:array-class->string:: MAT5:ARRAY-CLASS -> STRING
693;
694; Return the numeric code corresponding to the given array class
695;
696(define (MAT5:array-class->string x)
697  (cases MAT5:array-class x
698         (mxCELL_CLASS   ()    "mxCELL_CLASS")
699         (mxSTRUCT_CLASS ()    "mxSTRUCT_CLASS")
700         (mxOBJECT_CLASS ()    "mxOBJECT_CLASS")
701         (mxCHAR_CLASS   ()    "mxCHAR_CLASS")
702         (mxSPARSE_CLASS ()    "mxSPARSE_CLASS")
703         (mxDOUBLE_CLASS ()    "mxDOUBLE_CLASS")
704         (mxSINGLE_CLASS ()    "mxSINGLE_CLASS")
705         (mxINT8_CLASS   ()    "mxINT8_CLASS")
706         (mxUINT8_CLASS  ()    "mxUINT8_CLASS")
707         (mxINT16_CLASS  ()    "mxINT16_CLASS")
708         (mxUINT16_CLASS ()    "mxUINT16_CLASS")
709         (mxINT32_CLASS  ()    "mxINT32_CLASS")
710         (mxUINT32_CLASS ()    "mxUINT32_CLASS")
711         (else                 #f)))
712
713         
714; Procedure:
715; MAT5:array-class->wordMAT5:ARRAY-CLASS -> UINTEGER
716;
717; Returns the numeric code corresponding to the given array class
718;
719(define (MAT5:array-class->word x)
720  (cases MAT5:array-class x
721         (mxCELL_CLASS   ()    #x01)
722         (mxSTRUCT_CLASS ()    #x02)
723         (mxOBJECT_CLASS ()    #x03)
724         (mxCHAR_CLASS   ()    #x04)
725         (mxSPARSE_CLASS ()    #x05)
726         (mxDOUBLE_CLASS ()    #x06)
727         (mxSINGLE_CLASS ()    #x07)
728         (mxINT8_CLASS   ()    #x08)
729         (mxUINT8_CLASS  ()    #x09)
730         (mxINT16_CLASS  ()    #x0A)
731         (mxUINT16_CLASS ()    #x0B)
732         (mxINT32_CLASS  ()    #x0C)
733         (mxUINT32_CLASS ()    #x0D)
734         (else                 #f)))
735
736
737; Procedure:
738; MAT5:word->array-class:: UINTEGER -> MAT5:ARRAY-CLASS
739;
740;
741; Return the array class corresponding to the given numeric code, or
742; #f
743;
744(define (MAT5:word->array-class x)
745  (cond  ((= x #x01)   (mxCELL_CLASS))
746         ((= x #x02)   (mxSTRUCT_CLASS))
747         ((= x #x03)   (mxOBJECT_CLASS))
748         ((= x #x04)   (mxCHAR_CLASS))
749         ((= x #x05)   (mxSPARSE_CLASS))
750         ((= x #x06)   (mxDOUBLE_CLASS))
751         ((= x #x07)   (mxSINGLE_CLASS))
752         ((= x #x08)   (mxINT8_CLASS))
753         ((= x #x09)   (mxUINT8_CLASS))
754         ((= x #x0A)   (mxINT16_CLASS))
755         ((= x #x0B)   (mxUINT16_CLASS))
756         ((= x #x0C)   (mxINT32_CLASS))
757         ((= x #x0D)   (mxUINT32_CLASS))
758         (else         #f)))
759
760
761; Procedure:
762; MAT5:array-vector-ops:: MAT5:DATA-TYPE -> PROTO * VECTOR-SET * VECTOR-REF * VECTOR-LEN * VECTOR?
763;
764; Returns a set of routines used for the creation and manipulation of
765; homogenous numeric vectors and arrays. Used by the read-array-
766; functions below.
767;
768(define (MAT5:array-vector-ops data-type)
769  (cases MAT5:data-type data-type
770         (miINT8   ()       (values as8  s8vector-set!  s8vector-ref  s8vector-length  s8vector?))
771         (miUINT8  ()       (values au8  u8vector-set!  u8vector-ref  u8vector-length  u8vector?))
772         (miINT16  ()       (values as16 s16vector-set! s16vector-ref s16vector-length s16vector?))
773         (miUINT16 ()       (values au16 u16vector-set! u16vector-ref u16vector-length u16vector?))
774         (miINT32  ()       (values as32 s32vector-set! s32vector-ref s32vector-length s32vector?))
775         (miUINT32 ()       (values au32 u32vector-set! u32vector-ref u32vector-length u32vector?))
776;;       (miINT64  ()       (values as64 s64vector-set! s64vector-ref s64vector-length))
777;;       (miUINT64 ()       (values au64 u64vector-set! u64vector-ref u64vector-length))
778         (miUTF8   ()       (values au8   u8vector-ref   u8vector-ref  u8vector-length  u8vector?))
779         (miUTF16  ()       (values au16  u16vector-set! u16vector-ref u16vector-length u16vector?))
780         (miUTF32  ()       (values au32  u32vector-set! u32vector-ref u32vector-length u32vector?))
781         (miSINGLE ()       (values ar32  f32vector-set! f32vector-ref f32vector-length f32vector?))
782         (miDOUBLE ()       (values ar64  f64vector-set! f64vector-ref f64vector-length f64vector?))
783         (miMATRIX ()       (MAT5:error "nested arrays not permitted in numeric arrays"))
784         (else              (MAT5:error "unrecognized type " data-type))))
785
786; Procedure:
787; MAT5:array-vector-length:: MAT5:DATA-TYPE * UINTEGER -> UINTEGER
788;
789; Returns the length of the vector that would be necessary to hold all
790; values for an array of the specified size and type.
791;
792;   * data-type   array/vector type
793;   * data-size   array/vector size in bytes
794;
795(define (MAT5:array-vector-length data-type data-size)
796  (cases MAT5:data-type data-type
797         (miINT8   ()        data-size)
798         (miUINT8  ()        data-size)
799         (miINT16  ()        (/ data-size (MAT5:sizeof data-type)))
800         (miUINT16 ()        (/ data-size (MAT5:sizeof data-type)))
801         (miINT32  ()        (/ data-size (MAT5:sizeof data-type)))
802         (miUINT32 ()        (/ data-size (MAT5:sizeof data-type)))
803;;       (miINT64  ()        (/ data-size (MAT5:sizeof data-type)))
804;;       (miUINT64 ()        (/ data-size (MAT5:sizeof data-type)))
805         (miUTF8   ()        data-size)
806         (miUTF16  ()        (/ data-size (MAT5:sizeof data-type)))
807         (miUTF32  ()        (/ data-size (MAT5:sizeof data-type)))
808         (miSINGLE ()        (/ data-size (MAT5:sizeof data-type)))
809         (miDOUBLE ()        (/ data-size (MAT5:sizeof data-type)))
810         (miMATRIX ()       (MAT5:error "nested arrays not permitted in numeric arrays"))
811         (else              (MAT5:error "unrecognized type " data-type))))
812 
813
814; Procedure:
815; MAT5:array-vector-make:: MAT5:DATA-TYPE * UINTEGER -> VECTOR
816;
817; Returns a vector used to hold the values for an array of the
818; specified size and type.  Used by the read-array- functions below.
819;
820;   * data-type   array/vector type
821;   * data-size   array/vector size in bytes
822;
823(define (MAT5:array-vector-make data-type data-size)
824  (let ((len (MAT5:array-vector-length data-type data-size)))
825    (cases MAT5:data-type data-type
826         (miINT8   ()        (make-s8vector   len))
827         (miUINT8  ()        (make-u8vector   len))
828         (miINT16  ()        (make-s16vector  len))
829         (miUINT16 ()        (make-u16vector  len))
830         (miINT32  ()        (make-s32vector  len))
831         (miUINT32 ()        (make-u32vector  len))
832;;       (miINT64  ()        (make-s64vector  len))
833;;       (miUINT64 ()        (make-u64vector  len))
834         (miUTF8   ()        (make-u8vector   len))
835         (miUTF16  ()        (make-u16vector  len))
836         (miUTF32  ()        (make-u32vector  len))
837         (miSINGLE ()        (make-f32vector  len))
838         (miDOUBLE ()        (make-f64vector  len))
839         (miMATRIX ()       (MAT5:error "nested arrays not permitted in numeric arrays"))
840         (else              (MAT5:error "unrecognized type " data-type)))))
841 
842
843; Procedure:
844; MAT5:vector-make:: MAT5:DATA-TYPE * UINTEGER -> VECTOR
845;
846; Returns a vector used to hold the values of the specified type.
847;
848;   * data-type   vector type
849;   * len         vector length
850;
851(define (MAT5:vector-make data-type len)
852    (cases MAT5:data-type data-type
853         (miINT8   ()        (make-s8vector   len))
854         (miUINT8  ()        (make-u8vector   len))
855         (miINT16  ()        (make-s16vector  len))
856         (miUINT16 ()        (make-u16vector  len))
857         (miINT32  ()        (make-s32vector  len))
858         (miUINT32 ()        (make-u32vector  len))
859;;       (miINT64  ()        (make-s64vector  len))
860;;       (miUINT64 ()        (make-u64vector  len))
861         (miUTF8   ()        (make-u8vector   len))
862         (miUTF16  ()        (make-u16vector  len))
863         (miUTF32  ()        (make-u32vector  len))
864         (miSINGLE ()        (make-f32vector  len))
865         (miDOUBLE ()        (make-f64vector  len))
866         (miMATRIX ()       (MAT5:error "nested arrays not permitted in numeric arrays"))
867         (else              (MAT5:error "unrecognized type " data-type))))
868 
869
870;--------------------
871;  MAT5 cell routines
872;
873;  A MAT5 cell is represented as an R5RS non-homogenous array that
874;  contains elements of type MAT5:array
875;
876
877;  Structure: MAT5:cell
878;
879;  A representation of MAT5 cells.
880;
881;  * dims is the cell dimensions (see MAT5:dimensions?) 
882;  * data is an R5RS non-homogenous array that contains
883;    elements of type MAT5:array.
884(define-record MAT5:cell  dims  data)
885
886; Procedure:
887; init-MAT5:cell:: UINTEGER * ... -> MAT5:CELL
888;
889; Create a new MAT5:cell object with the specified dimensions. The
890; dimensions must all be positive integers, and at least one dimension
891; must be specified.
892;
893(define (init-MAT5:cell d . dimensions)
894  (let ((dims (cons d dimensions)))
895    (if (not (natural-list? dims))   (MAT5:error "invalid cell dimensions")
896        (let ((ar  (make-vector (car dims))))
897          (init-cell1 ar (cdr dims))
898          (make-MAT5:cell dims ar)))))
899
900; helper function for init-MAT5:cell above
901(define (init-cell1 v dims)
902  (if (not (null? dims))
903      (let loop ((i 0) (len (vector-length v)))
904        (if (not (zero? len))
905            (let ((v1  (make-vector (car dims))))
906              (vector-set! v i v1)
907              (init-cell1 v1 (cdr dims))
908              (loop (+ i 1) (- len 1)))))))
909
910; Procedure:
911; MAT5:cell-ref:: MAT5:CELL * UINTEGER ... -> VALUE
912;
913; Given a MAT5:cell object and an index, returns the value found at
914; that index in the cell. The index must be a list of positive
915; integers, and it must be within the bounds of the cell dimensions.
916;
917(define (MAT5:cell-ref x i . rest)
918  (let ((ar  (MAT5:cell-data x))
919        (il  (cons i rest)))
920    (if (not (positive-or-zero-list? il))
921        (MAT5:error "invalid cell dimensions")
922        (MAT5:cell-ref1 ar il))))
923 
924;; helper function for MAT5:cell-ref above
925;;  i is a list of dimensions
926(define (MAT5:cell-ref1 x i)
927  (cond ((vector? x)  (if (null? i) (MAT5:error "cell dimension mismatch")
928                          (MAT5:cell-ref1 (vector-ref x (car i)) (cdr i))))
929        (else    (if (null? i) x   (MAT5:error "cell dimension mismatch")))))
930
931; Procedure:
932; MAT5:cell-set!:: MAT5:CELL * VALUE * UINTEGER ... -> UNDEFINED
933;
934; Given a MAT5:cell object and an index, destructively replaces the
935; element at the given index of the cell with the given value.
936;
937(define (MAT5:cell-set! x  v  i . rest)
938  (let ((ar  (MAT5:cell-data x))
939        (il  (cons i rest)))
940    (if (not (positive-or-zero-list? il))
941        (MAT5:error "invalid cell dimensions")
942        (MAT5:cell-set1 ar v il))))
943 
944
945; helper function for MAT5:cell-set! above
946(define (MAT5:cell-set1 x v i)
947  (cond ((null? i)        (MAT5:error "invalid cell dimensions"))
948        ((null? (cdr i))  (if (and (vector? x) (not (vector? (vector-ref x (car i)))))
949                              (vector-set! x (car i) v)
950                              (MAT5:error "cell dimension mismatch")))
951        (else             (if (vector? x)
952                              (MAT5:cell-set1 (vector-ref x (car i)) v (cdr i))
953                              (MAT5:error "cell dimension mismatch")))))
954
955;---------------------------------------
956;  MAT5 vector/array conversion routines
957;
958
959
960; Procedure:
961; vector->MAT5:cell:: VECTOR * UINTEGER LIST [* ORDER] -> MAT5:CELL
962;
963; Given a vector, creates a new MAT5:cell object that consists of the
964; elements of the vector. The vector must be of length equal to the
965; total size of the array, and its elements are used to initialize the
966; cell in either row-major order (left to right and top to bottom), or
967; in column-major order (top to bottom and then left to right).
968;
969; The optional argument ORDER specifies the initialization order and
970; can be either 'row-major or 'col-major. The default is 'row-major.
971;
972; This is based on vector->array from SRFI-63 reference implementation
973; Copyright (C) 2001, 2003, 2005 Aubrey Jaffer
974(define (vector->MAT5:cell vect dimensions . rest)
975  (let-optionals  rest ((order 'row-major))
976   (let* ((vdx (vector-length vect)))
977     (if (not (eqv? vdx (apply * dimensions))) (MAT5:error "incompatible dimensions"))
978     (letrec ((ra (apply init-MAT5:cell  dimensions))
979              (v2ra  (lambda (dims idxs)
980                       (cond ((null? dims)
981                              (set! vdx (+ -1 vdx))
982                              (apply MAT5:cell-set! ra (vector-ref vect vdx) (reverse idxs)))
983                             (else
984                              (do ((idx (+ -1 (car dims)) (+ -1 idx)))
985                                  ((negative? idx) vect)
986                                (v2ra (cdr dims) (cons idx idxs)))))))
987              (v2ca   (lambda (dims idxs)
988                        (cond ((null? dims)
989                               (set! vdx (+ -1 vdx))
990                               (apply MAT5:cell-set! ra (vector-ref vect vdx) idxs))
991                              (else
992                               (do ((idx (+ -1 (car dims)) (+ -1 idx)))
993                                   ((negative? idx) vect)
994                                 (v2ca (cdr dims) (cons idx idxs))))))))
995       (if (eq? order 'row-major)
996                (v2ra dimensions '())
997                (v2ca (reverse dimensions) '()))
998       ra))))
999
1000
1001; Procedure:
1002; MAT5:array-foldi:: (INDEX * VALUE * AX -> AX) [* ORDER] -> MAT5:ARRAY * AX -> AX
1003;
1004; Iterator function for non-homogenous MAT5:array objects; that is,
1005; objects that are either MAT5:cell-array (MAT5 cell) or
1006; MAT5:structure (MAT5 structure).
1007;
1008; Analogous to the list iterator, this procedure repeatedly applies
1009; the given function to each element of the array, and accumulates the
1010; return value. The order of iteration is specified by the optional
1011; argument ORDER, which can be 'row-major (left to right and top to
1012; bottom) or 'col-major (top to bottom and then left to right). The
1013; default is 'row-major.
1014;
1015(define (MAT5:array-foldi  f . rest)
1016  (let-optionals  rest ((order 'row-major))
1017   (lambda (x ax)
1018     (let-values (((dims vdx arr elm-ref) 
1019                   (cases MAT5:array x
1020                          (MAT5:cell-array (name dims cell) 
1021                                           (let ((vdx (apply * dims)))
1022                                             (values dims vdx cell MAT5:cell-ref)))
1023                          (MAT5:structure (name dims field-names fields) 
1024                                          (let* ((flen  (length field-names))
1025                                                 (vdx   (* (apply * dims) flen))) 
1026                                            (values dims vdx fields MAT5:cell-ref)))
1027                          (MAT5:object (name dims class-name field-names fields)
1028                                       (let* ((flen  (length field-names))
1029                                              (vdx   (* (apply * dims) flen))) 
1030                                         (values dims vdx fields MAT5:cell-ref)))
1031                          (else   (MAT5:error "invalid MAT5 array")))))
1032                 (letrec ((ra2v  (lambda (dims idxs ax)
1033                                   (cond ((null? dims)
1034                                          (set! vdx (+ -1 vdx))
1035                                          (f vdx  (apply elm-ref arr (reverse idxs)) ax))
1036                                         (else
1037                                          (let ((dim (car dims)))
1038                                            (do ((idx 0 (+ 1 idx)))
1039                                                ((>= idx dim) ax)
1040                                              (set! ax (ra2v (cdr dims) (cons idx idxs) ax))))))))
1041                          (ca2v   (lambda (dims idxs ax)
1042                                    (cond ((null? dims)
1043                                           (set! vdx (+ -1 vdx))
1044                                           (f vdx  (apply elm-ref arr  idxs) ax))
1045                                          (else
1046                                           (let ((dim (car dims)))
1047                                             (do ((idx 0 (+ 1 idx)))
1048                                                 ((>= idx dim) ax)
1049                                             (set! ax (ca2v (cdr dims) (cons idx idxs) ax)))))))))
1050                   (if (eq? order 'row-major)
1051                       (ra2v dims '() ax)
1052                       (ca2v (reverse dims) '() ax)))))))
1053
1054
1055; Procedure:
1056; array-foldi
1057;
1058(define (array-foldi  f . rest)
1059  (let-optionals  rest ((order 'row-major))
1060   (lambda (arr ax)
1061     (if (not (array? arr)) (MAT5:error "invalid array"))
1062     (let ((dims (array-dimensions arr)))
1063       (letrec ((vdx   (apply * dims))
1064                (ra2v  (lambda (dims idxs ax)
1065                         (cond ((null? dims)
1066                                (set! vdx (+ -1 vdx))
1067                                (f vdx  (apply array-ref arr (reverse idxs)) ax))
1068                               (else
1069                                (let ((dim (car dims)))
1070                                  (do ((idx 0 (+ 1 idx)))
1071                                      ((>= idx dim) ax)
1072                                    (set! ax (ra2v (cdr dims) (cons idx idxs) ax))))))))
1073                (ca2v   (lambda (dims idxs ax)
1074                          (cond ((null? dims)
1075                                 (set! vdx (+ -1 vdx))
1076                                 (f vdx  (apply array-ref arr idxs) ax))
1077                                (else
1078                                 (let ((dim (car dims)))
1079                                   (do ((idx 0 (+ 1 idx)))
1080                                       ((>= idx dim) ax)
1081                                     (set! ax (ca2v (cdr dims) (cons idx idxs) ax)))))))))
1082         (if (eq? order 'row-major)
1083             (ra2v dims '() ax)
1084             (ca2v (reverse dims) '() ax)))))))
1085
1086
1087
1088;--------------------------
1089; MAT-file Reading Routines
1090;
1091
1092; Procedure:
1093; MAT5:read-header:: ENDIAN-PORT -> MAT5:HEADER
1094;
1095; Reads a MAT5 header from the given endian port. Returns a
1096; MAT5:header record.
1097;
1098(define (MAT5:read-header eport)
1099  (endian-port-setpos eport 0)
1100  (let* ((text     (endian-port-read-byte-vector  eport MAT5:header-text-size (MSB)))
1101         (subsys   (endian-port-read-byte-vector  eport MAT5:header-subsys-size (MSB)))
1102         (version  (endian-port-read-int2 eport (MSB)))
1103         (magic    (endian-port-read-int2 eport (MSB))))
1104    (cond ((= magic MAT5:lsb-magic) (endian-port-set-littlendian! eport))
1105          ((= magic MAT5:msb-magic) (endian-port-set-bigendian! eport))
1106          (else                     (MAT5:error "MAT-file magic number not found")))
1107    (if (not (or (= version MAT5:msb-version) (= version MAT5:lsb-version)))
1108        (MAT5:warning "unknown MAT-file version " version))
1109    (make-MAT5:header magic (byte-vector->zstring text) (byte-vector->zstring subsys) version eport)))
1110
1111
1112
1113; Procedure:
1114; read-data-element-header:: ENDIAN-PORT -> MAT5:DATA-TYPE * DATA-SIZE * BYTES
1115;
1116; Reads the header of a MAT5 data element.
1117;
1118; Returns (values data-type data-size bytes)
1119;
1120(define (read-data-element-header eport)
1121  (let ((type-word  (endian-port-read-int4 eport)))
1122    ;; Check for small data element format: when reading a MAT-file,
1123    ;; determine that we are processing a small data element by
1124    ;; comparing the value of the first two bytes of the tag with the
1125    ;; value zero. If these two bytes are not zero, the tag uses the
1126    ;; small data element format.
1127    (if (zero? (bitwise-and #xFFFF0000 type-word))
1128        (let ((data-type (MAT5:word->data-type type-word)))
1129          (values   data-type
1130                    (endian-port-read-int4 eport)
1131                    8))
1132        (let ((data-type (MAT5:word->data-type (bitwise-and #x0000FFFF type-word))))
1133          (values  data-type
1134                   (arithmetic-shift (bitwise-and #xFFFF0000 type-word) -16)
1135                   4)))))
1136
1137
1138; Procedure:
1139; read-data-element:: ENDIAN-PORT * DATA-TYPE * SIZE -> MAT5:DATA-ELEMENT
1140;
1141; Given an eport and MAT5 data type, read a word of that type from
1142; the eport. Returns a record of type MAT5:data-element.
1143;
1144; This function is parameterized over the routines for reading numeric
1145; array data, which can be either strict (read all data at once) or
1146; stream-based (read the data lazily column vector by column
1147; vector). See function read-array and the stream routines further
1148; down for understanding of the two interfaces.
1149;
1150(define (read-data-element read-sparse-array-data 
1151                           read-num-array-data)
1152  (let ((read-array (read-array read-sparse-array-data 
1153                                read-num-array-data)))
1154    (lambda (eport type size)
1155      (if (MAT5:numeric-type? type) 
1156          (read-num-data-element eport type size)
1157          (cases MAT5:data-type type
1158                 (miMATRIX ()     (let-values (((array bytes)  (read-array eport)))
1159                                              (make-MAT5:data-element type bytes array)))
1160                 (miCOMPRESSED () (read-compressed-data-element eport size))
1161                 
1162                 (else (MAT5:error eport "unrecognized type " type)))))))
1163
1164
1165; Given a compressed data element, uncompress it, create a temporary
1166; eport that points to the uncompressed data, then call
1167; MAT5:read-data-element with that eport
1168(define (read-compressed-data-element eport size)
1169  (let* ((h      (z3:decode-init))
1170         (zdata  (endian-port-read-byte-vector  eport size (MSB)))
1171         (dest   (open-output-string))
1172         (recv   (lambda (x) (write x dest))))
1173    (let loop ((z zdata))
1174      (let ((bytes  (z3:decode h recv z)))
1175        (if bytes 
1176            (loop (substring z bytes)))))
1177    (let* ((zport (open-input-string (get-output-string dest)))
1178           (ezport (port->endian-port zport)))
1179      (if (eq? (endian-port-byte-order eport) (MSB))
1180          (endian-port-set-bigendian! ezport)
1181          (endian-port-set-littlendian! ezport))
1182      (MAT5:read-data-element ezport))))
1183
1184; Procedure:
1185; read-num-data-element:: ENDIAN-PORT * MAT5:DATA-TYPE * SIZE -> MAT5:DATA-ELEMENT
1186;
1187; Reads an atomic numeric data element (i.e. one that is not a matrix).
1188;
1189(define (read-num-data-element eport type size)
1190  (cases MAT5:data-type type
1191         (miINT8   ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int1 eport)))
1192         (miUINT8  ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int1 eport)))
1193         (miINT16  ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int2 eport)))
1194         (miUINT16 ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int2 eport)))
1195         (miINT32  ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int4 eport)))
1196         (miUINT32 ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int4 eport)))
1197         (miSINGLE ()     (let-values (((nan fpword) (endian-port-read-ieee-float32 eport)))
1198                                      (make-MAT5:data-element type (MAT5:sizeof type)  (if nan +nan  fpword))))
1199         (miDOUBLE ()     (let-values (((nan fpword) (endian-port-read-ieee-float64 eport)))
1200                                      (make-MAT5:data-element type (MAT5:sizeof type)  (if nan +nan  fpword))))
1201         ;;      (miINT64  ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int8 eport)))
1202         ;;      (miUINT64 ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int8 eport)))
1203         (miUTF8   ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int1 eport)))
1204         (miUTF16  ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int2 eport)))
1205         (miUTF32  ()     (make-MAT5:data-element type (MAT5:sizeof type)  (endian-port-read-int4 eport)))
1206         (else            (MAT5:error eport "unrecognized type " type))))
1207         
1208
1209
1210; Procedure:
1211; read-array-flags:: ENDIAN-PORT -> MAT5:ARRAY-FLAGS * BYTES
1212;
1213; Reads array flags data element from the given eport, and returns the
1214; array flags and how many bytes were read.
1215;
1216(define (read-array-flags eport)
1217  (let-values
1218   (((data-type data-size header-bytes) (read-data-element-header eport)))
1219   (if (not data-type) (MAT5:error eport "invalid data element type"))
1220   (cases MAT5:data-type data-type
1221          (miUINT32 ()  (let* ((flags-word  (endian-port-read-int4 eport))
1222                               (flags       (integer->bit-vector (arithmetic-shift 
1223                                                                  (bitwise-and #x0000FF00 flags-word) -8)))
1224                               (class-word  (bitwise-and #x000000FF flags-word))
1225                               (class       (MAT5:word->array-class class-word))
1226                               (nzmax       (endian-port-read-int4 eport))
1227                               (pad-bytes   (align-eport-pos eport)))
1228                          (if (not class) (MAT5:error eport "invalid class: " class-word))
1229                          (values (make-MAT5:array-flags flags class nzmax) (+ header-bytes pad-bytes data-size))))
1230            (else         (MAT5:error eport "array flags data element is not of type UINT32; type is "
1231                                      data-type)))))
1232         
1233
1234; Procedure:
1235; read-array-dimensions:: ENDIAN-PORT -> UINTEGER LIST * BYTES
1236;
1237; Reads array dimensions data element from the given eport, and
1238; returns the dimensions as a list of positive integers, and how many
1239; bytes were read.
1240;
1241(define (read-array-dimensions eport)
1242  (let-values
1243   (((data-type data-size header-bytes) (read-data-element-header eport)))
1244   (if (not data-type) (MAT5:error eport "invalid data element type"))
1245   (let* ((count      (/ data-size (MAT5:sizeof data-type))))
1246     (cases MAT5:data-type data-type
1247            (miINT32 ()  (if (> count 1)
1248                             (let loop ((i 1) (lst (list (endian-port-read-int4 eport))))
1249                               (if (< i count)
1250                                   (loop (+ i 1) (cons (endian-port-read-int4 eport) lst))
1251                                   (let ((pad-bytes (align-eport-pos eport)))
1252                                     (values (reverse lst) (+ header-bytes pad-bytes data-size)))))
1253                             (MAT5:error eport count " array dimensions found; at least 2 array dimensions required")))
1254            (else         (MAT5:error eport "array dimension data element is not of type INT32"))))))
1255
1256       
1257; Procedure:
1258; read-array-name:: ENDIAN-PORT -> STRING * BYTES
1259;
1260;
1261; Reads array name data element from the given eport, and returns the
1262; name as a string, and how many bytes were read.
1263;
1264(define (read-array-name eport)
1265  (let-values
1266   (((data-type data-size header-bytes) (read-data-element-header eport)))
1267   (if (not data-type) (MAT5:error eport "invalid data element type"))
1268   (cases MAT5:data-type data-type
1269          (miINT8 ()   (let* ((bv  (endian-port-read-byte-vector eport data-size (MSB)))
1270                              (pad-bytes  (align-eport-pos eport)))
1271                         (values (byte-vector->zstring bv) (+ header-bytes pad-bytes data-size))))
1272          (else        (MAT5:error eport "array name data element is not of type INT8")))))
1273
1274
1275; Procedure:
1276; read-row-index:: ENDIAN-PORT * MAT5:ARRAY-FLAGS -> UINTEGER LIST * BYTES
1277;
1278; Reads a sparse array row-index data element from the given
1279; eport. Argument array-flags is used to determine the number of
1280; non-zero rows.
1281;
1282(define (read-row-index eport array-flags)
1283  (let-values
1284   (((data-type data-size header-bytes) (read-data-element-header eport)))
1285   (if (not data-type) (MAT5:error eport "invalid data element type"))
1286   (let*  ((count      (/ data-size (MAT5:sizeof data-type)))
1287           (nzmax      (MAT5:array-flags-nzmax array-flags)))
1288     (cases MAT5:data-type data-type
1289            (miINT32 ()   (if (= nzmax count)
1290                              (let loop ((i 1) (lst (list (endian-port-read-int4 eport))))
1291                                (if (< i count)
1292                                    (loop (+ i 1) (cons (endian-port-read-int4 eport) lst))
1293                                    (let ((pad-bytes  (align-eport-pos eport)))
1294                                      (values (reverse lst) (+ header-bytes pad-bytes data-size)))))
1295                              (MAT5:error eport "mismatch between ir count and nzmax: ir count = " 
1296                                          count " nzmax = " nzmax)))
1297            (else          (MAT5:error eport "array row index data element is not of type INT32"))))))
1298
1299
1300; Procedure:
1301; read-col-index:: ENDIAN-PORT * MAT5:ARRAY-FLAGS -> UINTEGER LIST * BYTES
1302;
1303; Reads a sparse array column-index data element from the given
1304; eport. Argument array-flags is used to determine the number of
1305; non-zero rows.
1306;
1307(define (read-col-index eport array-dims)
1308  (let-values
1309   (((data-type data-size header-bytes) (read-data-element-header eport)))
1310   (if (not data-type) (MAT5:error eport "invalid data element type"))
1311   (let* ((count      (/ data-size (MAT5:sizeof data-type))))
1312     (cases MAT5:data-type data-type
1313            (miINT32 ()   (if (= 1 (- count (second array-dims)))
1314                              (let loop ((i 1) (lst (list (endian-port-read-int4 eport))))
1315                                (if (< i count)
1316                                    (loop (+ i 1) (cons (endian-port-read-int4 eport) lst))
1317                                    (let ((pad-bytes  (align-eport-pos eport)))
1318                                      (values (reverse lst) (+ header-bytes pad-bytes data-size)))))
1319                              (MAT5:error eport "mismatch between jc count and second dimension: jc count = " count
1320                                          " dimensions = " array-dims)))
1321            (else (MAT5:error eport "array column index data element is not of type INT32"))))))
1322
1323
1324; Procedure:
1325; read-fieldname-len
1326;
1327;
1328; Read field name length (for structures and objects)
1329;
1330(define (read-fieldname-len eport)
1331  (let-values
1332   (((data-type data-size header-bytes)  (read-data-element-header eport)))
1333   (if (not data-type) (MAT5:error eport "invalid data element type"))
1334   (cases MAT5:data-type data-type
1335          (miINT32 ()     (let* ((fieldname-len    (endian-port-read-int4 eport))
1336                                 (pad-bytes (align-eport-pos eport)))
1337                            (values fieldname-len (+ header-bytes pad-bytes data-size))))
1338          (else           (MAT5:error eport "field name length data element is not of type INT32")))))
1339
1340
1341; Procedure:
1342; read-field-names:: ENDIAN-PORT * UINTEGER -> STRING LIST * BYTES
1343;
1344; Reads field names (for structures and objects). Returns the field
1345; names as a list of strings, and how many bytes were read.
1346;
1347(define (read-field-names eport fieldname-len)
1348  (let-values
1349   (((data-type data-size header-bytes)  (read-data-element-header eport)))
1350   (if (not data-type) (MAT5:error eport "invalid data element type"))
1351   (cases MAT5:data-type data-type
1352            (miINT8 ()      (let loop ((names '())   (len (/ data-size fieldname-len)))
1353                              (if (zero? len)
1354                                  (let ((pad-bytes (align-eport-pos eport))
1355                                        (bv->string (lambda (x) (byte-vector->zstring x))))
1356                                    (values (map bv->string (reverse names))  (+ header-bytes pad-bytes data-size)))
1357                                  (let  ((name  (endian-port-read-byte-vector eport fieldname-len (MSB))))
1358                                    (loop (cons name names) (- len 1))))))
1359            (else           (MAT5:error eport "field name data element is not of type INT8")))))
1360
1361; Procedure:
1362; read-fields:: ENDIAN-PORT * STRING LIST -> VALUE LIST * BYTES
1363;
1364; Given an endian port, and a list of field names, reads the field
1365; values for a MAT5 structure. This function is parameterized over
1366; the routines for reading numeric array data, which can be either
1367; strict (read all data at once) or stream-based (read the data lazily
1368; column vector by column vector). See function read-array and the
1369; stream routines further down for understanding of the two
1370; interfaces.
1371
1372(define (read-fields read-sparse-array-data 
1373                     read-num-array-data)
1374  (lambda (eport field-names)
1375    (let ((read-data-element (read-data-element read-sparse-array-data 
1376                                                read-num-array-data)))
1377      (let loop ((fields '()) (field-names field-names) (bytes 0))
1378        (if (null? field-names)  (values fields bytes)
1379            (let-values
1380             (((data-type data-size header-bytes) (read-data-element-header eport)))
1381             (if (not data-type) (MAT5:error eport "invalid data element type"))
1382             (let* ((data-element   (read-data-element eport data-type data-size))
1383                    (data           (MAT5:data-element-data data-element))
1384                    (data-bytes     (MAT5:data-element-bytes data-element))
1385                    (field          data))
1386;                   (field          `(,(car field-names) . ,data)))
1387               (loop (cons field fields) (cdr field-names) (+ bytes header-bytes data-bytes)))))))))
1388
1389; Procedure:
1390; read-struct-data:: ENDIAN-PORT * MAT5:ARRAY-FLAGS * MAT5:DIMENSIONS * STRING LIST ->
1391;                    MAT5:CELL * BYTES
1392;
1393; Reads the data for a MAT5 structure. The structure is represented
1394; as a MAT5:cell object, where each element is a list of field values.
1395;
1396(define (read-struct-data read-sparse-array-data 
1397                          read-num-array-data)
1398  (let ((read-fields (read-fields read-sparse-array-data 
1399                                  read-num-array-data)))
1400    (lambda (eport array-flags array-dims field-names)
1401      (let ((vector-data (make-vector (apply * array-dims))))
1402        (let loop ((i 0)  (len (apply * array-dims)) (bytes 0))
1403          (if (zero? len)
1404              (let ((pad-bytes  (align-eport-pos eport)))
1405                (values (vector->MAT5:cell vector-data array-dims 'col-major) 
1406                        (+ pad-bytes bytes)))
1407              (let-values (((fields fields-bytes)  (read-fields eport field-names)))
1408                          (vector-set! vector-data i (reverse fields))
1409                          (loop (+ i 1) (- len 1) (+ bytes fields-bytes)))))))))
1410
1411
1412; Procedure:
1413; read-cell-data:: ENDIAN-PORT * MAT5:ARRAY-FLAGS * MAT5:DIMENSIONS ->
1414;                  MAT5:CELL * BYTES
1415;
1416; Reads the data for a MAT5 cell.
1417;
1418(define (read-cell-data read-sparse-array-data 
1419                        read-num-array-data)
1420  (let ((read-data-element (read-data-element read-sparse-array-data 
1421                                              read-num-array-data)))
1422    (lambda (eport array-flags array-dims)
1423      (let ((vector-data (make-vector (apply * array-dims))))
1424        (let loop ((i 0)  (len (apply * array-dims)) (bytes 0))
1425          (if (zero? len)
1426              (let ((pad-bytes  (align-eport-pos eport)))
1427                (values (vector->MAT5:cell vector-data array-dims 'col-major) 
1428                        (+ pad-bytes bytes)))
1429              (let-values
1430               (((data-type data-size header-bytes) (read-data-element-header eport)))
1431               (if (not data-type) (MAT5:error eport "invalid data element type"))
1432               (let ((data-element (read-data-element eport data-type data-size)))
1433                 (vector-set! vector-data i (MAT5:data-element-data data-element))
1434                 (loop (+ i 1)  (- len 1) (+ bytes header-bytes data-size))))))))))
1435
1436
1437; Procedure:
1438; read-sparse-array-data:: ENDIAN-PORT * MAT5:ARRAY-FLAGS * ROW-INDEX * COL-INDEX ->
1439;                          VECTOR (of MAT5:data-element) * BYTES
1440;
1441;
1442; Reads the data for a sparse numeric array.
1443;
1444(define (read-sparse-array-data eport array-flags row-index col-index)
1445  (let-values
1446   (((data-type data-size header-bytes) (read-data-element-header eport)))
1447   (if (not data-type) (MAT5:error eport "invalid data element type"))
1448   (if (not (MAT5:numeric-type? data-type)) (MAT5:error eport "non-numeric sparse array data"))
1449   (if (not (= (MAT5:array-vector-length data-type data-size) 
1450               (* (length row-index) (length col-index))))
1451       (MAT5:error eport "incompatible dimensions: data length is " 
1452                   (MAT5:array-vector-length data-type data-size)
1453                   " sparse array dimensions are specified as " 
1454                   row-index " " col-index))
1455   (let-values   
1456    (((prototype vector-set! vector-ref vector-length vector?) 
1457      (MAT5:array-vector-ops data-type)))
1458    (let ((vector-data (MAT5:array-vector-make data-type data-size)))
1459      (let loop ((i 0) (len (vector-length vector-data)))
1460        (if (zero? len)
1461            (let ((vops     `((vector-ref . ,vector-ref)
1462                              (vector-length . ,vector-length)))
1463                  (pad-bytes  (align-eport-pos eport)))
1464              (values data-type
1465                      (vector->array vops vector-data (prototype) 
1466                                     (list (length row-index) (length col-index)) 'col-major) 
1467                      (+ pad-bytes header-bytes data-size)))
1468            (let ((data-element (read-num-data-element eport data-type data-size)))
1469              (vector-set! vector-data i (MAT5:data-element-data data-element))
1470              (loop (+ i 1)  (- len 1)))))))))
1471
1472; Procedure:
1473; read-num-array-data:: ENDIAN-PORT * MAT5:ARRAY-FLAGS * MAT5:DIMENSIONS ->
1474;                       ARRAY * BYTES
1475;
1476; Reads the data for a homogeneous numeric array.
1477(define (read-num-array-data eport array-flags array-dims)
1478  (let-values
1479   (((data-type data-size header-bytes) (read-data-element-header eport)))
1480   (if (not data-type) (MAT5:error eport "invalid data element type"))
1481   (if (not (MAT5:numeric-type? data-type)) (MAT5:error eport "non-numeric array data"))
1482   (if (not (= (MAT5:array-vector-length data-type data-size) (apply * array-dims)))
1483       (MAT5:error eport "incompatible dimensions: data length is " 
1484                   (MAT5:array-vector-length data-type data-size)
1485                   " array dimensions are specified as " array-dims))
1486   (let-values   
1487    (((prototype vector-set! vector-ref vector-length vector?) 
1488      (MAT5:array-vector-ops data-type)))
1489    (let ((vector-data (MAT5:array-vector-make data-type data-size)))
1490      (let loop ((i 0) (len (vector-length vector-data)))
1491        (if (zero? len)
1492            (let* ((vops      `((vector-ref . ,vector-ref)
1493                                (vector-length . ,vector-length)))
1494                   (pad-bytes  (align-eport-pos eport)))
1495              (values  data-type
1496                       (vector->array vops vector-data (prototype) array-dims 'col-major) 
1497                       (+ pad-bytes header-bytes data-size)))
1498            (let ((data-element (read-num-data-element eport data-type data-size)))
1499              (vector-set! vector-data i (MAT5:data-element-data data-element))
1500              (loop (+ i 1)  (- len 1)))))))))
1501
1502
1503; Procedure:
1504; read-array
1505;
1506(define (read-array read-sparse-array-data 
1507                    read-num-array-data)
1508  (lambda (eport)
1509    (let ((read-struct-data (read-struct-data read-sparse-array-data 
1510                                              read-num-array-data))
1511          (read-cell-data (read-cell-data read-sparse-array-data 
1512                                          read-num-array-data)))
1513         
1514    (let*-values 
1515     (((array-flags flags-bytes)  (read-array-flags eport))
1516      ((array-dims  dims-bytes)   (read-array-dimensions eport))
1517      ((array-name  name-bytes)   (read-array-name eport)))
1518     (cond  
1519     
1520     
1521      ;; read an object
1522      ((MAT5:object-class? array-flags) 
1523       (begin
1524         (let*-values (((class-name      class-name-bytes)      (read-array-name eport))
1525                       ((fieldname-len   fieldname-len-bytes)   (read-fieldname-len eport))
1526                       ((field-names     field-names-bytes)     (read-field-names eport fieldname-len))
1527                       ((fields          fields-bytes)          (read-struct-data eport array-flags array-dims 
1528                                                                                  field-names)))
1529                      (values  (MAT5:object array-name array-dims class-name field-names fields)
1530                               (+ flags-bytes dims-bytes name-bytes class-name-bytes 
1531                                  fieldname-len-bytes field-names-bytes fields-bytes)))))
1532     
1533      ;; read a structure
1534      ((MAT5:structure-class? array-flags) 
1535       (begin
1536         (let*-values (((fieldname-len   fieldname-len-bytes)   (read-fieldname-len eport))
1537                       ((field-names     field-names-bytes)     (read-field-names eport fieldname-len))
1538                       ((fields          fields-bytes)          (read-struct-data eport array-flags array-dims
1539                                                                                  field-names)))
1540                     
1541                      (values  (MAT5:structure array-name array-dims field-names fields)
1542                               (+ flags-bytes dims-bytes name-bytes fieldname-len-bytes
1543                                  field-names-bytes fields-bytes)))))
1544     
1545      ;; read a cell array
1546      ((MAT5:cell-class? array-flags) 
1547       (begin
1548         (let*-values (((cell   cell-bytes)   
1549                        (read-cell-data eport array-flags array-dims)))
1550                      (values  (MAT5:cell-array array-name array-dims cell)
1551                               (+ flags-bytes dims-bytes name-bytes cell-bytes)))))
1552     
1553      ;; read a sparse array
1554      ((MAT5:sparse-class? array-flags) 
1555       (begin
1556         (if (not (= (length array-dims) 2)) 
1557             (MAT5:error eport "read-array supports only two-dimensional sparse arrays"))
1558         (let*-values (((row-index   row-bytes)    (read-row-index eport array-flags))
1559                       ((col-index   col-bytes)    (read-col-index eport array-dims))
1560                       ((data-type real-part   real-bytes)   
1561                        (read-sparse-array-data eport array-flags row-index col-index))
1562                       ((dummy imag-part   imag-bytes)   
1563                        (if (MAT5:complex-array? array-flags)
1564                            (read-sparse-array-data eport array-flags row-index col-index)
1565                            (values #f #f 0))))
1566                      (values  (MAT5:sparse-array
1567                                array-name data-type array-dims row-index col-index real-part imag-part)
1568                               (+ flags-bytes dims-bytes name-bytes row-bytes col-bytes real-bytes imag-bytes)))))
1569     
1570      ;; read a homogeneous numeric array
1571      (else (let*-values (((data-type real-part real-bytes)   
1572                           (read-num-array-data eport array-flags array-dims))
1573                          ((dummy imag-part   imag-bytes)   
1574                           (if (MAT5:complex-array? array-flags)
1575                               (read-num-array-data eport array-flags array-dims)  (values #f #f 0))))
1576                         (values  (MAT5:num-array array-name data-type array-dims real-part imag-part)
1577                                  (+ flags-bytes dims-bytes name-bytes real-bytes imag-bytes)))))))))
1578
1579(define strict-read-data-element
1580  (read-data-element  read-sparse-array-data
1581                      read-num-array-data))
1582
1583
1584; Procedure:
1585; MAT5:read-data-element:: ENDIAN-PORT -> MAT5:DATA-ELEMENT
1586;
1587; Reads a MAT5 data element from the given endian port.
1588;
1589(define (MAT5:read-data-element eport)
1590  (let-values  (((data-type data-size header-bytes)  (read-data-element-header eport)))
1591               (if (not data-type) (MAT5:error eport "unrecognized type word " data-type))
1592               (let ((data
1593                      (let  loop  ((i 0) (lst '()))
1594                        (if (< i data-size)
1595                            (let ((word  (strict-read-data-element eport data-type data-size)))
1596                              (loop (+ i (MAT5:data-element-bytes word)) (cons word lst)))
1597                            (if (= i data-size)
1598                                (reverse lst)
1599                                (MAT5:error eport "data element size mismatch: i = " i " data-size = " data-size))))))
1600                 (align-eport-pos eport)
1601                 data)))
1602
1603 
1604;--------------------------
1605;
1606; MAT-file Writing Routines
1607;
1608
1609; Procedure:
1610; write-pad-bytes:: ENDIAN-PORT ->  BYTES
1611;
1612; Writes pad bytes to a file, so that it is aligned on a 64-bit
1613; boundary. Returns the number of bytes written.
1614;
1615(define (write-pad-bytes  eport)
1616  (let* ((pos  (endian-port-pos eport))
1617         (pos1  (let loop ((pos1 pos))
1618                   (if (not (zero? (modulo pos1 8))) 
1619                       (loop (+ pos1 1)) pos1))))
1620    (if (> pos1 pos)
1621        (let* ((dx   (- pos1 pos))
1622               (bv   (make-byte-vector dx 0)))
1623          (endian-port-write-byte-vector eport bv)
1624          dx)  0)))
1625
1626; Procedure:
1627; MAT5:write-header:: ENDIAN-PORT * STRING * STRING -> UNDEFINED
1628;
1629; Writes a MAT5 header to the given endian port. Arguments TEXT and
1630; SUBSYS are strings. If they are longer than their maximum permitted
1631; lengths (116 and 8, respectively), they will be truncated.
1632;
1633(define (MAT5:write-header eport text subsys)
1634  (let* ((text     (let ((src   (cond ((string? text)  (string-pad-right text MAT5:header-text-size  #\nul))
1635                                      (else     (MAT5:error eport "text argument is of invalid type: " text)))))
1636                         (string->byte-vector src)))
1637         (subsys   (let ((src   (cond ((string? subsys)  (string-pad-right subsys MAT5:header-subsys-size  #\nul))
1638                                      (else     (MAT5:error eport "subsys argument is of invalid type: " subsys))))) 
1639                         (string->byte-vector src)))
1640         (magic    MAT5:msb-magic))
1641    (endian-port-setpos eport 0)
1642    (endian-port-write-byte-vector eport text (MSB))
1643    (endian-port-write-byte-vector eport subsys (MSB))
1644    (endian-port-write-int2 eport MAT5:msb-version (MSB))
1645    (endian-port-write-int2 eport magic)
1646    (write-pad-bytes eport)))
1647
1648
1649; Procedure:
1650; write-num-data-element:: ENDIAN-PORT * MAT5:DATA-TYPE * VALUE -> BYTES
1651;
1652; Writes a data atom (a number or a char) to the given endian
1653; port. Returns the number of bytes written.
1654;
1655(define (write-num-data-element eport type data)
1656  (cases MAT5:data-type type
1657         (miINT8   ()     (endian-port-write-int1 eport  data))
1658         (miUINT8  ()     (endian-port-write-int1 eport  data))
1659         (miINT16  ()     (endian-port-write-int2 eport  data))
1660         (miUINT16 ()     (endian-port-write-int2 eport  data))
1661         (miINT32  ()     (endian-port-write-int4 eport  data))
1662         (miUINT32 ()     (endian-port-write-int4 eport  data))
1663         (miSINGLE ()     (endian-port-write-ieee-float32 eport  data))
1664         (miDOUBLE ()     (endian-port-write-ieee-float64 eport  data))
1665         ;;      (miINT64  ()     (endian-port-write-int8 eport  data) 8)
1666         ;;      (miUINT64 ()     (endian-port-write-int8 eport  data) 8)
1667         (miUTF8   ()     (endian-port-write-int1 eport  data))
1668         (miUTF16  ()     (endian-port-write-int2 eport  data))
1669         (miUTF32  ()     (endian-port-write-int4 eport  data))
1670         (else            (MAT5:error eport "unrecognized type " type))))
1671
1672; Procedure:
1673; write-data-element:: ENDIAN-PORT * MAT5:DATA-TYPE * VALUE -> BYTES
1674;
1675; Writes a data atom or an array/cell/structure to the given endian
1676; port.
1677(define (write-data-element write-sparse-array-data
1678                            write-num-array-data)
1679  (let ((write-array (write-array write-sparse-array-data 
1680                                  write-num-array-data)))
1681    (lambda (eport type data . rest)
1682      (let-optionals  rest ((include-header-bytes #f))
1683        (if (MAT5:numeric-type? type)
1684            (write-num-data-element eport type data)
1685            (cases MAT5:data-type type
1686                   (miMATRIX ()     (write-array eport  data include-header-bytes))
1687                   (else            (MAT5:error eport "unrecognized type " type))))))))
1688 
1689
1690; Procedure:
1691; write-data-element-header
1692;
1693;
1694(define (write-data-element-header eport data-type data-size . rest)
1695  (let-optionals 
1696   rest ((small #f))
1697   (let ((type-word (MAT5:data-type->word data-type)))
1698     (if (and small (<= data-size 4))
1699         (let* ((bytes  (endian-port-write-int2 eport data-size))
1700                (bytes  (+ bytes (endian-port-write-int2 eport type-word))))
1701           bytes)
1702         (let* ((bytes  (endian-port-write-int4 eport type-word))
1703                (bytes  (+ bytes (endian-port-write-int4 eport data-size))))
1704           bytes)))))
1705
1706
1707; Procedure:
1708; write-array-flags
1709;
1710;
1711(define (write-array-flags eport flags class nzmax . rest)
1712  (let-optionals 
1713   rest  ((small #f))
1714   (let* ((data-type   (miUINT32))
1715          (data-size   (* 2 (MAT5:sizeof data-type)))
1716          (flags-word  (bitwise-ior  (if (bit-vector-ref flags mxLOGICAL_FLAG) #b00000010 0)
1717                                     (if (bit-vector-ref flags mxGLOBAL_FLAG)  #b00000100 0)
1718                                     (if (bit-vector-ref flags mxCOMPLEX_FLAG) #b00001000 0)))
1719          (flags-word  (bitwise-ior (arithmetic-shift flags-word 8)
1720                                    (MAT5:array-class->word class)))
1721          (bytes       (write-data-element-header eport data-type  data-size small))
1722          (bytes       (+ bytes (endian-port-write-int4 eport flags-word)))
1723          (bytes       (+ bytes (endian-port-write-int4 eport nzmax)))
1724          (bytes       (+ bytes (write-pad-bytes eport))))
1725     bytes)))
1726     
1727         
1728; Procedure:
1729; write-array-dimensions
1730;
1731;
1732(define (write-array-dimensions eport dims . rest)
1733  (let-optionals rest  ((small #f))
1734   (if (null? dims) (MAT5:error "empty dimension list"))
1735   (if (not (MAT5:dimensions? dims)) (MAT5:error "invalid dimension list"))
1736   (let* ((data-type   (miINT32))
1737          (data-size   (* (length dims) (MAT5:sizeof data-type)))
1738          (bytes       (write-data-element-header eport data-type  data-size small))
1739          (bytes       (+ bytes 
1740                          (let loop ((dims dims) (bytes 0))
1741                            (cond ((not (null? dims))
1742                                   (let* ((dim    (car dims))
1743                                          (x      (if (eq? dim '??) 0 dim))
1744                                          (bytes  (+ bytes (endian-port-write-int4 eport x))))
1745                                     (loop (cdr dims) bytes)))
1746                                  (else bytes)))))
1747          (bytes       (+ bytes (write-pad-bytes eport))))
1748     bytes)))
1749
1750
1751; Procedure:
1752; write-array-name
1753;
1754;
1755(define (write-array-name eport name . rest)
1756  (let-optionals 
1757   rest  ((small #f))
1758   (let* ((bv          (string->byte-vector name))
1759          (data-type   (miINT8))
1760          (data-size   (* (byte-vector-length bv) (MAT5:sizeof data-type)))
1761          (bytes       (write-data-element-header eport data-type  data-size small))
1762          (bytes       (+ bytes (endian-port-write-byte-vector eport bv (MSB))))
1763          (bytes       (+ bytes (write-pad-bytes eport))))
1764     bytes)))
1765   
1766
1767; Procedure:
1768; write-row-index
1769;
1770;
1771(define (write-row-index eport row-index . rest)
1772  (let-optionals 
1773   rest  ((small #f))
1774   (if (null? row-index) (MAT5:error "empty row-index list"))
1775   (let* ((data-type   (miINT32))
1776          (data-size   (* (length row-index) (MAT5:sizeof data-type)))
1777          (bytes       (write-data-element-header eport data-type  data-size small))
1778          (bytes       (+ bytes (let loop ((row-index row-index) (bytes 0))
1779                                  (cond ((not (null? row-index))
1780                                         (let ((bytes (+ bytes (endian-port-write-int4 eport (car row-index)))))
1781                                           (loop (cdr row-index) bytes)))))))
1782          (bytes       (+ bytes (write-pad-bytes eport))))
1783     bytes)))
1784
1785
1786; Procedure:
1787; write-col-index
1788;
1789;
1790(define (write-col-index eport col-index . rest)
1791  (let-optionals 
1792   rest  ((small #f))
1793   (if (null? col-index) (MAT5:error "empty col-index list"))
1794   (let* ((data-type   (miINT32))
1795          (data-size   (* (length col-index) (MAT5:sizeof data-type)))
1796          (bytes       (write-data-element-header eport data-type  data-size small))
1797          (bytes       (+ bytes (let loop ((col-index col-index) (bytes 0))
1798                                  (cond ((not (null? col-index))
1799                                         (let ((bytes (+ bytes (endian-port-write-int4 eport (car col-index)))))
1800                                           (loop (cdr col-index) bytes)))))))
1801          (bytes       (+ bytes (write-pad-bytes eport))))
1802     bytes)))
1803
1804
1805; Procedure:
1806; write-fieldname-len
1807;
1808;
1809(define (write-fieldname-len eport len . rest)
1810  (let-optionals 
1811   rest  ((small #f))
1812   (let* ((data-type   (miINT32))
1813          (data-size   (MAT5:sizeof data-type))
1814          (bytes       (write-data-element-header eport data-type  data-size small))
1815          (bytes       (+ bytes (endian-port-write-int4 eport len)))
1816          (bytes       (+ bytes (write-pad-bytes eport))))
1817     bytes)))
1818   
1819
1820; Procedure:
1821; write-field-names
1822;
1823;
1824(define (write-field-names eport fieldname-len field-names . rest)
1825  (let-optionals 
1826   rest  ((small #f))
1827   (if (null? field-names) (MAT5:error "empty field-names list"))
1828   (let* ((data-type   (miINT8))
1829          (data-size   (* (length field-names) (* fieldname-len (MAT5:sizeof data-type))))
1830          (bytes       (write-data-element-header eport data-type  data-size small))
1831          (bytes       (+ bytes (let loop ((field-names field-names) (bytes 0))
1832                                  (let* ((field-name  (string-pad-right (car field-names) fieldname-len #\nul))
1833                                         (bv          (string->byte-vector field-name)))
1834                                    (let ((bytes  (+ bytes (endian-port-write-byte-vector eport bv (MSB)))))
1835                                      (if (null? (cdr field-names))
1836                                          bytes
1837                                          (loop (cdr field-names) bytes)))))))
1838          (bytes       (+  bytes (write-pad-bytes eport))))
1839     bytes)))
1840
1841; Procedure:
1842; write-fields
1843
1844(define (write-fields write-sparse-array-data 
1845                      write-num-array-data)
1846  (lambda (eport fields)
1847    (let ((write-data-element  (write-data-element write-sparse-array-data
1848                                                   write-num-array-data)))
1849      (let loop ((fields fields) (bytes 0))
1850        (if (null? fields)  (+ bytes (write-pad-bytes eport))
1851            (loop (cdr fields) (+ bytes (write-data-element eport (miMATRIX) (car fields) #t))))))))
1852
1853
1854; Procedure:
1855; write-struct-data
1856;
1857; the same as write-cell-data, only we expect each element of the cell
1858; to be a list of fields (MAT5 data elements), instead of individual
1859; data elements.
1860;
1861(define (write-struct-data  write-sparse-array-data
1862                            write-num-array-data)
1863  (let ((write-fields  (write-fields write-sparse-array-data
1864                                     write-num-array-data)))
1865    (lambda (eport x)
1866      (let* ((write-el    (lambda (i x ax) (+ ax (write-fields eport x))))
1867             (data-size   ((MAT5:array-foldi write-el 'col-major) x 0))
1868             (bytes       (+ data-size (write-pad-bytes eport))))
1869        bytes))))
1870   
1871
1872
1873; Procedure:
1874; write-cell-data
1875;
1876;
1877(define (write-cell-data  write-sparse-array-data
1878                          write-num-array-data)
1879  (let ((write-data-element  (write-data-element write-sparse-array-data
1880                                                 write-num-array-data)))
1881    (lambda (eport x)
1882      (let* ((write-el    (lambda (i x ax) (+ ax (write-data-element eport (miMATRIX) x #t))))
1883             (data-size   ((MAT5:array-foldi write-el 'col-major) x 0))
1884             (bytes       (+ data-size (write-pad-bytes eport))))
1885        bytes))))
1886   
1887
1888; Procedure:
1889; write-sparse-array-data
1890;
1891;
1892(define (write-sparse-array-data eport data-type x)
1893  (write-num-array-data eport data-type x))
1894   
1895
1896; Procedure:
1897; write-num-array-data
1898;
1899;
1900(define (write-num-array-data eport data-type x)
1901  (let ((begin-pos    (endian-port-pos eport))
1902        (bytes        (write-data-element-header eport data-type 0)))
1903    (let* ((write-el    (lambda (i x ax) 
1904                          (+ ax (write-num-data-element eport data-type x))))
1905           (data-size   ((array-foldi write-el 'col-major) x 0))
1906           (bytes       (+ bytes data-size (write-pad-bytes eport)))
1907           (end-pos     (endian-port-pos eport)))
1908      (endian-port-setpos eport begin-pos)
1909      (write-data-element-header eport data-type data-size)
1910      (endian-port-setpos eport end-pos)
1911      (values bytes (array-dimensions x)))))
1912   
1913
1914; Procedure:
1915; write-array
1916;
1917;
1918(define (write-array write-sparse-array-data 
1919                     write-num-array-data)
1920  (lambda (eport x . rest)
1921    (let-optionals  rest ((include-header-bytes #f))
1922    (cases MAT5:array x
1923           (MAT5:object (name dims class-name field-names fields)
1924                        (let* ((flags        (integer->bit-vector 0))
1925                               (begin-pos    (endian-port-pos eport))
1926                               (header-bytes (write-data-element-header eport (miMATRIX) 0))
1927                               (bytes        (write-array-flags eport flags (mxOBJECT_CLASS) 0 #t))
1928                               (bytes        (+ bytes (write-array-dimensions eport dims #t)))
1929                               (bytes        (+ bytes (write-array-name eport name #t)))
1930                               (bytes        (+ bytes (write-array-name eport class-name #t)))
1931                               (bytes        (+ bytes (write-fieldname-len eport MAT5:field-name-length #t)))
1932                               (bytes        (+ bytes (write-field-names eport MAT5:field-name-length
1933                                                                         field-names #t)))
1934                               (bytes        (+ bytes ((write-struct-data write-sparse-array-data
1935                                                                          write-num-array-data) eport x)))
1936                               (end-pos      (endian-port-pos eport)))
1937                          (endian-port-setpos eport begin-pos)
1938                          (write-data-element-header eport (miMATRIX) bytes)
1939                          (endian-port-setpos eport end-pos)
1940                          (if include-header-bytes (+ bytes header-bytes)
1941                              bytes)))
1942           
1943           
1944           (MAT5:structure (name dims field-names fields)
1945                           (let* ((flags        (integer->bit-vector 0))
1946                                  (begin-pos    (endian-port-pos eport))
1947                                  (header-bytes (write-data-element-header eport (miMATRIX) 0))
1948                                  (bytes        (write-array-flags eport flags (mxSTRUCT_CLASS) 0 #t))
1949                                  (bytes        (+ bytes (write-array-dimensions eport dims #t)))
1950                                  (bytes        (+ bytes (write-array-name eport name #t)))
1951                                  (bytes        (+ bytes (write-fieldname-len eport MAT5:field-name-length #t)))
1952                                  (bytes        (+ bytes (write-field-names eport MAT5:field-name-length
1953                                                                            field-names #t)))
1954                                  (bytes        (+ bytes ((write-struct-data write-sparse-array-data
1955                                                                             write-num-array-data) eport x)))
1956                                  (end-pos      (endian-port-pos eport)))
1957                             (endian-port-setpos eport begin-pos)
1958                             (write-data-element-header eport (miMATRIX) bytes)
1959                             (endian-port-setpos eport end-pos)
1960                             (if include-header-bytes (+ bytes header-bytes)
1961                                 bytes)))
1962           
1963           (MAT5:cell-array (name dims cell)
1964                            (let* ((flags        (integer->bit-vector 0))
1965                                   (begin-pos    (endian-port-pos eport))
1966                                   (header-bytes (write-data-element-header eport (miMATRIX) 0))
1967                                   (bytes        (write-array-flags eport flags (mxCELL_CLASS) 0 #t))
1968                                   (bytes        (+ bytes (write-array-dimensions eport dims #t)))
1969                                   (bytes        (+ bytes (write-array-name eport name #t)))
1970                                   (bytes        (+ bytes ((write-cell-data write-sparse-array-data
1971                                                                            write-num-array-data) 
1972                                                           eport x)))
1973                                   (end-pos      (endian-port-pos eport)))
1974                              (endian-port-setpos eport begin-pos)
1975                              (write-data-element-header eport (miMATRIX) bytes)
1976                              (endian-port-setpos eport end-pos)
1977                              (if include-header-bytes (+ bytes header-bytes)
1978                                  bytes)))
1979
1980
1981           (MAT5:sparse-array (name data-type dims row-index col-index real imag) 
1982                              (let* ((flags        (integer->bit-vector 0))
1983                                     (flags        (if imag (bit-vector-set! flags mxCOMPLEX_FLAG #t) flags))
1984                                     (nzmax        (* (length row-index) (length col-index)))
1985                                     (begin-pos    (endian-port-pos eport))
1986                                     (header-bytes (write-data-element-header eport (miMATRIX) 0))
1987                                     (bytes        (write-array-flags eport flags (mxSPARSE_CLASS) nzmax #t))
1988                                     (bytes        (+ bytes (write-array-dimensions eport dims #t)))   
1989                                     (bytes        (+ bytes (write-array-name eport name #t)))
1990                                     (bytes        (+ bytes (write-row-index eport row-index)))
1991                                     (bytes        (+ bytes (write-col-index eport col-index)))
1992                                     (bytes        (+ bytes (write-sparse-array-data eport data-type real)))
1993                                     (bytes        (if imag 
1994                                                       (+ bytes (write-sparse-array-data eport data-type imag))
1995                                                       bytes))
1996                                     (end-pos      (endian-port-pos eport)))
1997                                (endian-port-setpos eport begin-pos)
1998                                (write-data-element-header eport (miMATRIX) bytes)
1999                                (endian-port-setpos eport end-pos)
2000                                (if include-header-bytes (+ bytes header-bytes)
2001                                    bytes)))
2002           
2003           
2004           (MAT5:num-array (name data-type dims real imag)
2005                           (let* ((flags        (integer->bit-vector 0))
2006                                  (flags        (if imag (bit-vector-set! flags mxCOMPLEX_FLAG #t) flags))
2007                                  (class        (MAT5:data-type->array-class data-type))
2008                                  (begin-pos    (endian-port-pos eport))
2009                                  (header-bytes (write-data-element-header eport (miMATRIX) 0))
2010                                  (bytes        (write-array-flags eport flags class 0 #t))
2011                                  (dims-pos     (endian-port-pos eport))
2012                                  (bytes        (+ bytes (write-array-dimensions eport dims #t)))       
2013                                  (bytes        (+ bytes (write-array-name eport name #t))))
2014                             (let-values (((data-size real-dims) (write-num-array-data eport data-type real)))
2015                                         (let* ((bytes (+ bytes data-size))
2016                                                (bytes (if imag 
2017                                                           (let-values 
2018                                                            (((data-size imag-dims)
2019                                                              (write-num-array-data eport data-type imag)))
2020                                                            (+ bytes data-size))
2021                                                           bytes))
2022                                                (end-pos      (endian-port-pos eport)))
2023                                           (endian-port-setpos eport begin-pos)
2024                                           (write-data-element-header eport (miMATRIX) bytes)
2025                                           (endian-port-setpos eport dims-pos)
2026                                           (write-array-dimensions eport real-dims #t)
2027                                           (endian-port-setpos eport end-pos)
2028                                           (if include-header-bytes (+ bytes header-bytes)
2029                                               bytes)))))
2030           
2031           (else (MAT5:error "invalid array type"))))))
2032
2033
2034; Procedure:
2035; write-vector
2036;
2037;
2038(define (write-vector eport type vops vect)
2039  (let* ((vector-length  (alist-ref 'vector-length vops))
2040         (vector-ref     (alist-ref 'vector-ref vops))
2041         (len            (vector-length vect))
2042         (size           (* (MAT5:sizeof type) len))
2043         (bytes          (write-data-element-header eport type size))
2044         (bytes          (+ bytes (let loop ((i 0)  (bytes 0))
2045                                    (if (< i len) bytes
2046                                        (let ((v (vector-ref vect i)))
2047                                          (loop (+ i 1)  (+ bytes (write-num-data-element eport type v)))))))))
2048    (write-pad-bytes  eport)
2049    bytes))
2050
2051
2052; Procedure:
2053; write-string
2054;
2055;
2056(define (write-string eport  str)   
2057  (let* ((type  (miINT8))
2058         (size  (* (MAT5:sizeof type) (length str)))
2059         (bytes  (write-data-element-header eport type size))
2060         (bytes  (+ bytes (let loop ((chars (string->list str))  (bytes 0))
2061                            (if (null? chars) bytes
2062                                (let ((bytes (+ bytes (write-num-data-element eport type (car chars)))))
2063                                  (loop (cdr chars) bytes)))))))
2064    bytes))
2065
2066
2067(define s8vops   `((vector-ref .    ,s8vector-ref)
2068                   (vector-length . ,s8vector-length)))
2069(define u8vops   `((vector-ref .    ,u8vector-ref)
2070                   (vector-length . ,u8vector-length)))
2071(define s16vops  `((vector-ref .   ,s16vector-ref)
2072                   (vector-length . ,s16vector-length)))
2073(define u16vops  `((vector-ref .   ,u16vector-ref)
2074                   (vector-length . ,u16vector-length)))
2075(define s32vops  `((vector-ref .   ,s32vector-ref)
2076                   (vector-length . ,s32vector-length)))
2077(define u32vops  `((vector-ref .   ,u32vector-ref)
2078                   (vector-length . ,u32vector-length)))
2079;;(define s64vops  `((vector-ref .   ,s64vector-ref)
2080;;                 (vector-length . ,s64vector-length)))
2081;;(define u64vops  `((vector-ref .   ,u64vector-ref)
2082;;                 (vector-length . ,u64vector-length)))
2083(define f32vops  `((vector-ref .   ,f32vector-ref)
2084                   (vector-length . ,f32vector-length)))
2085(define f64vops  `((vector-ref .   ,f64vector-ref)
2086                   (vector-length . ,f64vector-length)))
2087
2088
2089
2090(define strict-write-data-element
2091  (write-data-element  write-sparse-array-data
2092                       write-num-array-data))
2093
2094
2095; Procedure:
2096; MAT5:write-data-element:: ENDIAN-PORT * VALUE -> BYTES
2097;
2098;
2099; Writes a MAT5 data element to the given endian port.
2100;
2101(define (MAT5:write-data-element eport data-element)
2102  (cond ((MAT5:array? data-element)  (strict-write-data-element eport (miMATRIX) data-element))
2103       
2104        ((string? data-element)      (write-string eport data-element))
2105       
2106        ((s8vector? data-element)    (write-vector eport (miINT8)   s8vops data-element))
2107
2108        ((u8vector? data-element)    (write-vector eport (miUINT8)  u8vops data-element))
2109
2110        ((s16vector? data-element)   (write-vector eport (miINT16)  s16vops data-element))
2111
2112        ((u16vector? data-element)   (write-vector eport (miUINT16) u16vops data-element))
2113
2114        ((s32vector? data-element)   (write-vector eport (miINT32)  s32vops data-element))
2115
2116        ((u32vector? data-element)   (write-vector eport (miUINT32) u32vops data-element))
2117
2118;;      ((s64vector? data-element)   (write-vector eport (miINT64)  s64vops data-element))
2119
2120;;      ((u64vector? data-element)   (write-vector eport (miUINT64) u64vops data-element))
2121
2122        ((f32vector? data-element)   (write-vector eport (miSINGLE) f32vops data-element))
2123
2124        ((f64vector? data-element)   (write-vector eport (miDOUBLE) f64vops data-element))
2125
2126        (else        (MAT5:error "element " data-element " is of unknown type"))))
2127
2128
2129;---------------------------
2130;
2131; MAT-file stream interface
2132;
2133;
2134
2135; Procedure:
2136; MAT5:read-stream:: ENDIAN-PORT -> STREAM
2137;
2138; Given an endian port, creates an SRFI-40 stream that has the MAT5
2139; header as its first element, followed by all the MAT-File data
2140; elements that can be read from the endian port. Any numeric arrays
2141; that are in the file will be read lazily; i.e. their DATA field will
2142; be another SRFI-40 stream containing the column vectors of the
2143; numeric array.
2144;
2145(define (MAT5:read-stream eport)
2146  (let ((header (MAT5:read-header eport)))
2147    (let loop ((x header) (pos (endian-port-pos eport)))
2148      (stream-delay
2149       (begin (endian-port-setpos eport pos)
2150              (if (endian-port-eof? eport)
2151                  (stream-cons x stream-null)
2152                  (let* ((y    (MAT5:stream-read-data-element eport))
2153                         (pos  (endian-port-pos eport)))
2154                    (stream-cons x (loop y pos)))))))))
2155
2156
2157(define (stream-read-sparse-array-data eport array-flags row-index col-index)
2158  (let-values
2159   (((data-type data-size header-bytes) (read-data-element-header eport)))
2160   (if (not data-type) (MAT5:error eport "invalid data element type"))
2161   (if (not (MAT5:numeric-type? data-type)) (MAT5:error eport "non-numeric sparse array data"))
2162   (if (not (= (MAT5:array-vector-length data-type data-size) 
2163               (* (length row-index) (length col-index))))
2164       (MAT5:error eport "incompatible dimensions: data length is " 
2165                   (MAT5:array-vector-length data-type data-size)
2166                   " array dimensions are specified as " 
2167                   row-index " " col-index))
2168   (let-values   
2169    (((prototype vector-set! vector-ref vector-length vector?) 
2170      (MAT5:array-vector-ops data-type)))
2171    (let ((beg-pos (endian-port-pos eport))
2172          (vector-data (MAT5:vector-make data-type (length col-index))))
2173      (endian-port-setpos eport (+ beg-pos data-size))
2174      (let ((pad-bytes  (align-eport-pos eport)))
2175        (values data-type
2176                (let loop ((len (MAT5:array-vector-length data-type data-size))
2177                           (pos beg-pos))
2178                  (stream-delay
2179                   (if (zero? len)
2180                       (let ((pad-bytes  (align-eport-pos eport)))
2181                         stream-null)
2182                       (begin (endian-port-setpos eport pos)
2183                              (let vector-loop ((i 0) (n (vector-length vector-data)))
2184                                (if (zero? n) 
2185                                    (stream-cons vector-data (loop (- len i) (endian-port-pos eport)))
2186                                    (let ((data-element (read-num-data-element eport data-type data-size)))
2187                                      (vector-set! vector-data i (MAT5:data-element-data data-element))
2188                                      (vector-loop (+ i 1)  (- n 1)))))))))
2189                (+ pad-bytes header-bytes data-size)))))))
2190
2191(define (stream-read-num-array-data eport array-flags array-dims)
2192  (let-values
2193   (((data-type data-size header-bytes) (read-data-element-header eport)))
2194   (if (not data-type) (MAT5:error eport "invalid data element type"))
2195   (if (not (MAT5:numeric-type? data-type)) (MAT5:error eport "non-numeric array data"))
2196   (if (not (= (MAT5:array-vector-length data-type data-size) 
2197               (apply * array-dims)))
2198       (MAT5:error eport "incompatible dimensions: data length is " 
2199                   (MAT5:array-vector-length data-type data-size)
2200                   " array dimensions are specified as " array-dims))
2201   (let-values   
2202    (((prototype vector-set! vector-ref vector-length vector?) 
2203      (MAT5:array-vector-ops data-type)))
2204    (let ((beg-pos (endian-port-pos eport))
2205          (vector-data (MAT5:vector-make data-type (first array-dims))))
2206      (endian-port-setpos eport (+ beg-pos data-size))
2207      (let ((pad-bytes  (align-eport-pos eport)))
2208        (values data-type
2209                (let loop ((len (MAT5:array-vector-length data-type data-size))
2210                           (pos beg-pos))
2211                  (stream-delay
2212                   (if (zero? len)
2213                       (let ((pad-bytes  (align-eport-pos eport))) stream-null)
2214                       (begin (endian-port-setpos eport pos)
2215                              (let vector-loop ((i 0) (n (vector-length vector-data)))
2216                                (if (zero? n) 
2217                                    (begin
2218                                      (stream-cons vector-data (loop (- len i) (endian-port-pos eport))))
2219                                    (let ((data-element (read-num-data-element eport data-type data-size)))
2220                                      (vector-set! vector-data i (MAT5:data-element-data data-element))
2221                                      (vector-loop (+ i 1)  (- n 1)))))))))
2222                (+ pad-bytes header-bytes data-size)))))))
2223
2224
2225(define stream-read-data-element
2226  (read-data-element stream-read-sparse-array-data
2227                     stream-read-num-array-data))
2228
2229
2230; Procedure:
2231; MAT5:stream-read-data-element:: ENDIAN-PORT -> STREAM
2232;
2233; Given an endian port, reads MAT-File data elements from the
2234; port. Any numeric arrays that are in the file will be read lazily;
2235; i.e. their DATA field will be another SRFI-40 stream containing the
2236; column vectors of the numeric array.
2237;
2238(define (MAT5:stream-read-data-element eport)
2239  (let-values  (((data-type data-size header-bytes)  (read-data-element-header eport)))
2240               (if (not data-type) (MAT5:error eport "unrecognized type word " data-type))
2241               (let ((data
2242                      (let  loop  ((i 0) (lst '()))
2243                        (if (< i data-size)
2244                            (let ((word  (stream-read-data-element eport data-type data-size)))
2245                              (loop (+ i (MAT5:data-element-bytes word)) (cons word lst)))
2246                            (if (= i data-size)
2247                                (reverse lst)
2248                                (MAT5:error eport "data element size mismatch"))))))
2249                 (align-eport-pos eport)
2250                 data)))
2251
2252
2253
2254; Procedure:
2255; MAT5:write-stream:: ENDIAN-PORT * STRING * STRING * STREAM -> UNDEFINED
2256;
2257; Given an endian port, and an SRFI-40 stream of MAT5:data-element
2258; objects, write a MAT-File header and all elements in the stream to
2259; the endian port. If any numeric arrays are contained in the input
2260; stream, it is expected that their data elements are encoded as
2261; SRFI-40 streams containing the column vectors of the array.
2262;
2263(define (MAT5:write-stream eport text subsys strm)
2264  (begin
2265    (MAT5:write-header eport text subsys)
2266    (stream-for-each (lambda (x) (MAT5:stream-write-data-element eport x)) strm)))
2267
2268
2269; Procedure:
2270; stream-write-sparse-array-data
2271;
2272;
2273(define (stream-write-sparse-array-data eport data-type x)
2274  (stream-write-num-array-data eport data-type x))
2275
2276; Procedure:
2277; stream-write-num-array-data
2278;
2279;
2280(define (stream-write-num-array-data eport data-type x)
2281  (let ((dim++    (lambda (dims) (cons (+ 1 (car dims)) (cdr dims))))
2282        (col-len  0))
2283   
2284    (let-values   
2285     (((prototype vector-set! vector-ref vector-length vector?) 
2286       (MAT5:array-vector-ops data-type)))
2287     (letrec ((write-col  (lambda (elm ax) 
2288                            (let-values (((data-size dims top?) (split3 ax)))
2289                             (cond
2290                              ((stream? elm) 
2291                               (if top?
2292                                   (let-values 
2293                                    (((data-size dims) 
2294                                      (split2 (stream-fold write-col 
2295                                                           (list data-size (list 0 (+ 1 (last dims))) #f) elm))))
2296                                               (list data-size dims #t))
2297                                   (stream-fold write-col (list data-size (cons 0 (dim++ dims)) #f) elm)))
2298                             
2299                              ((vector? elm)
2300                               (let ((len (vector-length elm)))
2301                                 (set! col-len len)
2302                                 (let loop ((i 0)  (bytes 0))
2303                                   (if (>= i len)  (list (+ data-size bytes) (dim++ dims) #f)
2304                                       (let ((v (vector-ref elm i)))
2305                                         (loop (+ i 1)  (+ bytes (write-num-data-element eport data-type v))))))))
2306                              (else (MAT5:error "invalid array data type")))))))
2307     (let* ((begin-pos    (endian-port-pos eport))
2308            (bytes        (write-data-element-header eport data-type 0)))
2309       (let-values (((data-size dims)  (split2 (stream-fold write-col (list 0 '(0) #t) x))))
2310                   (let ((bytes     (+ bytes data-size (write-pad-bytes eport)))
2311                         (end-pos   (endian-port-pos eport)))
2312                     (endian-port-setpos eport begin-pos)
2313                     (write-data-element-header eport data-type data-size)
2314                     (endian-port-setpos eport end-pos)
2315                     (values bytes (cons col-len dims)))))))))
2316   
2317
2318(define stream-write-data-element
2319  (write-data-element stream-write-sparse-array-data
2320                      stream-write-num-array-data))
2321
2322
2323; Procedure:
2324; MAT5:stream-write-data-element:: ENDIAN-PORT * VALUE -> BYTES
2325;
2326; Writes a MAT5 data element to the given endian port. If the input
2327; element is a MAT5 array, any numeric array data must be encoded as
2328; an SRFI-40 stream, where each element of the stream is a column
2329; vector of the array.
2330(define (MAT5:stream-write-data-element eport data-element)
2331  (cond ((MAT5:array? data-element)  (stream-write-data-element eport (miMATRIX) data-element))
2332       
2333        ((string? data-element)      (write-string eport data-element))
2334       
2335        ((s8vector? data-element)    (write-vector eport (miINT8)   s8vops data-element))
2336
2337        ((u8vector? data-element)    (write-vector eport (miUINT8)  u8vops data-element))
2338
2339        ((s16vector? data-element)   (write-vector eport (miINT16)  s16vops data-element))
2340
2341        ((u16vector? data-element)   (write-vector eport (miUINT16) u16vops data-element))
2342
2343        ((s32vector? data-element)   (write-vector eport (miINT32)  s32vops data-element))
2344
2345        ((u32vector? data-element)   (write-vector eport (miUINT32) u32vops data-element))
2346
2347;;      ((s64vector? data-element)   (write-vector eport (miINT64)  s64vops data-element))
2348
2349;;      ((u64vector? data-element)   (write-vector eport (miUINT64) u64vops data-element))
2350
2351        ((f32vector? data-element)   (write-vector eport (miSINGLE) f32vops data-element))
2352
2353        ((f64vector? data-element)   (write-vector eport (miDOUBLE) f64vops data-element))
2354
2355        (else        (MAT5:error "element " data-element " is of unknown type"))))
Note: See TracBrowser for help on using the repository browser.