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 /*                         Author :  Paul Taylor                         */
34 /*                         Date   :  April 1995                          */
35 /*-----------------------------------------------------------------------*/
36 /*                        Matrix Class for floats                        */
37 /*                                                                       */
38 /*=======================================================================*/
39 
40 #include <cstdlib>
41 #include <cstdio>
42 #include <fstream>
43 #include <cmath>
44 #include <climits>
45 using namespace std;
46 
47 #include "EST_String.h"
48 #include "EST_types.h"
49 #include "EST_FileType.h"
50 #include "EST_Option.h"
51 #include "EST_FMatrix.h"
52 #include "EST_cutils.h"  // for swap functions
53 #include "EST_Token.h"
54 
55 
56 /* EST_FVector may used as EST_Val */
57 VAL_REGISTER_CLASS(fvector,EST_FVector)
58 
59 /* EST_FMatrix may used as EST_Val */
60 VAL_REGISTER_CLASS(fmatrix,EST_FMatrix)
61 
62 EST_String EST_FMatrix::default_file_type = "est_ascii";
63 
EST_FMatrix(const EST_FMatrix & a,int b)64 EST_FMatrix::EST_FMatrix(const EST_FMatrix &a, int b)
65 :EST_TSimpleMatrix<float>(a.num_rows(), a.num_columns())
66 {
67     float vv = 0.0;
68     if (b < 0)
69 	return;
70     if (b == 0)
71 	fill(vv);
72 }
73 
operator +=(const EST_FMatrix & a)74 EST_FMatrix & EST_FMatrix::operator+=(const EST_FMatrix &a)
75 {
76     int i, j;
77     if (a.num_columns() != num_columns())
78     {
79 	cerr <<"Matrix addition error: bad number of columns\n";
80 	return *this;
81     }
82     if (a.num_rows() != num_rows())
83     {
84 	cerr <<"Matrix addition error: bad number of rows\n";
85 	return *this;
86     }
87     for (i = 0; i < num_rows(); ++i)
88 	for (j = 0; j < num_columns(); ++j)
89 	    a_no_check(i, j) += a.a_no_check(i,j);
90 
91     return *this;
92 }
93 
operator -=(const EST_FMatrix & a)94 EST_FMatrix & EST_FMatrix::operator-=(const EST_FMatrix &a)
95 {
96     int i, j;
97     if (a.num_columns() != num_columns())
98     {
99 	cerr <<"Matrix subtraction error: bad number of columns\n";
100 	return *this;
101     }
102     if (a.num_rows() != num_rows())
103     {
104 	cerr <<"Matrix subtraction error: bad number of rows\n";
105 	return *this;
106     }
107     for (i = 0; i < num_rows(); ++i)
108 	for (j = 0; j < num_columns(); ++j)
109 	    a_no_check(i, j) -= a.a_no_check(i,j);
110 
111     return *this;
112 }
113 
operator *=(const float f)114 EST_FMatrix & EST_FMatrix::operator*=(const float f)
115 {
116 
117     int i,j;
118     for (i = 0; i < num_rows(); ++i)
119 	for (j = 0; j < num_columns(); ++j)
120 	    a_no_check(i, j) *= f;
121 
122     return *this;
123 }
124 
operator /=(const float f)125 EST_FMatrix & EST_FMatrix::operator/=(const float f)
126 {
127 
128     int i,j;
129     for (i = 0; i < num_rows(); ++i)
130 	for (j = 0; j < num_columns(); ++j)
131 	    a_no_check(i, j) /= f;
132 
133     return *this;
134 }
135 
operator +(const EST_FMatrix & a,const EST_FMatrix & b)136 EST_FMatrix operator+(const EST_FMatrix &a, const EST_FMatrix &b)
137 {
138     EST_FMatrix ab;
139     int i, j;
140     if (a.num_columns() != b.num_columns())
141     {
142 	cerr <<"Matrix addition error: bad number of columns\n";
143 	return ab;
144     }
145     if (a.num_rows() != b.num_rows())
146     {
147 	cerr <<"Matrix addition error: bad number of rows\n";
148 	return ab;
149     }
150     ab.resize(a.num_rows(), a.num_columns());
151     for (i = 0; i < a.num_rows(); ++i)
152 	for (j = 0; j < a.num_columns(); ++j)
153 	    ab.a_no_check(i, j) = a.a_no_check(i, j) + b.a_no_check(i, j);
154 
155     return ab;
156 }
157 
operator -(const EST_FMatrix & a,const EST_FMatrix & b)158 EST_FMatrix operator-(const EST_FMatrix &a,const EST_FMatrix &b)
159 {
160     EST_FMatrix ab;
161     int i, j;
162 
163     if (a.num_columns() != b.num_columns())
164     {
165 	cerr <<"Matrix subtraction error: bad number of columns:" <<
166 	    a.num_columns() << " and " << b.num_columns() << endl;
167 	return ab;
168     }
169     if (a.num_rows() != b.num_rows())
170     {
171 	cerr <<"Matrix subtraction error: bad number of rows\n";
172 	return ab;
173     }
174     ab.resize(a.num_rows(), a.num_columns());
175     for (i = 0; i < a.num_rows(); ++i)
176 	for (j = 0; j < a.num_columns(); ++j)
177 	    ab.a_no_check(i, j) = a.a_no_check(i, j) - b.a_no_check(i, j);
178 
179     return ab;
180 }
181 
operator *(const EST_FMatrix & a,const float x)182 EST_FMatrix operator*(const EST_FMatrix &a, const float x)
183 {
184     EST_FMatrix b(a, 0);
185     int i, j;
186 
187     for (i = 0; i < a.num_rows(); ++i)
188 	for (j = 0; j < a.num_columns(); ++j)
189 	    b.a_no_check(i,j) = a.a_no_check(i,j) * x;
190 
191     return b;
192 }
193 
operator !=(const EST_FVector & fv1,const EST_FVector & fv2)194 int operator !=(const EST_FVector &fv1,
195 		const EST_FVector &fv2)
196 {
197     int i;
198     if(fv1.length() != fv2.length())
199 	return FALSE;
200     for(i=0;i<fv1.length();i++)
201 	if(fv1.a_no_check(i) != fv2.a_no_check(i))
202 	    return FALSE;
203 
204     return TRUE;
205 }
206 
operator *(const EST_FMatrix & a,const EST_FVector & v)207 EST_FVector operator*(const EST_FMatrix &a, const EST_FVector &v)
208 {
209     // treat the vector as a column vector
210     // multiply each row of the matrix in turn by the vector
211 
212     EST_FVector b;
213     b.resize(a.num_rows());
214 
215     if(a.num_columns() != v.n())
216     {
217 	cerr <<"Matrix-vector multiplication error: matrix rows != vector size"
218 	     << endl;
219 	return b;
220     }
221 
222     int i, j;
223     for (i = 0; i < a.num_rows(); ++i){
224 	b[i] = 0.0;
225 	for (j = 0; j < a.num_columns(); ++j)
226 	    b.a_no_check(i) += a.a_no_check(i,j) * v.a_no_check(j);
227     }
228     return b;
229 }
230 
operator +(const EST_FVector & a,const EST_FVector & b)231 EST_FVector operator+(const EST_FVector &a, const EST_FVector &b)
232 {
233     EST_FVector ab;
234     int i;
235     if (a.length() != b.length())
236     {
237 	cerr <<"Vector addition error: mismatched lengths\n";
238 	return ab;
239     }
240 
241     ab.resize(a.length());
242     for (i = 0; i < a.length(); ++i)
243 	ab.a_no_check(i) = a.a_no_check(i) + b.a_no_check(i);
244 
245     return ab;
246 }
247 
operator -(const EST_FVector & a,const EST_FVector & b)248 EST_FVector operator-(const EST_FVector &a, const EST_FVector &b)
249 {
250     EST_FVector ab;
251     int i;
252     if (a.length() != b.length())
253     {
254 	cerr <<"Vector subtraction error: mismatched lengths\n";
255 	return ab;
256     }
257 
258     ab.resize(a.length());
259     for (i = 0; i < a.length(); ++i)
260 	ab.a_no_check(i) = a.a_no_check(i) - b.a_no_check(i);
261 
262     return ab;
263 }
264 
265 
operator *(const EST_FVector & v,const EST_FMatrix & a)266 EST_FVector operator*(const EST_FVector &v,const EST_FMatrix &a)
267 {
268     // treat the vector as a row vector
269     // multiply the vector by each column of the matrix in turn
270 
271     EST_FVector b;
272     b.resize(a.num_columns());
273 
274     if(a.num_columns() != v.n())
275     {
276 	cerr <<"Matrix-vector multiplication error: matrix rows != vector size"
277 	     << endl;
278 	return b;
279     }
280 
281     int i, j;
282     for (j = 0; j < a.num_columns(); ++j){
283 	b[j] = 0.0;
284 	for (i = 0; i < a.num_rows(); ++i)
285 	    b.a_no_check(i) += a.a_no_check(i,j) * v.a_no_check(j);
286     }
287     return b;
288 }
289 
290 
291 #if 0
292 EST_FMatrix operator/(const EST_FMatrix &a, float x)
293 {
294     return (a * (1/x));
295 }
296 #endif
297 
operator *(const EST_FMatrix & a,const EST_FMatrix & b)298 EST_FMatrix operator*(const EST_FMatrix &a, const EST_FMatrix &b)
299 {
300     EST_FMatrix ab;
301     multiply(a,b,ab);
302     return ab;
303 }
304 
multiply(const EST_FMatrix & a,const EST_FMatrix & b,EST_FMatrix & ab)305 void multiply(const EST_FMatrix &a, const EST_FMatrix &b, EST_FMatrix &ab)
306 {
307 
308     if (a.num_columns() != b.num_rows())
309     {
310 	cerr <<"Matrix multiply error: a.num_columns() != b.num_rows()\n";
311 	return;
312     }
313 
314     ab.resize(a.num_rows(), b.num_columns());
315     int i, j, k, n;
316     n = a.num_columns();	// could also be b.num_rows()
317 
318     for (i = 0; i < a.num_rows(); ++i)
319 	for (k = 0; k < b.num_columns(); ++k)
320 	{
321 	    ab.a_no_check(i, k) = 0.0;
322 	    for (j = 0; j < n; ++j)
323 		ab.a_no_check(i, k) +=
324 		    a.a_no_check(i, j) * b.a_no_check(j, k);
325 	}
326 }
327 
copyin(float ** inx,int rows,int cols)328 void EST_FMatrix::copyin(float **inx, int rows, int cols)
329 {
330     int i, j;
331 
332     resize(rows, cols);
333 
334     for (i = 0; i < rows; ++i)
335 	for (j = 0; j < cols; ++j)
336 	    a_no_check(i,j) = inx[i][j];
337 
338 }
339 
save(const EST_String & filename,const EST_String & type)340 EST_write_status EST_FMatrix::save(const EST_String &filename,
341 				   const EST_String &type)
342 {
343 
344     if ((type == "est_ascii") || (type == "est_binary"))
345 	return est_save(filename,type);
346     else
347     {   // the old stuff raw unheadered
348 	int i, j;
349 	ostream *outf;
350 	if (filename == "-")
351 	    outf = &cout;
352 	else
353 	    outf = new ofstream(filename);
354 
355 	if (!(*outf))
356 	{
357 	    cerr << "FMatrix: can't open file \"" << filename
358 		<<"\" for writing" << endl;
359 	    return misc_write_error;
360 	}
361 
362 	for (i = 0; i < num_rows(); ++i)
363 	{
364 	    for (j = 0; j < num_columns(); ++j)
365 		*outf << a_no_check(i,j) << " ";
366 	    *outf << endl;
367 	}
368 
369 	if (outf != &cout)
370 	    delete outf;
371 
372 	return write_ok;
373     }
374 }
375 
est_save(const EST_String & filename,const EST_String & type)376 EST_write_status EST_FMatrix::est_save(const EST_String &filename,
377 				       const EST_String &type)
378 {
379     // Binary save with short header for byte swap and sizes
380     int i,j;
381     FILE *fd;
382     if (filename == "-")
383 	fd = stdout;
384     else if ((fd = fopen(filename, "wb")) == NULL)
385     {
386 	cerr << "EST_FMatrix: binsave: failed to open \"" << filename <<
387 	    "\" for writing" << endl;
388 	return misc_write_error;
389     }
390 
391     fprintf(fd,"EST_File fmatrix\n");
392     fprintf(fd,"version 1\n");
393     if (type == "est_binary")
394     {
395 	fprintf(fd,"DataType binary\n");
396 	if (EST_LITTLE_ENDIAN)
397 	    fprintf(fd,"ByteOrder LittleEndian\n");
398 	else
399 	    fprintf(fd,"ByteOrder BigEndian\n");
400     }
401     else
402 	fprintf(fd,"DataType ascii\n");
403 
404     fprintf(fd,"rows %d\n",num_rows());
405     fprintf(fd,"columns %d\n",num_columns());
406 
407     fprintf(fd,"EST_Header_End\n");
408 
409     if (type == "est_binary")
410     {
411 	for (i = 0; i < num_rows(); ++i)
412 	    for (j=0; j < num_columns(); j++)
413 		if (fwrite(&a_no_check(i,j),sizeof(float),1,fd) != 1)
414 		{
415 		    cerr << "EST_FMatrix: binsave: failed to write row "
416 			<< i << " column " << j
417 			    << " to \"" << filename << "\"" << endl;
418 		    return misc_write_error;
419 		}
420     }
421     else
422     {   // est_ascii
423 	for (i = 0; i < num_rows(); ++i)
424 	{
425 	    for (j=0; j < num_columns(); j++)
426 		fprintf(fd,"%f ",a_no_check(i,j));
427 	    fprintf(fd,"\n");
428 	}
429     }
430 
431     if (fd != stdout)
432 	fclose(fd);
433 
434     return write_ok;
435 }
436 
est_load(const EST_String & filename)437 EST_read_status EST_FMatrix::est_load(const EST_String &filename)
438 {
439     // ascii/binary load with short header for byte swap and sizes
440     int i,j,k;
441     int rows, cols, swap;
442     EST_TokenStream ts;
443     EST_read_status r;
444     EST_EstFileType t;
445     EST_Option hinfo;
446     bool ascii;
447 
448     if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
449     {
450 	cerr << "FMatrix: can't open fmatrix input file "
451 	    << filename << endl;
452 	return misc_read_error;
453     }
454     if ((r = read_est_header(ts, hinfo, ascii, t)) != format_ok)
455 	return r;
456     if (t != est_file_fmatrix)
457 	return misc_read_error;
458     if (hinfo.ival("version") != 1)
459     {
460 	cerr << "FMatrix load: " << ts.pos_description() <<
461 	    " wrong version of fmatrix format expected 1 but found " <<
462 		hinfo.ival("version") << endl;
463 	return misc_read_error;
464     }
465     rows = hinfo.ival("rows");
466     cols = hinfo.ival("columns");
467     resize(rows,cols);
468 
469     if (ascii)
470     {   // an ascii file
471 	for (i = 0; i < num_rows(); ++i)
472 	{
473 	    for (j = 0; j < num_columns(); ++j)
474 		a_no_check(i,j) = atof(ts.get().string());
475 	    if (!ts.eoln())
476 	    {
477 		cerr << "FMatrix load: " << ts.pos_description() <<
478 		    " missing end of line at end of row " << i << endl;
479 		return misc_read_error;
480 	    }
481 	}
482     }
483     else
484     {   // a binary file
485 	float *buff;
486 	if ((EST_BIG_ENDIAN && (hinfo.sval("ByteOrder")=="LittleEndian")) ||
487 	    (EST_LITTLE_ENDIAN && (hinfo.sval("ByteOrder") == "BigEndian")))
488 	    swap = TRUE;
489 	else
490 	    swap = FALSE;
491 
492 	buff = walloc(float,rows*cols);
493 	// A single read is *much* faster than multiple reads
494 	if (ts.fread(buff,sizeof(float),rows*cols) != rows*cols)
495 	{
496 	    cerr << "EST_FMatrix: binload: short file in \""
497 		<< filename << "\"" << endl;
498 	    return misc_read_error;
499 	}
500 	if (swap)
501 	    swap_bytes_float(buff,rows*cols);
502 	for (k = i = 0; i < num_rows(); ++i)
503 	    for (j = 0; j < num_columns(); ++j)
504 		a_no_check(i,j) = buff[k++];
505 	wfree(buff);
506     }
507 
508     ts.close();
509 
510     return read_ok;
511 }
512 
load(const EST_String & filename)513 EST_read_status EST_FMatrix::load(const EST_String &filename)
514 {
515     EST_read_status r;
516 
517     if ((r = est_load(filename)) == format_ok)
518 	return r;
519     else if (r == wrong_format)
520     {   // maybe its an ancient ascii file
521 	EST_TokenStream ts, tt;
522 	EST_StrList sl;
523 	int i, j, n_rows=0, n_cols=0;
524 	EST_String t;
525 	EST_Litem *p;
526 	if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
527 	{
528 	    cerr << "Can't open fmatrix file " << filename << endl;
529 	    return misc_read_error;
530 	}
531 	// set up the character constant values for this stream
532 	ts.set_SingleCharSymbols(";");
533 
534 	// first read in as list
535 	for (n_rows = 0; !ts.eof(); ++n_rows)
536 	    sl.append(ts.get_upto_eoln().string());
537 
538 	if (n_rows > 0)
539 	{
540 	    tt.open_string(sl.first());
541 	    for (n_cols = 0; !tt.eof(); ++n_cols)
542 		tt.get().string();
543 	}
544 
545 	// resize track and copy values in
546 	resize(n_rows, n_cols);
547 
548 	for (p = sl.head(), i = 0; p != 0; ++i, p = p->next())
549 	{
550 	    tt.open_string(sl(p));
551 	    for (j = 0; !tt.eof(); ++j)
552 		a_no_check(i,j) = atof(tt.get().string());
553 	    if (j != n_cols)
554 	    {
555 		cerr << "Wrong number of points in row " << i << endl;
556 		cerr << "Expected " << n_cols << " got " << j << endl;
557 		return misc_read_error;
558 	    }
559 	}
560 	return format_ok;
561     }
562     else
563 	return r;
564 }
565 
566 
operator +=(const EST_FVector & s)567 EST_FVector & EST_FVector::operator+=(const EST_FVector &s)
568 {
569     int i;
570     if(n() != s.n()){
571 	cerr << "Cannot elementwise add vectors of differing lengths"
572 	    << endl;
573 	return *this;
574     }
575 
576     for (i = 0; i < n(); ++i)
577 	(*this)[i] += s(i);
578 
579 
580     return *this;
581 }
582 
583 
operator *=(const EST_FVector & s)584 EST_FVector& EST_FVector::operator*=(const EST_FVector &s)
585 {
586     if(n() != s.n()){
587 	cerr << "Cannot elementwise multiply vectors of differing lengths"
588 	    << endl;
589 	return *this;
590     }
591 
592     for (int i = 0; i < n(); ++i)
593 	(*this)[i] *= s(i);
594 
595     return *this;
596 }
597 
operator *=(const float f)598 EST_FVector& EST_FVector::operator*=(const float f)
599 {
600     for (int i = 0; i < n(); ++i)
601 	(*this)[i] *= f;
602 
603     return *this;
604 }
605 
606 
operator /=(const float f)607 EST_FVector& EST_FVector::operator/=(const float f)
608 {
609     for (int i = 0; i < n(); ++i)
610 	(*this)[i] /= f;
611 
612     return *this;
613 }
614 
615 
616 
est_load(const EST_String & filename)617 EST_read_status EST_FVector::est_load(const EST_String &filename)
618 {
619   // ascii/binary load with short header for byte swap and sizes
620   int i,k;
621   int l, swap;
622   EST_TokenStream ts;
623   EST_read_status r;
624   EST_EstFileType t;
625   EST_Option hinfo;
626   bool ascii;
627 
628   if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
629     {
630       cerr << "FVector: can't open FVector input file "
631 	   << filename << endl;
632       return misc_read_error;
633     }
634   if ((r = read_est_header(ts, hinfo, ascii, t)) != format_ok)
635     return r;
636     if (t != est_file_fvector)
637       return misc_read_error;
638     if (hinfo.ival("version") != 1)
639       {
640 	cerr << "FVector load: " << ts.pos_description() <<
641 	  " wrong version of FVector format expected 1 but found " <<
642 	  hinfo.ival("version") << endl;
643 	return misc_read_error;
644       }
645     l = hinfo.ival("length");
646     resize(l);
647 
648     if (ascii)
649       {   // an ascii file
650 	for (i = 0; i < length(); ++i)
651 	  a_no_check(i) = atof(ts.get().string());
652       }
653     else
654       {   // a binary file
655        float *buff;
656 	if ((EST_BIG_ENDIAN && (hinfo.sval("ByteOrder")=="LittleEndian")) ||
657 	    (EST_LITTLE_ENDIAN && (hinfo.sval("ByteOrder") == "BigEndian")))
658 	  swap = TRUE;
659 	else
660 	  swap = FALSE;
661 
662 	buff = walloc(float,l);
663 	// A single read is *much* faster than multiple reads
664 	if (ts.fread(buff,sizeof(float),l) != l)
665 	  {
666 	    cerr << "EST_FVector: binload: short file in \""
667 		 << filename << "\"" << endl;
668 	    return misc_read_error;
669 	  }
670 	if (swap)
671 	  swap_bytes_float(buff,l);
672 	for (k = i = 0; i < length(); ++i)
673 	  a_no_check(i) = buff[k++];
674 	wfree(buff);
675       }
676 
677     ts.close();
678     return read_ok;
679 }
680 
load(const EST_String & filename)681 EST_read_status EST_FVector::load(const EST_String &filename)
682 {
683 
684     EST_read_status r;
685 
686     if ((r = est_load(filename)) == format_ok)
687 	return r;
688     else if (r == wrong_format)
689     {   // maybe its an ancient ascii file
690       EST_TokenStream ts;
691       EST_String s;
692       int i;
693 
694       i = 0;
695 
696       if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
697 	{
698 	  cerr << "can't open vector input file " << filename << endl;
699 	  return misc_read_error;
700 	}
701       ts.set_SingleCharSymbols(";");
702 
703       while (!ts.eof())
704 	{
705 	  ts.get();
706 	  ++i;
707 	}
708       resize(i);
709 
710       ts.close();
711       if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
712 	{
713 	  cerr << "can't open vector input file " << filename << endl;
714 	  return misc_read_error;
715 	}
716 
717       for (i = 0; !ts.eof(); ++i)
718 	{
719 	  s = ts.get().string();
720 	  (*this)[i] = (float)(atof(s));  // actually returns double
721 	}
722       ts.close();
723       return format_ok;
724     }
725     else
726       return r;
727 
728     return format_ok;
729 
730 }
731 
732 /*
733 
734 EST_read_status EST_FVector::load(EST_String &filename)
735 {
736     EST_TokenStream ts;
737     EST_String s;
738     int i;
739 
740     i = 0;
741 
742     if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
743     {
744 	cerr << "can't open vector input file " << filename << endl;
745 	return misc_read_error;
746     }
747     ts.set_SingleCharSymbols(";");
748 
749     while (!ts.eof())
750     {
751 	ts.get();
752 	++i;
753     }
754     resize(i);
755 
756     ts.close();
757     if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
758     {
759 	cerr << "can't open vector input file " << filename << endl;
760 	return misc_read_error;
761     }
762 
763     for (i = 0; !ts.eof(); ++i)
764     {
765 	s = ts.get().string();
766 	(*this)[i] = atof(s);
767     }
768     ts.close();
769     return format_ok;
770 }
771 */
772 
773 //  EST_read_status EST_DVector::load(EST_String &filename)
774 //  {
775 //      EST_TokenStream ts;
776 //      EST_String s;
777 //      int i;
778 
779 //      i = 0;
780 
781 //      if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
782 //      {
783 //  	cerr << "can't open vector input file " << filename << endl;
784 //  	return misc_read_error;
785 //      }
786 //      ts.set_SingleCharSymbols(";");
787 
788 //      while (!ts.eof())
789 //      {
790 //  	ts.get();
791 //  	++i;
792 //      }
793 //      resize(i);
794 
795 //      ts.close();
796 //      if (((filename == "-") ? ts.open(cin) : ts.open(filename)) != 0)
797 //      {
798 //  	cerr << "can't open vector input file " << filename << endl;
799 //  	return misc_read_error;
800 //      }
801 
802 //      for (i = 0; !ts.eof(); ++i)
803 //      {
804 //  	s = ts.get().string();
805 //  	(*this)[i] = atof(s);  // actually returns double
806 //      }
807 //      ts.close();
808 //      return format_ok;
809 //  }
810 
811 
operator *(const EST_FVector & v1,const EST_FVector & v2)812 float operator*(const EST_FVector &v1, const EST_FVector &v2)
813 {
814     // dot product
815 
816     float b=0;
817 
818     if(v1.length() != v2.length())
819     {
820 	cerr <<"Vector dot product error: differing vector size"
821 	     << endl;
822 	return b;
823     }
824 
825     int i;
826     for (i = 0; i < v1.length(); ++i)
827 	b += v1.a_no_check(i) * v2.a_no_check(i);
828 
829     return b;
830 }
831 
832 
save(const EST_String & filename,const EST_String & type)833 EST_write_status EST_FVector::save(const EST_String &filename,
834 				   const EST_String &type)
835 {
836 
837     if ((type == "est_ascii") || (type == "est_binary"))
838 	return est_save(filename,type);
839     else
840     {   // the old stuff raw unheadered
841 	int i;
842 	ostream *outf;
843 	if (filename == "-")
844 	    outf = &cout;
845 	else
846 	    outf = new ofstream(filename);
847 
848 	outf->precision(25);
849 	if (!(*outf))
850 	{
851 	    cerr << "FVector: can't open file \"" << filename
852 		<<"\" for writing" << endl;
853 	    return misc_write_error;
854 	}
855 
856 	for (i = 0; i < length(); ++i)
857 	    *outf << a_no_check(i) << " ";
858 	*outf << endl;
859 
860 	if (outf != &cout)
861 	    delete outf;
862 
863 	return write_ok;
864     }
865 }
866 
est_save(const EST_String & filename,const EST_String & type)867 EST_write_status EST_FVector::est_save(const EST_String &filename,
868 				      const EST_String &type)
869 {
870     // Binary save with short header for byte swap and sizes
871     int i;
872     FILE *fd;
873     if (filename == "-")
874 	fd = stdout;
875     else if ((fd = fopen(filename, "wb")) == NULL)
876     {
877 	cerr << "EST_FVector: binsave: failed to open \"" << filename <<
878 	    "\" for writing" << endl;
879 	return misc_write_error;
880     }
881 
882     fprintf(fd,"EST_File fvector\n");
883     fprintf(fd,"version 1\n");
884     if (type == "est_binary")
885     {
886 	fprintf(fd,"DataType binary\n");
887 	if (EST_LITTLE_ENDIAN)
888 	    fprintf(fd,"ByteOrder LittleEndian\n");
889 	else
890 	    fprintf(fd,"ByteOrder BigEndian\n");
891     }
892     else
893 	fprintf(fd,"DataType ascii\n");
894 
895     fprintf(fd,"length %d\n",length());
896     fprintf(fd,"EST_Header_End\n");
897 
898     if (type == "est_binary")
899     {
900 	for (i = 0; i < length(); ++i)
901 	    if (fwrite(&a_no_check(i),sizeof(float),1,fd) != 1)
902 	    {
903 		cerr << "EST_FVector: binsave: failed to write item "
904 		     << i << " to \"" << filename << "\"" << endl;
905 		return misc_write_error;
906 	    }
907     }
908     else
909     {   // est_ascii
910 	for (i = 0; i < length(); ++i)
911 	    fprintf(fd,"%.25f ",a_no_check(i));
912 	fprintf(fd,"\n");
913     }
914 
915     if (fd != stdout)
916 	fclose(fd);
917 
918     return write_ok;
919 }
920 
921