Fix test in statistics module.
[scilab.git] / scilab / modules / matio / macros / ReadmiMatrix.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2002-2004 - INRIA - Vincent COUVERT
3 // Copyright (C) ???? - INRIA - Serge STEER
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING, which
7 // you should have received as part of this distribution.  The terms
8 // are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
10
11 function [value,ArrayName]=ReadmiMatrix(fd)
12     // Read a variable in a Matlab binary file
13     // This function has been developed following the 'MAT-File Format' description:
14     // www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/matfile_format.pdf
15     // Copyright INRIA
16     // Authors: SS, VC
17
18     [DataType,NumberOfBytes,Compressed]=ReadTag(fd);
19     if meof(fd) then value=[],ArrayName="",return,end
20     if DataType<>miMatrix then
21         error(msprintf(gettext("Found Datatype=%d, expecting %d."),DataType,miMatrix));
22     end
23     if NumberOfBytes==0 then value=[],return,end
24     [Flags,Class,NnzMax]=ReadArrayFlags(fd);
25     DimensionArray=ReadDimensionArray(fd);
26     ArrayName=ReadArrayName(fd)
27     select Class
28     case DoubleClass
29         value=double(ReadSimpleElement(fd,prod(DimensionArray),Class))
30         if Flags(1) then
31             value=double(value)+%i*double(ReadSimpleElement(fd,prod(DimensionArray)))
32         end
33         value=matrix(value,DimensionArray)
34     case SingleClass
35         value=ReadSimpleElement(fd,prod(DimensionArray),Class)
36         if Flags(1) then
37             value=double(value)+%i*double(ReadSimpleElement(fd,prod(DimensionArray)))
38         end
39         value=matrix(value,DimensionArray)
40     case Int8Class
41         value=int8(ReadSimpleElement(fd,prod(DimensionArray),Class))
42         if Flags(1) then
43             value=double(value)+%i*double(ReadSimpleElement(fd,prod(DimensionArray)))
44         end
45         value=matrix(value,DimensionArray)
46     case Uint8Class
47         value=uint8(ReadSimpleElement(fd,prod(DimensionArray),Class))
48         if Flags(1) then
49             value=double(value)+%i*double(ReadSimpleElement(fd,prod(DimensionArray)))
50         end
51         value=matrix(value,DimensionArray)
52     case Int16Class
53         value=int16(ReadSimpleElement(fd,prod(DimensionArray),Class))
54         if Flags(1) then
55             value=double(value)+%i*double(ReadSimpleElement(fd,prod(DimensionArray)))
56         end
57         value=matrix(value,DimensionArray)
58     case Uint16Class
59         value=uint16(ReadSimpleElement(fd,prod(DimensionArray),Class))
60         if Flags(1) then
61             value=double(value)+%i*double(ReadSimpleElement(fd,prod(DimensionArray)))
62         end
63         value=matrix(value,DimensionArray)
64     case Int32Class
65         value=int32(ReadSimpleElement(fd,prod(DimensionArray),Class))
66         if Flags(1) then
67             value=double(value)+%i*double(ReadSimpleElement(fd,prod(DimensionArray)))
68         end
69         value=matrix(value,DimensionArray)
70     case Uint32Class
71         value=uint32(ReadSimpleElement(fd,prod(DimensionArray),Class))
72         if Flags(1) then
73             value=double(value)+%i*double(ReadSimpleElement(fd,prod(DimensionArray)))
74         end
75         value=matrix(value,DimensionArray)
76     case CellClass
77
78         entries=list()
79         for k=1:prod(DimensionArray)
80             entries(k)=ReadmiMatrix(fd)
81         end
82         value=mlist(["ce","dims","entries"],int32(DimensionArray),entries)
83     case CharClass
84         value=matrix(ReadSimpleElement(fd,prod(DimensionArray)),DimensionArray(1),-1)
85         t=[];for v=value',t=[t;stripblanks(ascii(double(v)))];end
86         value=t
87     case StructClass
88         FieldNameLength=double(ReadSimpleElement(fd,1))
89         FieldNames=matrix(ReadSimpleElement(fd),FieldNameLength,-1)
90         NumberOfFields=size(FieldNames,2)
91         Fnams=[];Fields=list();
92         for k=1:NumberOfFields
93             l=find(FieldNames(:,k)==0,1)-1;
94             Fnams=[Fnams,stripblanks(ascii(double(FieldNames(1:l,k))))];
95             Fields(k)=list();
96         end
97
98         if prod(DimensionArray)==1 then
99             for k=1:NumberOfFields
100                 Fields(k)=ReadmiMatrix(fd);
101             end
102         else
103             Fk=list();for i=1:size(DimensionArray,"*"),Fk(i)=[];end
104             for k=1:NumberOfFields,Fields(k)=Fk,end
105             for i=1:prod(DimensionArray)
106                 for k=1:NumberOfFields
107                     Fields(k)(i)=ReadmiMatrix(fd);
108                 end
109             end
110         end
111         //Form Scilab representation
112         value=mlist(["st" "dims" Fnams],int32(DimensionArray),Fields(:))
113     case ObjectClass
114         ClassName=stripblanks(ascii(double(ReadSimpleElement(fd))))
115         FieldNameLength=double(ReadSimpleElement(fd,1))
116         FieldNames=matrix(ReadSimpleElement(fd),FieldNameLength,-1)
117         NumberOfFields=size(FieldNames,2)
118         Fields=list();Fnams=[]
119         for k=1:NumberOfFields
120             l=find(FieldNames(:,k)==0,1)-1
121             Fnams=[Fnams,stripblanks(ascii(double(FieldNames(1:l,k))))]
122             Fields(k)=ReadmiMatrix(fd)
123         end
124         //Form Scilab representation
125         value=tlist([ClassName, Fnams],Fields(:))
126         select ClassName
127         case "inline" then
128             value=Object2Inline(value)
129         case "ss" then
130             value=Object2SS(value)
131         case "tf" then
132             value=Object2tf(value)
133         end
134     case SparseClass then
135         RowIndex=double(ReadSimpleElement(fd,NnzMax))
136         ColumnIndex=double(ReadSimpleElement(fd,DimensionArray(2)+1))
137         value=double(ReadSimpleElement(fd))
138         if Flags(1) then
139             value=value+%i*double(ReadSimpleElement(fd))
140         end
141
142         //Form Scilab representation
143         ptr=ColumnIndex(2:$)-ColumnIndex(1:$-1);
144         col=[];cc=1;
145         for ic=1:size(ptr,"*")
146             col=[col;cc(ones(ptr(ic),1))];cc=cc+1;
147         end
148         //in some cases the initial value of ne is  bigger than necessary
149         ne=min(size(RowIndex,"*"),size(col,"*"));
150         RowIndex=RowIndex(1:ne);col=col(1:ne);
151         if RowIndex<>[] then RowIndex=RowIndex(:)+1,end
152         value=sparse([col(:),RowIndex],value(:),DimensionArray([2 1])).'
153     else
154         error(gettext("Unknown Class."));
155     end
156 endfunction
157
158 function [DataType,NumberOfBytes,Compressed]=ReadTag(fd)
159     //--TAG
160     //Copyright INRIA
161     //Author Serge Steer
162     p1=mtell(fd)
163
164     t=mget(2,md_s,fd);
165     if t==[] then //EOF
166         DataType=0;NumberOfBytes=0,Compressed=%f
167     else
168         if endian=="l" then t=t([2 1]),end
169         Compressed=t(1)<>0;
170         if Compressed then // compressed data element format
171             NumberOfBytes=t(1)
172             DataType=t(2)
173         else
174             mseek(p1,fd)
175             DataType=mget(1,md_i,fd);
176             NumberOfBytes=mget(1,md_i,fd);
177         end
178     end
179 endfunction
180
181
182
183 function [Flags,Class,NnzMax]=ReadArrayFlags(fd)
184     //Copyright INRIA
185     //Author Serge Steer
186     [DataType,NumberOfBytes,Compressed]=ReadTag(fd)
187     B=mget(4,"uc",fd);
188     if endian=="l" then B=B([4 3 2 1]),end
189     Class=B(4)
190     Flags=byte2bits(B(3));Flags=Flags(4:-1:2)
191     NnzMax=mget(1,md_i,fd)
192 endfunction
193
194 function dims=ReadDimensionArray(fd)
195     //Copyright INRIA
196     //Author Serge Steer
197     dims=double(ReadSimpleElement(fd))
198 endfunction
199
200 function ArrayName=ReadArrayName(fd)
201     //Copyright INRIA
202     //Author Serge Steer
203     ArrayName=ascii(double(ReadSimpleElement(fd)))
204 endfunction
205
206 function value=ReadSimpleElement(fd,NumberOfValues,Class)
207     //Copyright INRIA
208     //Author Serge Steer
209     pse=mtell(fd)
210     [DataType,NumberOfBytes,Compressed]=ReadTag(fd)
211     select DataType
212     case miDOUBLE
213         if argn(2)==1 then NumberOfValues=NumberOfBytes/8,end
214         value=mget(NumberOfValues,md_d,fd)
215     case miSINGLE
216         if argn(2)==1 then NumberOfValues=NumberOfBytes/4,end
217         value=mget(NumberOfValues,md_f,fd)
218     case miINT8
219         if argn(2)==1 then NumberOfValues=NumberOfBytes,end
220         value=mgeti(NumberOfValues,"c",fd)
221     case miUINT8
222         if argn(2)==1 then NumberOfValues=NumberOfBytes,end
223         value=mgeti(NumberOfValues,"uc",fd)
224     case miINT16
225         if argn(2)==1 then NumberOfValues=NumberOfBytes/2,end
226         value=mgeti(NumberOfValues,md_s,fd)
227     case miUINT16
228         if argn(2)==1 then NumberOfValues=NumberOfBytes/2,end
229         value=mget(NumberOfValues,"u"+md_s,fd)
230     case miUINT32
231         if argn(2)==1 then NumberOfValues=NumberOfBytes/4,end
232         value=mgeti(NumberOfValues,"u"+md_i,fd)
233     case miINT32
234         if argn(2)==1 then NumberOfValues=NumberOfBytes/4,end
235         value=mgeti(NumberOfValues,md_i,fd)
236     case miUINT64
237         if argn(2)==1 then NumberOfValues=NumberOfBytes/8,end
238         value=mget(NumberOfValues,"u"+md_l,fd)
239     case miINT64
240         if argn(2)==1 then NumberOfValues=NumberOfBytes/8,end
241         value=mget(NumberOfValues,md_l,fd)
242     case miMatrix
243         mseek(pse,fd)
244         [value,ArrayName]=ReadmiMatrix(fd)
245     else
246         error(msprintf(gettext("Not implemented DataType: %d."),DataType));
247     end
248     padding()
249
250 endfunction
251
252
253 function padding()
254     // skip padding data
255     //----------------------------------------------
256     //Copyright INRIA
257     //Author Serge Steer
258
259     //data fields are aligned on double words
260     np=modulo(8-modulo(mtell(fd),8),8)
261     if np>0 then mget(np,"uc",fd),end
262 endfunction
263
264 function showbin(n,pi)
265     //for debugging purpose
266     //----------------------------------------------
267     //Copyright INRIA
268     //Author Serge Steer
269
270     p=mtell(fd)
271     if argn(2)==2 then mseek(pi,fd),end
272     x=string(matrix(mgeti(8*n,"uc",fd),8,-1)')
273     t=emptystr(n,1)+"|"
274     for k=1:4
275         t=t+part(x(:,k),1:max(length(x(:,k)))+1)
276     end
277     t=t+"|"
278     for k=5:8
279         t=t+part(x(:,k),1:max(length(x(:,k)))+1)
280     end
281     t=t+"|"
282     write(%io(2),t,"(a)")
283     mseek(p,fd)
284 endfunction
285
286
287 function [head,version,swap]=matfile_header(fd)
288     //get the mat file header information
289     //Copyright INRIA
290     //Author Serge Steer
291
292     head=ascii(mget(124,"uc",fd))
293     version=mget(2,"uc",fd)
294     //Magic number endian coding
295     IM_MI=mget(2,"uc",fd);
296     if and(IM_MI==[73,77]) then // little endian file
297         swap="l"
298     elseif and(IM_MI==[77,73]) then // big endian file
299         swap="b"
300     else
301         mclose(fd);
302         // This line has to be mofified according to message in 'loadmatfile' function
303         error(gettext("Invalid level 5 binary MAT-file!."));
304     end
305 endfunction
306
307 function LoadMatConstants()
308     //set constants. This function should be exec'ed
309     //Copyright INRIA
310     //Author Serge Steer
311
312     miINT8=1
313     miUINT8=2
314     miINT16=3
315     miUINT16=4
316     miINT32=5
317     miUINT32=6
318     miSINGLE=7
319     //
320     miDOUBLE=9
321     //
322     //
323     miINT64=12
324     miUINT64=13
325     miMatrix=14
326
327     CellClass=1
328     StructClass=2
329     ObjectClass=3
330     CharClass=4
331     SparseClass=5
332     DoubleClass=6
333     SingleClass=7
334     Int8Class=8
335     Uint8Class=9
336     Int16Class=10
337     Uint16Class=11
338     Int32Class=12
339     Uint32Class=13
340
341     //--set various reading format
342     md_i="i"+endian;md_d="d"+endian;md_s="s"+endian;md_l="l"+endian;md_f="f"+endian;
343
344 endfunction
345
346 function value=Object2Inline(value)
347     //convert inline object to scilab function
348     //Copyright INRIA
349     //Author Serge Steer
350
351     deff("ans=value("+strcat(stripblanks(value.args),",")+")",value.expr,"n")
352     comp(value,1);code=macr2lst(value)
353     load SCI/modules/m2sci/macros/lib
354     killed=[];quote="''";dquote="""";batch=%f
355     [value,trad]=m2sci(code,"value",%f,%f)
356     value($)="endfunction"
357     //define the final version
358     execstr(value)
359 endfunction
360
361 function res=Object2SS(res)
362     //convert ss object to scilab 'lss'
363     //Copyright INRIA
364     //Author Serge Steer
365     A=res.a;if type(A)==17 then A=A.entries(1),end
366     B=res.b;if type(B)==17 then B=B.entries(1),end
367     C=res.c;if type(C)==17 then C=C.entries(1),end
368     D=res.d;if type(D)==17 then D=D.entries(1),end
369     E=res.e;if type(E)==17 then E=E.entries(1),end
370     st_nam=res.StateName
371     props=res.lti
372     dt=props.Ts;if dt==0 then dt="c",end
373     res=syslin(dt,A,B,C,D)
374     res($+1)=props
375     res(1)($+1)="Properties"
376 endfunction
377
378 function res=Object2tf(res)
379     //convert tf object to scilab 'r'
380     //Copyright INRIA
381     //Author Serge Steer
382     v=res.Variable
383     dims=double(res.num.dims) //res.num.dims may be an integer array
384     props=res.lti
385     num=[];den=[];
386     for k=1:prod(dims)
387         num=[num;poly(res.num.entries(k)($:-1:1),v,"c")];
388         den=[den;poly(res.den.entries(k)($:-1:1),v,"c")];
389     end
390     num=matrix(num,dims)
391     den=matrix(den,dims)
392     dt=props.Ts;if dt==0 then dt="c",end
393     res=syslin(dt,num,den)
394     res(1)($+1)="Properties"
395     res($+1)=props
396 endfunction
397
398 function fd=open_matfile(fil)
399     //Copyright INRIA
400     //Author Serge Steer
401     fil=stripblanks(fil)
402     fd=mopen(fil,"rb",0)
403 endfunction
404
405 function b=int2bytes(i)
406     //Copyright INRIA
407     //Author Serge Steer
408     it=inttype(i);it1=modulo(it,10)
409     if it1==1 then
410         b=i(:)
411     else
412         s=iconvert(2^(4*(0:it1-1)),it)
413         d=i;b=s(:);
414         for k=1:it1
415             x=s(it1-k+1);b(k) = d/x; d = d-b(k)*x;
416         end
417     end
418 endfunction
419
420 function b=byte2bits(i)
421     //Copyright INRIA
422     //Author Serge Steer
423     b=(iconvert(i,11)&iconvert(2^(0:3),11))<>uint8(0)
424 endfunction
425
426 function I=columnfirstorder(d)
427     nd=size(d,"*")
428     if nd==2 then
429         I=matrix(matrix(1:prod(d),d)',1,-1)
430     else
431         dd=prod(d(3:$))
432         I1=matrix(1:prod(d),d(1),d(2),dd)
433         I=[]
434         for k=1:dd
435             I=[I matrix(I1(:,:,k)',1,-1)]
436         end
437     end
438 endfunction