source: project/mat5-lib/trunk/mat5-lib.scm @ 4857

Last change on this file since 4857 was 4857, checked in by Ivan Raikov, 13 years ago

License upgrade to GPL v3.

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