8a507bc72f9f2acac72768c7ba7c7c2f719eefba
[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-en.txt
10
11 function [value,ArrayName]=ReadmiMatrix(fd)
12 // Read a variable in a Matlab binary file
13 // This function has been developped 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 informations
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