Edinburgh Speech Tools  2.4-release
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends Pages
EST_DMatrix.cc
1 /*************************************************************************/
2 /* */
3 /* Centre for Speech Technology Research */
4 /* University of Edinburgh, UK */
5 /* Copyright (c) 1995,1996 */
6 /* All Rights Reserved. */
7 /* */
8 /* Permission is hereby granted, free of charge, to use and distribute */
9 /* this software and its documentation without restriction, including */
10 /* without limitation the rights to use, copy, modify, merge, publish, */
11 /* distribute, sublicense, and/or sell copies of this work, and to */
12 /* permit persons to whom this work is furnished to do so, subject to */
13 /* the following conditions: */
14 /* 1. The code must retain the above copyright notice, this list of */
15 /* conditions and the following disclaimer. */
16 /* 2. Any modifications must be clearly marked as such. */
17 /* 3. Original authors' names are not deleted. */
18 /* 4. The authors' names are not used to endorse or promote products */
19 /* derived from this software without specific prior written */
20 /* permission. */
21 /* */
22 /* THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK */
23 /* DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING */
24 /* ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT */
25 /* SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE */
26 /* FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES */
27 /* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN */
28 /* AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, */
29 /* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF */
30 /* THIS SOFTWARE. */
31 /* */
32  /*************************************************************************/
33  /* */
34  /* Author : Simon King */
35  /* Date : February 1999 */
36  /* --------------------------------------------------------------------- */
37  /* Double matrix class - copied from DMatrix ! */
38  /* */
39  /*************************************************************************/
40 
41 #include <cstdlib>
42 #include <cstdio>
43 #include <fstream>
44 #include <cmath>
45 #include <climits>
46 #include "EST_String.h"
47 #include "EST_types.h"
48 #include "EST_FileType.h"
49 #include "EST_Option.h"
50 #include "EST_DMatrix.h"
51 #include "EST_cutils.h" // for swap functions
52 #include "EST_Token.h"
53 #include "rateconv.h"
54 
55 EST_String EST_DMatrix::default_file_type = "est_ascii";
56 
58 :EST_TSimpleMatrix<double>(a.num_rows(), a.num_columns())
59 {
60  double vv = 0.0;
61  if (b < 0)
62  return;
63  if (b == 0)
64  fill(vv);
65 }
66 
68 {
69  int i, j;
70  if (a.num_columns() != num_columns())
71  {
72  cerr <<"Matrix addition error: bad number of columns\n";
73  return *this;
74  }
75  if (a.num_rows() != num_rows())
76  {
77  cerr <<"Matrix addition error: bad number of rows\n";
78  return *this;
79  }
80  for (i = 0; i < num_rows(); ++i)
81  for (j = 0; j < num_columns(); ++j)
82  a_no_check(i, j) += a.a_no_check(i,j);
83 
84  return *this;
85 }
86 
88 {
89  int i, j;
90  if (a.num_columns() != num_columns())
91  {
92  cerr <<"Matrix subtraction error: bad number of columns\n";
93  return *this;
94  }
95  if (a.num_rows() != num_rows())
96  {
97  cerr <<"Matrix subtraction error: bad number of rows\n";
98  return *this;
99  }
100  for (i = 0; i < num_rows(); ++i)
101  for (j = 0; j < num_columns(); ++j)
102  a_no_check(i, j) -= a.a_no_check(i,j);
103 
104  return *this;
105 }
106 
108 {
109 
110  int i,j;
111  for (i = 0; i < num_rows(); ++i)
112  for (j = 0; j < num_columns(); ++j)
113  a_no_check(i, j) *= f;
114 
115  return *this;
116 }
117 
119 {
120 
121  int i,j;
122  for (i = 0; i < num_rows(); ++i)
123  for (j = 0; j < num_columns(); ++j)
124  a_no_check(i, j) /= f;
125 
126  return *this;
127 }
128 
129 EST_DMatrix operator+(const EST_DMatrix &a, const EST_DMatrix &b)
130 {
131  EST_DMatrix ab;
132  int i, j;
133  if (a.num_columns() != b.num_columns())
134  {
135  cerr <<"Matrix addition error: bad number of columns\n";
136  return ab;
137  }
138  if (a.num_rows() != b.num_rows())
139  {
140  cerr <<"Matrix addition error: bad number of rows\n";
141  return ab;
142  }
143  ab.resize(a.num_rows(), a.num_columns());
144  for (i = 0; i < a.num_rows(); ++i)
145  for (j = 0; j < a.num_columns(); ++j)
146  ab.a_no_check(i, j) = a.a_no_check(i, j) + b.a_no_check(i, j);
147 
148  return ab;
149 }
150 
151 EST_DMatrix operator-(const EST_DMatrix &a,const EST_DMatrix &b)
152 {
153  EST_DMatrix ab;
154  int i, j;
155 
156  if (a.num_columns() != b.num_columns())
157  {
158  cerr <<"Matrix subtraction error: bad number of columns:" <<
159  a.num_columns() << " and " << b.num_columns() << endl;
160  return ab;
161  }
162  if (a.num_rows() != b.num_rows())
163  {
164  cerr <<"Matrix subtraction error: bad number of rows\n";
165  return ab;
166  }
167  ab.resize(a.num_rows(), a.num_columns());
168  for (i = 0; i < a.num_rows(); ++i)
169  for (j = 0; j < a.num_columns(); ++j)
170  ab.a_no_check(i, j) = a.a_no_check(i, j) - b.a_no_check(i, j);
171 
172  return ab;
173 }
174 
175 EST_DMatrix operator*(const EST_DMatrix &a, const double x)
176 {
177  EST_DMatrix b(a, 0);
178  int i, j;
179 
180  for (i = 0; i < a.num_rows(); ++i)
181  for (j = 0; j < a.num_columns(); ++j)
182  b.a_no_check(i,j) = a.a_no_check(i,j) * x;
183 
184  return b;
185 }
186 
187 EST_DVector operator*(const EST_DMatrix &a, const EST_DVector &v)
188 {
189  // treat the vector as a column vector
190  // multiply each row of the matrix in turn by the vector
191 
192  EST_DVector b;
193  b.resize(a.num_rows());
194 
195  if(a.num_columns() != v.n())
196  {
197  cerr <<"Matrix-vector multiplication error: matrix rows != vector size"
198  << endl;
199  return b;
200  }
201 
202  int i, j;
203  for (i = 0; i < a.num_rows(); ++i){
204  b[i] = 0.0;
205  for (j = 0; j < a.num_columns(); ++j)
206  b.a_no_check(i) += a.a_no_check(i,j) * v.a_no_check(j);
207  }
208  return b;
209 }
210 
211 EST_DVector operator+(const EST_DVector &a, const EST_DVector &b)
212 {
213  EST_DVector ab;
214  int i;
215  if (a.length() != b.length())
216  {
217  cerr <<"Vector addition error: mismatched lengths\n";
218  return ab;
219  }
220 
221  ab.resize(a.length());
222  for (i = 0; i < a.length(); ++i)
223  ab.a_no_check(i) = a.a_no_check(i) + b.a_no_check(i);
224 
225  return ab;
226 }
227 
228 EST_DVector operator-(const EST_DVector &a, const EST_DVector &b)
229 {
230  EST_DVector ab;
231  int i;
232  if (a.length() != b.length())
233  {
234  cerr <<"Vector subtraction error: mismatched lengths\n";
235  return ab;
236  }
237 
238  ab.resize(a.length());
239  for (i = 0; i < a.length(); ++i)
240  ab.a_no_check(i) = a.a_no_check(i) - b.a_no_check(i);
241 
242  return ab;
243 }
244 
245 
246 EST_DVector operator*(const EST_DVector &v,const EST_DMatrix &a)
247 {
248  // treat the vector as a row vector
249  // multiply the vector by each column of the matrix in turn
250 
251  EST_DVector b;
252  b.resize(a.num_columns());
253 
254  if(a.num_columns() != v.n())
255  {
256  cerr <<"Matrix-vector multiplication error: matrix rows != vector size"
257  << endl;
258  return b;
259  }
260 
261  int i, j;
262  for (j = 0; j < a.num_columns(); ++j){
263  b[j] = 0.0;
264  for (i = 0; i < a.num_rows(); ++i)
265  b.a_no_check(i) += a.a_no_check(i,j) * v.a_no_check(j);
266  }
267  return b;
268 }
269 
270 
271 #if 0
272 EST_DMatrix operator/(const EST_DMatrix &a, double x)
273 {
274  return (a * (1/x));
275 }
276 #endif
277 
278 EST_DMatrix operator*(const EST_DMatrix &a, const EST_DMatrix &b)
279 {
280  EST_DMatrix ab;
281  multiply(a,b,ab);
282  return ab;
283 }
284 
285 void multiply(const EST_DMatrix &a, const EST_DMatrix &b, EST_DMatrix &ab)
286 {
287 
288  if (a.num_columns() != b.num_rows())
289  {
290  cerr <<"Matrix multiply error: a.num_columns() != b.num_rows()\n";
291  return;
292  }
293 
294  ab.resize(a.num_rows(), b.num_columns());
295  int i, j, k, n;
296  n = a.num_columns(); // could also be b.num_rows()
297 
298  for (i = 0; i < a.num_rows(); ++i)
299  for (k = 0; k < b.num_columns(); ++k)
300  {
301  ab.a_no_check(i, k) = 0.0;
302  for (j = 0; j < n; ++j)
303  ab.a_no_check(i, k) +=
304  a.a_no_check(i, j) * b.a_no_check(j, k);
305  }
306 }
307 
308 void EST_DMatrix::copyin(double **inx, int rows, int cols)
309 {
310  int i, j;
311 
312  resize(rows, cols);
313 
314  for (i = 0; i < rows; ++i)
315  for (j = 0; j < cols; ++j)
316  a_no_check(i,j) = inx[i][j];
317 
318 }
319 
320 EST_write_status EST_DMatrix::save(const EST_String &filename,
321  const EST_String &type)
322 {
323 
324  if ((type == "est_ascii") || (type == "est_binary"))
325  return est_save(filename,type);
326  else
327  { // the old stuff raw unheadered
328  int i, j;
329  ostream *outf;
330  if (filename == "-")
331  outf = &cout;
332  else
333  outf = new ofstream(filename);
334 
335  outf->precision(25);
336  if (!(*outf))
337  {
338  cerr << "DMatrix: can't open file \"" << filename
339  <<"\" for writing" << endl;
340  return misc_write_error;
341  }
342 
343  for (i = 0; i < num_rows(); ++i)
344  {
345  for (j = 0; j < num_columns(); ++j)
346  *outf << a_no_check(i,j) << " ";
347  *outf << endl;
348  }
349 
350  if (outf != &cout)
351  delete outf;
352 
353  return write_ok;
354  }
355 }
356 
357 EST_write_status EST_DMatrix::est_save(const EST_String &filename,
358  const EST_String &type)
359 {
360  // Binary save with short header for byte swap and sizes
361  int i,j;
362  FILE *fd;
363  if (filename == "-")
364  fd = stdout;
365  else if ((fd = fopen(filename, "wb")) == NULL)
366  {
367  cerr << "EST_DMatrix: binsave: failed to open \"" << filename <<
368  "\" for writing" << endl;
369  return misc_write_error;
370  }
371 
372  fprintf(fd,"EST_File dmatrix\n");
373  fprintf(fd,"version 1\n");
374  if (type == "est_binary")
375  {
376  fprintf(fd,"DataType binary\n");
377  if (EST_LITTLE_ENDIAN)
378  fprintf(fd,"ByteOrder LittleEndian\n");
379  else
380  fprintf(fd,"ByteOrder BigEndian\n");
381  }
382  else
383  fprintf(fd,"DataType ascii\n");
384 
385  fprintf(fd,"rows %d\n",num_rows());
386  fprintf(fd,"columns %d\n",num_columns());
387 
388  fprintf(fd,"EST_Header_End\n");
389 
390  if (type == "est_binary")
391  {
392  for (i = 0; i < num_rows(); ++i)
393  for (j=0; j < num_columns(); j++)
394  if (fwrite(&a_no_check(i,j),sizeof(double),1,fd) != 1)
395  {
396  cerr << "EST_DMatrix: binsave: failed to write row "
397  << i << " column " << j
398  << " to \"" << filename << "\"" << endl;
399  return misc_write_error;
400  }
401  }
402  else
403  { // est_ascii
404  for (i = 0; i < num_rows(); ++i)
405  {
406  for (j=0; j < num_columns(); j++)
407  fprintf(fd,"%.25f ",a_no_check(i,j));
408  fprintf(fd,"\n");
409  }
410  }
411 
412  if (fd != stdout)
413  fclose(fd);
414 
415  return write_ok;
416 }
417 
418 EST_read_status EST_DMatrix::est_load(const EST_String &filename)
419 {
420 
421  // ascii/binary load with short header for byte swap and sizes
422  int i,j,k;
423  int rows, cols, swap;
424  EST_TokenStream ts;
425  EST_read_status r;
426  EST_EstFileType t;
427  EST_Option hinfo;
428  bool ascii;
429 
430  if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
431  {
432  cerr << "DMatrix: can't open DMatrix input file "
433  << filename << endl;
434  return misc_read_error;
435  }
436  if ((r = read_est_header(ts, hinfo, ascii, t)) != format_ok)
437  return r;
438  if (t != est_file_dmatrix)
439  return misc_read_error;
440  if (hinfo.ival("version") != 1)
441  {
442  cerr << "DMatrix load: " << ts.pos_description() <<
443  " wrong version of DMatrix format expected 1 but found " <<
444  hinfo.ival("version") << endl;
445  return misc_read_error;
446  }
447  rows = hinfo.ival("rows");
448  cols = hinfo.ival("columns");
449  resize(rows,cols);
450 
451  if (ascii)
452  { // an ascii file
453  for (i = 0; i < num_rows(); ++i)
454  {
455  for (j = 0; j < num_columns(); ++j)
456  a_no_check(i,j) = atof(ts.get().string());
457  if (!ts.eoln())
458  {
459  cerr << "DMatrix load: " << ts.pos_description() <<
460  " missing end of line at end of row " << i << endl;
461  return misc_read_error;
462  }
463  }
464  }
465  else
466  { // a binary file
467  double *buff;
468  if ((EST_BIG_ENDIAN && (hinfo.sval("ByteOrder")=="LittleEndian")) ||
469  (EST_LITTLE_ENDIAN && (hinfo.sval("ByteOrder") == "BigEndian")))
470  swap = TRUE;
471  else
472  swap = FALSE;
473 
474  buff = walloc(double,rows*cols);
475  // A single read is *much* faster than multiple reads
476  if (ts.fread(buff,sizeof(double),rows*cols) != rows*cols)
477  {
478  cerr << "EST_DMatrix: binload: short file in \""
479  << filename << "\"" << endl;
480  return misc_read_error;
481  }
482  if (swap)
483  swap_bytes_double(buff,rows*cols);
484  for (k = i = 0; i < num_rows(); ++i)
485  for (j = 0; j < num_columns(); ++j)
486  a_no_check(i,j) = buff[k++];
487  wfree(buff);
488  }
489 
490  ts.close();
491  return read_ok;
492 }
493 
494 EST_read_status EST_DMatrix::load(const EST_String &filename)
495 {
496  EST_read_status r;
497 
498  if ((r = est_load(filename)) == format_ok)
499  return r;
500  else if (r == wrong_format)
501  { // maybe its an ancient ascii file
502  EST_TokenStream ts, tt;
503  EST_StrList sl;
504  int i, j, n_rows=0, n_cols=0;
505  EST_String t;
506  EST_Litem *p;
507  if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
508  {
509  cerr << "Can't open DMatrix file " << filename << endl;
510  return misc_read_error;
511  }
512  // set up the character constant values for this stream
513  ts.set_SingleCharSymbols(";");
514 
515  // first read in as list
516  for (n_rows = 0; !ts.eof(); ++n_rows)
517  sl.append(ts.get_upto_eoln().string());
518 
519  if (n_rows > 0)
520  {
521  tt.open_string(sl.first());
522  for (n_cols = 0; !tt.eof(); ++n_cols)
523  tt.get().string();
524  }
525 
526  // resize track and copy values in
527  resize(n_rows, n_cols);
528 
529  for (p = sl.head(), i = 0; p != 0; ++i, p = p->next())
530  {
531  tt.open_string(sl(p));
532  for (j = 0; !tt.eof(); ++j)
533  a_no_check(i,j) = atof(tt.get().string());
534  if (j != n_cols)
535  {
536  cerr << "Wrong number of points in row " << i << endl;
537  cerr << "Expected " << n_cols << " got " << j << endl;
538  return misc_read_error;
539  }
540  }
541  return format_ok;
542  }
543  else
544  return r;
545 
546  return format_ok;
547 }
548 
549 EST_read_status EST_DVector::est_load(const EST_String &filename)
550 {
551  // ascii/binary load with short header for byte swap and sizes
552  int i,k;
553  int l, swap;
554  EST_TokenStream ts;
555  EST_read_status r;
556  EST_EstFileType t;
557  EST_Option hinfo;
558  bool ascii;
559 
560  if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
561  {
562  cerr << "DVector: can't open DVector input file "
563  << filename << endl;
564  return misc_read_error;
565  }
566  if ((r = read_est_header(ts, hinfo, ascii, t)) != format_ok)
567  return r;
568  if (t != est_file_dvector)
569  return misc_read_error;
570  if (hinfo.ival("version") != 1)
571  {
572  cerr << "DVector load: " << ts.pos_description() <<
573  " wrong version of DVector format expected 1 but found " <<
574  hinfo.ival("version") << endl;
575  return misc_read_error;
576  }
577  l = hinfo.ival("length");
578  resize(l);
579 
580  if (ascii)
581  { // an ascii file
582  for (i = 0; i < length(); ++i)
583  a_no_check(i) = atof(ts.get().string());
584  }
585  else
586  { // a binary file
587  double *buff;
588  if ((EST_BIG_ENDIAN && (hinfo.sval("ByteOrder")=="LittleEndian")) ||
589  (EST_LITTLE_ENDIAN && (hinfo.sval("ByteOrder") == "BigEndian")))
590  swap = TRUE;
591  else
592  swap = FALSE;
593 
594  buff = walloc(double,l);
595  // A single read is *much* faster than multiple reads
596  if (ts.fread(buff,sizeof(double),l) != l)
597  {
598  cerr << "EST_DVector: binload: short file in \""
599  << filename << "\"" << endl;
600  return misc_read_error;
601  }
602  if (swap)
603  swap_bytes_double(buff,l);
604  for (k = i = 0; i < length(); ++i)
605  a_no_check(i) = buff[k++];
606  wfree(buff);
607  }
608 
609  ts.close();
610  return read_ok;
611 }
612 
613 EST_read_status EST_DVector::load(const EST_String &filename)
614 {
615 
616  EST_read_status r;
617 
618  if ((r = est_load(filename)) == format_ok)
619  return r;
620  else if (r == wrong_format)
621  { // maybe its an ancient ascii file
622  EST_TokenStream ts;
623  EST_String s;
624  int i;
625 
626  i = 0;
627 
628  if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
629  {
630  cerr << "can't open vector input file " << filename << endl;
631  return misc_read_error;
632  }
633  ts.set_SingleCharSymbols(";");
634 
635  while (!ts.eof())
636  {
637  ts.get();
638  ++i;
639  }
640  resize(i);
641 
642  ts.close();
643  if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
644  {
645  cerr << "can't open vector input file " << filename << endl;
646  return misc_read_error;
647  }
648 
649  for (i = 0; !ts.eof(); ++i)
650  {
651  s = ts.get().string();
652  (*this)[i] = atof(s); // actually returns double
653  }
654  ts.close();
655  return format_ok;
656  }
657  else
658  return r;
659 
660  return format_ok;
661 
662 }
663 
664 
666 {
667  int i;
668  if(n() != s.n()){
669  cerr << "Cannot elementwise add vectors of differing lengths"
670  << endl;
671  return *this;
672  }
673 
674  for (i = 0; i < n(); ++i)
675  (*this)[i] += s(i);
676 
677 
678  return *this;
679 }
680 
681 
683 {
684  if(n() != s.n()){
685  cerr << "Cannot elementwise multiply vectors of differing lengths"
686  << endl;
687  return *this;
688  }
689 
690  for (int i = 0; i < n(); ++i)
691  (*this)[i] *= s(i);
692 
693  return *this;
694 }
695 
697 {
698  for (int i = 0; i < n(); ++i)
699  (*this)[i] *= f;
700 
701  return *this;
702 }
703 
704 double operator*(const EST_DVector &v1, const EST_DVector &v2)
705 {
706  if(v1.length() != v2.length())
707  {
708  cerr << "Can't do vector dot prod - differing vector sizes !" << endl;
709  return 0;
710  }
711  double p=0;
712  for (int i = 0; i < v1.length(); ++i)
713  p += v1.a_no_check(i) * v2.a_no_check(i);
714 
715  return p;
716 }
717 
718 
720 {
721  for (int i = 0; i < n(); ++i)
722  (*this)[i] /= f;
723 
724  return *this;
725 }
726 
727 
728 EST_write_status EST_DVector::save(const EST_String &filename,
729  const EST_String &type)
730 {
731 
732  if ((type == "est_ascii") || (type == "est_binary"))
733  return est_save(filename,type);
734  else
735  { // the old stuff raw unheadered
736  int i;
737  ostream *outf;
738  if (filename == "-")
739  outf = &cout;
740  else
741  outf = new ofstream(filename);
742 
743  outf->precision(25);
744  if (!(*outf))
745  {
746  cerr << "DVector: can't open file \"" << filename
747  <<"\" for writing" << endl;
748  return misc_write_error;
749  }
750 
751  for (i = 0; i < length(); ++i)
752  *outf << a_no_check(i) << " ";
753  *outf << endl;
754 
755  if (outf != &cout)
756  delete outf;
757 
758  return write_ok;
759  }
760 }
761 
762 EST_write_status EST_DVector::est_save(const EST_String &filename,
763  const EST_String &type)
764 {
765  // Binary save with short header for byte swap and sizes
766  int i;
767  FILE *fd;
768  if (filename == "-")
769  fd = stdout;
770  else if ((fd = fopen(filename, "wb")) == NULL)
771  {
772  cerr << "EST_DVector: binsave: failed to open \"" << filename <<
773  "\" for writing" << endl;
774  return misc_write_error;
775  }
776 
777  fprintf(fd,"EST_File dvector\n");
778  fprintf(fd,"version 1\n");
779  if (type == "est_binary")
780  {
781  fprintf(fd,"DataType binary\n");
782  if (EST_LITTLE_ENDIAN)
783  fprintf(fd,"ByteOrder LittleEndian\n");
784  else
785  fprintf(fd,"ByteOrder BigEndian\n");
786  }
787  else
788  fprintf(fd,"DataType ascii\n");
789 
790  fprintf(fd,"length %d\n",length());
791  fprintf(fd,"EST_Header_End\n");
792 
793  if (type == "est_binary")
794  {
795  for (i = 0; i < length(); ++i)
796  if (fwrite(&a_no_check(i),sizeof(double),1,fd) != 1)
797  {
798  cerr << "EST_DVector: binsave: failed to write item "
799  << i << " to \"" << filename << "\"" << endl;
800  return misc_write_error;
801  }
802  }
803  else
804  { // est_ascii
805  for (i = 0; i < length(); ++i)
806  fprintf(fd,"%.25f ",a_no_check(i));
807  fprintf(fd,"\n");
808  }
809 
810  if (fd != stdout)
811  fclose(fd);
812 
813  return write_ok;
814 }
815