« back to all changes in this revision
Viewing changes to lib/vb_vector.cpp
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Michael Hanke
- Date: 2010-06-06 11:33:11 UTC
- Revision ID: james.westby@ubuntu.com-20100606113311-v3c13imdkkd5n7ae
Tags: upstream-1.8.5~svn1172
Import upstream version 1.8.5~svn1172
- files added:
- brainBrowser
- client
- client/icons
- crunch
- db_import
- dbmisc
- dbmisc/networkdemo
- dbutil
- dist
- dist/elements
- dist/elements/filters
- dist/elements/noisemodels
- dist/elements/templates
- dist/etc
- dist/etc/jobtypes
- dist/etc/servers
- docs
- env
- fileformats
- gdscript
- getput
- lib
- munge
- qa
- qtglm
- qtvlsm
- qtvlsm/icons
- resample
- scheduler
- scripts
- scripts/examples
- server
- stand_alone
- stats
- tests
- tests/glm1
- utils
- vbsequence
- vbsequence/data
- vbsequence/kimberg-test
- vbsequence/kimberg-test/synmmn
- vbsequence/test-data
- vbview
- vbview/icons
- vbwidgets
added
removed
lib/vb_vector.cpp
1
2
// vb_vector.h
3
// VoxBo vector class
4
// Copyright (c) 1998-2010 by The VoxBo Development Team
5
6
// This file is part of VoxBo
7
//
8
// VoxBo is free software: you can redistribute it and/or modify it
9
// under the terms of the GNU General Public License as published by
10
// the Free Software Foundation, either version 3 of the License, or
11
// (at your option) any later version.
12
//
13
// VoxBo is distributed in the hope that it will be useful, but
14
// WITHOUT ANY WARRANTY; without even the implied warranty of
15
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16
// General Public License for more details.
17
//
18
// You should have received a copy of the GNU General Public License
19
// along with VoxBo. If not, see <http://www.gnu.org/licenses/>.
20
//
21
// For general information on VoxBo, including the latest complete
22
// source code and binary distributions, manual, and associated files,
23
// see the VoxBo home page at: http://www.voxbo.org/
24
//
25
// original version written by Kosh Banerjee
26
// changes/additions by Dan Kimberg
27
28
using namespace std;
29
30
#include "vbio.h"
31
#include <gsl/gsl_fit.h>
32
#include <gsl/gsl_statistics_double.h>
33
34
gsl_error_handler_t *VB_Vector::currentGSLErrorHandler = NULL;
35
36
VB_Vector::VB_Vector()
37
{
38
this->init(false, vb_double, "ref1");
39
this->theVector = NULL;
40
}
41
42
VB_Vector::VB_Vector(const size_t len)
43
{
44
this->init(false, vb_double, "ref1");
45
init(len);
46
}
47
48
void
49
VB_Vector::clear()
50
{
51
if (this->valid)
52
gsl_vector_free(this->theVector);
53
this->init(false, vb_double, "ref1");
54
this->theVector=NULL;
55
this->valid=false;
56
}
57
58
void
59
VB_Vector::print()
60
{
61
printf("vector:\n");
62
for (int i=0; i<size(); i++)
63
printf(" %010d: %.8f\n",i,getElement(i));
64
}
65
66
void
67
VB_Vector::init(const size_t len) throw (GenericExcep)
68
{
69
if (this->valid)
70
gsl_vector_free(this->theVector);
71
this->valid = false;
72
if (len>0) {
73
this->theVector=gsl_vector_calloc(len);
74
if (this->theVector != NULL)
75
this->valid = true;
76
}
77
}
78
79
void
80
VB_Vector::init(const bool validFlag,const VB_datatype dType,const string signature)
81
{
82
init(validFlag,dType,findFileFormat(signature));
83
}
84
85
void VB_Vector::init(const bool validFlag, const VB_datatype dType, const VBFF fType)
86
{
87
this->valid = validFlag;
88
this->dataType = dType;
89
this->fileFormat = fType;
90
}
91
92
/*********************************************************************
93
* This method allocates memory to the input gsl_matrix struct. *
94
*********************************************************************/
95
gsl_matrix *VB_Vector::initMatrix(const size_t rows, const size_t cols) const throw (GenericExcep)
96
{
97
98
/*********************************************************************
99
* Allocating the required memory to m. All elements of m will be *
100
* initialized to 0. *
101
*********************************************************************/
102
gsl_matrix *m = gsl_matrix_calloc(rows, cols);
103
104
/*********************************************************************
105
* If m is null, then we were unsuccessful in allocating memory for *
106
* it. Therefore, a GenericExcep is thrown. *
107
*********************************************************************/
108
if (m==NULL) {
109
110
/*********************************************************************
111
* errorMsg[] will hold an appropriate error message. *
112
*********************************************************************/
113
char errorMsg[OPT_STRING_LENGTH];
114
memset(errorMsg, 0, OPT_STRING_LENGTH);
115
116
/*********************************************************************
117
* Populating errorMsg[] with the appropriate error message. *
118
*********************************************************************/
119
sprintf(errorMsg, "The requested matrix size [%d x %d] could not be allocated.",
120
(int)rows, (int)cols);
121
122
/*********************************************************************
123
* Now throwing the exception. *
124
*********************************************************************/
125
throw GenericExcep(__LINE__, __FILE__, __FUNCTION__,
126
errorMsg);
127
128
} // if
129
130
/*********************************************************************
131
* Now returning the address of the allocated gsl_matrix struct. *
132
*********************************************************************/
133
return m;
134
135
} // gsl_matrix *VB_Vector::initMatrix(const size_t rows, const size_t cols) const throw (GenericExcep)
136
137
/*********************************************************************
138
* This private method is used as a wrapper around gsl_vector_memcpy()*
139
* to enable throwing an exception (if necessary). *
140
*********************************************************************/
141
void VB_Vector::GSLVectorMemcpy(gsl_vector *dest, const gsl_vector *src) const throw (GenericExcep)
142
{
143
144
/*********************************************************************
145
* gsl_vector_memcpy() is called to copy src to dest . It's return *
146
* value will be processed by VB_Vector::checkGSLStatus() and if the *
147
* return value is non-zero, then VB_Vector::checkGSLStatus() will *
148
* throw (and catch) an exception with an appropriate error message *
149
* (as returned by gsl_strerror()). *
150
*********************************************************************/
151
VB_Vector::checkGSLStatus(gsl_vector_memcpy(dest, src), __LINE__, __FILE__, __FUNCTION__);
152
153
} // void VB_Vector::GSLVectorMemcpy(gsl_vector *dest, const gsl_vector *src) const throw (GenericExcep)
154
155
/*********************************************************************
156
* This private method throws an exception if the input index is not *
157
* within the range of the length of this instance of VB_Vector. *
158
*********************************************************************/
159
void VB_Vector::checkVectorRange(const size_t index, const int lineNumber, const char *fileName, const char *prettyFunctionName) const throw (GenericExcep)
160
{
161
162
/*********************************************************************
163
* If the index is out of range of the vector length, then an *
164
* exception is thrown. *
165
*********************************************************************/
166
if (index >= this->getLength())
167
{
168
169
/*********************************************************************
170
* errorMsg[] will hold the error message that the index is out of *
171
* range. *
172
*********************************************************************/
173
char errorMsg[OPT_STRING_LENGTH];
174
memset(errorMsg, 0, OPT_STRING_LENGTH);
175
176
/*********************************************************************
177
* Populating errorMsg[] with the appropriate error message. *
178
*********************************************************************/
179
sprintf(errorMsg, "The index [%d] is not in the vector range [0, %d].",
180
(int)index, (int)(this->getLength() - 1));
181
182
/*********************************************************************
183
* Now throwing the exception. *
184
*********************************************************************/
185
throw GenericExcep(lineNumber, fileName, prettyFunctionName, errorMsg);
186
187
} // if
188
189
} // void VB_Vector::checkVectorRange(const size_t index, const int lineNumber, const char *fileName, const char *prettyFunctionName) const throw (GenericExcep)
190
191
/*********************************************************************
192
* This private static method throws an exception if the lengths of *
193
* the 2 input gsl_vector structs are not equal. *
194
*********************************************************************/
195
void VB_Vector::checkVectorLengths(const gsl_vector *V1, const gsl_vector *V2, const int lineNumber, const char *fileName, const char *prettyFunctionName) throw (GenericExcep)
196
{
197
198
/*********************************************************************
199
* We now check to see that we don't have a NULL gsl_vector. If either*
200
* gsl_vector argument is NULL, then an exception is thrown with an *
201
* appropriate error message. *
202
*********************************************************************/
203
if ( V1==NULL || V2==NULL )
204
{
205
ostringstream errorMsg;
206
errorMsg << "Have a NULL gsl_vector in [" << __FUNCTION__ << "].";
207
throw GenericExcep(lineNumber, fileName, prettyFunctionName, errorMsg.str());
208
} // if
209
210
/*********************************************************************
211
* If the vector lengths do no match, then an exception is thrown. *
212
*********************************************************************/
213
if (V1->size != V2->size)
214
{
215
216
/*********************************************************************
217
* errorMsg[] will hold the error message that the vector lengths are *
218
* not equal. *
219
*********************************************************************/
220
char errorMsg[OPT_STRING_LENGTH];
221
memset(errorMsg, 0, OPT_STRING_LENGTH);
222
223
/*********************************************************************
224
* Populating errorMsg[] with the appropriate error message. *
225
*********************************************************************/
226
sprintf(errorMsg, "Unequal vector lengths: [%d] and [%d]",
227
(int)(V1->size), (int)(V2->size));
228
229
/*********************************************************************
230
* Now throwing the exception. *
231
*********************************************************************/
232
throw GenericExcep(lineNumber, fileName, prettyFunctionName, errorMsg);
233
234
} // if
235
236
} // void VB_Vector::checkVectorLengths(const gsl_vector *V1, const gsl_vector *V2, const int lineNumber, const char *fileName, const char *prettyFunctionName) throw (GenericExcep)
237
238
/*********************************************************************
239
* This private static method throws an exception if the two input *
240
* size_t variables, len1 and len2, are not equal. *
241
*********************************************************************/
242
void VB_Vector::checkVectorLengths(const size_t len1, const size_t len2, const int lineNumber, const char *fileName, const char *prettyFunctionName) throw (GenericExcep)
243
{
244
245
/*********************************************************************
246
* If the vectors length do no match, then an exception is thrown. *
247
*********************************************************************/
248
if (len1 != len2)
249
{
250
251
/*********************************************************************
252
* errorMsg[] will hold the error message that the vector lengths are *
253
* not equal. *
254
*********************************************************************/
255
char errorMsg[OPT_STRING_LENGTH];
256
memset(errorMsg, 0, OPT_STRING_LENGTH);
257
258
/*********************************************************************
259
* Populating errorMsg[] with the appropriate error message. *
260
*********************************************************************/
261
sprintf(errorMsg, "Unequal vector lengths: [%d] and [%d]",
262
(int)len1, (int)len2);
263
264
/*********************************************************************
265
* Now throwing the exception. *
266
*********************************************************************/
267
throw GenericExcep(lineNumber, fileName, prettyFunctionName, errorMsg);
268
269
} // if
270
271
} // void VB_Vector::checkVectorLengths(const size_t len1, const size_t len2, const int lineNumber, const char *fileName, const char *prettyFunctionName) throw (GenericExcep)
272
273
/*********************************************************************
274
* This private static method throws a GenericExcep if the input *
275
* gsl_vector struct is null. *
276
*********************************************************************/
277
void VB_Vector::vectorNull(const gsl_vector *v) throw (GenericExcep)
278
{
279
280
/*********************************************************************
281
* If the input gsl_vector struct is null, then the exception is *
282
* thrown. *
283
*********************************************************************/
284
if (v==NULL)
285
{
286
throw GenericExcep(__LINE__, __FILE__, __FUNCTION__,
287
"ERROR: Unable to allocate memory for VB_Vector.");
288
} // if
289
290
} // void VB_Vector::vectorNull(const gsl_vector *v) throw (GenericExcep)
291
292
/*********************************************************************
293
* This private static method throws a GenericExcep and catches it *
294
* if the input variable status is non-zero. This method should be *
295
* called after a GSL function returns its status. *
296
*********************************************************************/
297
void VB_Vector::checkGSLStatus(const int status, const int lineNumber, const char *fileName, const char *prettyFunctionName)
298
{
299
300
/*********************************************************************
301
* If status is non-zero, then try/catch blocks are used to throw a *
302
* GenericExcep and then process it. *
303
*********************************************************************/
304
if (status)
305
{
306
try
307
{
308
throw GenericExcep(lineNumber, fileName, prettyFunctionName,
309
gsl_strerror(status));
310
} // try
311
catch (GenericExcep &e)
312
{
313
e.whatNoExit(lineNumber, fileName, prettyFunctionName);
314
} // catch
315
} // if
316
317
} // void VB_Vector::checkGSLStatus(const int status, const int lineNumber, const char *fileName, const char *prettyFunctionName)
318
319
/*********************************************************************
320
* This private static method throws a GenericExcep and catches it *
321
* if the input variable norm is zero. *
322
*********************************************************************/
323
void VB_Vector::zeroVector(const double norm, const int lineNumber, const char *fileName, const char *prettyFunctionName)
324
{
325
326
/*********************************************************************
327
* If norm is zero, then try/catch blocks are used to throw a *
328
* GenericExcep and then process it. *
329
*********************************************************************/
330
if (norm == 0.0)
331
{
332
try
333
{
334
throw GenericExcep(lineNumber, fileName, prettyFunctionName,
335
"The vector is a zero vector.");
336
} // try
337
catch (GenericExcep &e)
338
{
339
e.whatNoExit(lineNumber, fileName, prettyFunctionName);
340
} // catch
341
} // if
342
343
} // void VB_Vector::zeroVector(const double norm, const int lineNumber, const char *fileName, const char *prettyFunctionName)
344
345
/*********************************************************************
346
* This static method throws and catches an exception of type *
347
* GenericExcep. *
348
*********************************************************************/
349
void VB_Vector::createException(const char *errorMsg, const int lineNumber, const char *fileName, const char *prettyFunctionName)
350
{
351
352
/*********************************************************************
353
* Now throwing and catching the exception. *
354
*********************************************************************/
355
try
356
{
357
throw GenericExcep(lineNumber, fileName, prettyFunctionName, errorMsg);
358
} // try
359
catch (GenericExcep &e)
360
{
361
e.whatNoExit(lineNumber, fileName, prettyFunctionName);
362
} // catch
363
364
} // void VB_Vector::createException(const char *errorMsg, const int lineNumber, const char *fileName, const char *prettyFunctionName)
365
366
/*********************************************************************
367
* This static method throws and catches an exception of type *
368
* GenericExcep. *
369
*********************************************************************/
370
void VB_Vector::createException(const string& errorMsg, const int lineNumber, const string& fileName, const string& prettyFunctionName)
371
{
372
373
/*********************************************************************
374
* Now throwing the exception. *
375
*********************************************************************/
376
VB_Vector::createException(errorMsg.c_str(), lineNumber, fileName.c_str(), prettyFunctionName.c_str());
377
378
} // void VB_Vector::createException(const string& errorMsg, const int lineNumber, const string& fileName, const string& prettyFunctionName)
379
380
/*********************************************************************
381
* This static method will call VB_Vector::createException() if any *
382
* elements of the input gsl_struct are Inf or Nan. *
383
*********************************************************************/
384
void VB_Vector::checkFiniteness(const gsl_vector *v, const int lineNumber, const char *fileName, const char *prettyFunctionName)
385
{
386
387
/*********************************************************************
388
* The following for loop is used to examine each element of the *
389
* input gsl_vector struct. *
390
*********************************************************************/
391
for (size_t i = 0; i < v->size; i++)
392
{
393
394
/*********************************************************************
395
* If the input double to gsl_finite() is an Inf or a Nan, then *
396
* gsl_finite() will return 0. Otherwise, 1 is returned. *
397
*********************************************************************/
398
if (!gsl_finite(v->data[i]))
399
{
400
401
/*********************************************************************
402
* errorMsg[] will hold the error message that an Inf or Nan was *
403
* encountered. *
404
*********************************************************************/
405
char errorMsg[OPT_STRING_LENGTH];
406
memset(errorMsg, 0, OPT_STRING_LENGTH);
407
408
/*********************************************************************
409
* Populating errorMsg[] with the appropriate error message. *
410
*********************************************************************/
411
sprintf(errorMsg, "The vector element at index [%d] is an Inf or a NaN.", (int)i);
412
413
/*********************************************************************
414
* Now VB_Vector::createException() is called to trhow and catch a *
415
* GenericExcep so that an appropriate error message gets displayed. *
416
*********************************************************************/
417
VB_Vector::createException(errorMsg, lineNumber, fileName, prettyFunctionName);
418
419
/*********************************************************************
420
* We now break from the for loop since we have found at least 1 Inf *
421
* or Nan. *
422
*********************************************************************/
423
break;
424
425
} // if
426
427
} // for i
428
429
} // void VB_Vector::checkFiniteness(const gsl_vector *v, const int lineNumber, const char *fileName, const char *prettyFunctionName)
430
431
/*********************************************************************
432
* This constructor uses the input variable len as the length of the *
433
* the vector and the input variable data as the elements of the *
434
* vector. NOTE: Is is expected that data is an array of doubles whose*
435
* length is len. *
436
*********************************************************************/
437
VB_Vector::VB_Vector(const double *data, const size_t len)
438
{
439
440
/*********************************************************************
441
* Setting this->valid, this->dataType, and this->fileFormat. *
442
*********************************************************************/
443
this->init(false, vb_double, "ref1");
444
445
/*********************************************************************
446
* We call this->init() and catch its exception, if one is thrown. *
447
*********************************************************************/
448
try
449
{
450
this->init(len);
451
} // try
452
catch (GenericExcep &e)
453
{
454
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
455
} // catch
456
457
/*********************************************************************
458
* Now copying the values from data to this->theVector->data. *
459
*********************************************************************/
460
memcpy(this->theVector->data, data, this->theVector->size * sizeof(double));
461
462
} // VB_Vector::VB_Vector(const double *data, const size_t len)
463
464
/*********************************************************************
465
* This constructor builds a VB_Vector from the input gsl_vector *
466
* struct. *
467
*********************************************************************/
468
VB_Vector::VB_Vector(const gsl_vector *V2)
469
{
470
471
/*********************************************************************
472
* Setting this->fileName, this->valid, this->dataType, and *
473
* this->fileFormat. Since we don't have an actual file name, a *
474
* temporary file name is used. *
475
*********************************************************************/
476
char tmpFileName[STRINGLEN];
477
memset(tmpFileName, 0, STRINGLEN);
478
strcpy(tmpFileName, "./tmp-");
479
this->init(false, vb_double, "ref1");
480
481
/*********************************************************************
482
* Setting this->theVector to NULL. *
483
*********************************************************************/
484
this->theVector = NULL;
485
486
/*********************************************************************
487
* We call this->init() and catch its exception, if one is thrown. *
488
*********************************************************************/
489
try
490
{
491
492
/*********************************************************************
493
* Allocating a gsl_vector struct and all the elements will be set to *
494
* 0. *
495
*********************************************************************/
496
this->init(V2->size);
497
498
} // try
499
catch (GenericExcep &e)
500
{
501
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
502
} // catch
503
504
/*********************************************************************
505
* Now copying the values from V2.data to this->theVector->data. The *
506
* try/catch blocks are used to process the exception, if one is *
507
* thrown. *
508
*********************************************************************/
509
try
510
{
511
this->GSLVectorMemcpy(this->theVector, V2);
512
} // try
513
catch (GenericExcep &e)
514
{
515
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
516
} // catch
517
518
} // VB_Vector::VB_Vector(const gsl_vector *V2)
519
520
/*********************************************************************
521
* This constructor builds a VB_Vector from the input gsl_vector *
522
* struct. *
523
*********************************************************************/
524
VB_Vector::VB_Vector(const gsl_vector& V2)
525
{
526
527
/*********************************************************************
528
* Setting this->valid, this->dataType, and this->fileFormat. *
529
*********************************************************************/
530
this->init(false, vb_double, "ref1");
531
532
/*********************************************************************
533
* Setting this->theVector to NULL. *
534
*********************************************************************/
535
this->theVector = NULL;
536
537
/*********************************************************************
538
* We call this->init() and catch its exception, if one is thrown. *
539
*********************************************************************/
540
try
541
{
542
543
/*********************************************************************
544
* Allocating a gsl_vector struct and all the elements will be set to *
545
* 0. *
546
*********************************************************************/
547
this->init(V2.size);
548
} // try
549
catch (GenericExcep &e)
550
{
551
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
552
} // catch
553
554
/*********************************************************************
555
* Now copying the values from V2.data to this->theVector->data. The *
556
* try/catch blocks are used to process the exception, if one is *
557
* thrown. *
558
*********************************************************************/
559
try
560
{
561
this->GSLVectorMemcpy(this->theVector, &V2);
562
} // try
563
catch (GenericExcep &e)
564
{
565
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
566
} // catch
567
568
} // VB_Vector::VB_Vector(const gsl_vector& V2)
569
570
/*********************************************************************
571
* This is the copy constructor. *
572
*********************************************************************/
573
VB_Vector::VB_Vector(const VB_Vector& V2)
574
{
575
576
/*********************************************************************
577
* Setting this->valid, this->dataType, and this->fileFormat. *
578
*********************************************************************/
579
this->init(false, V2.dataType,V2.fileFormat);
580
581
/*********************************************************************
582
* Setting this->fileName to V2.fileName. *
583
*********************************************************************/
584
this->fileName = V2.fileName;
585
586
/*********************************************************************
587
* We call this->init() and catch its exception, if one is thrown. *
588
*********************************************************************/
589
try
590
{
591
592
/*********************************************************************
593
* If V2.theVector is not null, then a gsl_vector struct is allocated *
594
* and all the elements of this->theVector->data[] will be set to 0. *
595
* Otherwise, this->theVector is set to NULL. *
596
*********************************************************************/
597
if (V2.theVector)
598
{
599
this->init(V2.theVector->size);
600
601
/*********************************************************************
602
* Now copying the values from V2.theVector to this->theVector. The *
603
* try/catch blocks are used to process the exception, if one is *
604
* thrown. *
605
*********************************************************************/
606
try
607
{
608
this->GSLVectorMemcpy(this->theVector, V2.theVector);
609
} // try
610
catch (GenericExcep &e)
611
{
612
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
613
} // catch
614
615
} // if
616
else
617
{
618
this->theVector = NULL;
619
} // else
620
621
} // try
622
catch (GenericExcep &e)
623
{
624
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
625
} // catch
626
627
} // VB_Vector::VB_Vector(const VB_Vector& V2)
628
629
/*********************************************************************
630
* This constructor instantiates a VB_Vector from the input VBVector *
631
* pointer. *
632
*********************************************************************/
633
VB_Vector::VB_Vector(const VB_Vector *V)
634
{
635
636
/*********************************************************************
637
* Setting this->valid, this->dataType, and this->fileFormat. *
638
*********************************************************************/
639
this->init(false, V->dataType,V->fileFormat);
640
641
/*********************************************************************
642
* Setting this->fileName to V->fileName. *
643
*********************************************************************/
644
this->fileName = V->fileName;
645
646
/*********************************************************************
647
* We call this->init() and catch its exception, if one is thrown. *
648
*********************************************************************/
649
try
650
{
651
652
/*********************************************************************
653
* If V->theVector is not null, then a gsl_vector struct is allocated *
654
* and all the elements of this->theVector->data[] will be set to 0. *
655
* Otherwise, this->theVector is set to NULL. *
656
*********************************************************************/
657
if (V->theVector)
658
{
659
this->init(V->theVector->size);
660
} // if
661
else
662
{
663
this->theVector = NULL;
664
} // else
665
666
} // try
667
catch (GenericExcep &e)
668
{
669
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
670
} // catch
671
672
/*********************************************************************
673
* Now copying the values from V->theVector to this->theVector. The *
674
* try/catch blocks are used to process the exception, if one is *
675
* thrown. *
676
*********************************************************************/
677
try
678
{
679
this->GSLVectorMemcpy(this->theVector, V->theVector);
680
} // try
681
catch (GenericExcep &e)
682
{
683
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
684
} // catch
685
686
} // VB_Vector::VB_Vector(const VB_Vector *V)
687
688
/*********************************************************************
689
* Constructor to read in a vector from a VoxBo vector file. *
690
*********************************************************************/
691
VB_Vector::VB_Vector(const string& vecFile)
692
{
693
694
/*********************************************************************
695
* Setting this->valid, this->dataType, and this->fileFormat. *
696
*********************************************************************/
697
this->init(false, vb_double, "ref1");
698
699
/*********************************************************************
700
* Reading the VoxBo vector file. *
701
*********************************************************************/
702
if (this->ReadFile(vecFile))
703
{
704
ostringstream errorMsg;
705
errorMsg << "[" << __FUNCTION__ << "]: Unable to read the file [" << vecFile << "].";
706
printErrorMsg(VB_WARNING, errorMsg.str());
707
} // if
708
709
} // VB_Vector::VB_Vector(const string& vecFile)
710
711
/*********************************************************************
712
* Constructor to read in a vector from a VoxBo vector file. *
713
*********************************************************************/
714
VB_Vector::VB_Vector(const char *vecFile)
715
{
716
717
/*********************************************************************
718
* Setting this->valid, this->dataType, and this->fileFormat. *
719
*********************************************************************/
720
this->init(false, vb_double, "ref1");
721
722
/*********************************************************************
723
* Setting this->fileName. *
724
*********************************************************************/
725
this->fileName = vecFile;
726
727
/*********************************************************************
728
* Reading the VoxBo vector file. *
729
*********************************************************************/
730
if (this->ReadFile(vecFile))
731
{
732
ostringstream errorMsg;
733
errorMsg << "[" << __FUNCTION__ << "]: Unable to read the file [" << vecFile << "].";
734
printErrorMsg(VB_WARNING, errorMsg.str());
735
} // if
736
737
} // VB_Vector::VB_Vector(const char *vecFile)
738
739
/*********************************************************************
740
* This constructor creates an instance of VB_Vector from a *
741
* reference to a vector <double>. *
742
*********************************************************************/
743
VB_Vector::VB_Vector(const vector<double>& theVector)
744
{
745
746
/*********************************************************************
747
* Now setting this->valid (to false), this->dataType, and *
748
* this->fileFormat. Also, this->fileName will be set to temporary file *
749
* name. *
750
*********************************************************************/
751
this->init(false, vb_double, "ref1");
752
753
/*********************************************************************
754
* We call this->init() and catch its exception, if one is thrown. *
755
*********************************************************************/
756
try
757
{
758
this->init(theVector.size());
759
} // try
760
catch (GenericExcep &e)
761
{
762
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
763
} // catch
764
765
/*********************************************************************
766
* Now the elements from theVector are copied to this->theVector. *
767
*********************************************************************/
768
copy(theVector.begin(), theVector.end(), this->theVector->data);
769
770
} // VB_Vector::VB_Vector(const vector<double>& theVector)
771
772
/*********************************************************************
773
* This constructor creates an instance of VB_Vector from a *
774
* pointer to a vector <double>. *
775
*********************************************************************/
776
VB_Vector::VB_Vector(const vector<double> *theVector)
777
{
778
779
/*********************************************************************
780
* Now setting this->valid (to false), this->dataType, and *
781
* this->fileFormat. Also, this->fileName will be set to temporary *
782
* file name. *
783
*********************************************************************/
784
this->init(false, vb_double, "ref1");
785
786
/*********************************************************************
787
* We call this->init() and catch its exception, if one is thrown. *
788
*********************************************************************/
789
try
790
{
791
this->init(theVector->size());
792
} // try
793
catch (GenericExcep &e)
794
{
795
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
796
} // catch
797
798
/*********************************************************************
799
* Now the elements from theVector are copied to this->theVector. *
800
*********************************************************************/
801
copy(theVector->begin(), theVector->end(), this->theVector->data);
802
803
} // VB_Vector::VB_Vector(const vector<double> *theVector)
804
805
/*********************************************************************
806
* This constructor builds a VB_Vector from a pointer to a Vec. *
807
*********************************************************************/
808
VB_Vector::VB_Vector(const Vec *theVector)
809
{
810
811
/*********************************************************************
812
* Now setting this->valid (to false), this->dataType, and *
813
* this->fileFormat. Also, this->fileName will be set to temporary *
814
* file name. *
815
*********************************************************************/
816
this->init(false, vb_double, "ref1");
817
818
/*********************************************************************
819
* We call this->init() and catch its exception, if one is thrown. *
820
*********************************************************************/
821
try
822
{
823
this->init(theVector->size());
824
} // try
825
catch (GenericExcep &e)
826
{
827
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
828
} // catch
829
830
/*********************************************************************
831
* Now the elements from theVector are copied to this->theVector. *
832
*********************************************************************/
833
memcpy(this->theVector->data, theVector->data, sizeof(double) * theVector->size());
834
835
} // VB_Vector::VB_Vector(const Vec *theVector)
836
837
/*********************************************************************
838
* This constructor builds a VB_Vector from a reference to a Vec. *
839
*********************************************************************/
840
VB_Vector::VB_Vector(const Vec& theVector)
841
{
842
843
/*********************************************************************
844
* Now setting this->valid (to false), this->dataType, and *
845
* this->fileFormat. Also, this->fileName will be set to temporary *
846
* file name. *
847
*********************************************************************/
848
this->init(false, vb_double, "ref1");
849
850
/*********************************************************************
851
* We call this->init() and catch its exception, if one is thrown. *
852
*********************************************************************/
853
try
854
{
855
this->init(theVector.size());
856
} // try
857
catch (GenericExcep &e)
858
{
859
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
860
} // catch
861
862
/*********************************************************************
863
* Now the elements from theVector are copied to this->theVector. *
864
*********************************************************************/
865
memcpy(this->theVector->data, theVector.data, sizeof(double) * theVector.size());
866
867
} // VB_Vector::VB_Vector(const Vec& theVector)
868
869
/*********************************************************************
870
* This constructor is used to quickly generate a VB_Vector object *
871
* from the time series stored at the specified index. *
872
*********************************************************************/
873
VB_Vector::VB_Vector(const Tes& theTes, const unsigned long tSeriesIndex)
874
{
875
876
/*********************************************************************
877
* Setting this->valid, this->dataType, and this->fileFormat. *
878
*********************************************************************/
879
this->init(false, vb_double, "ref1");
880
881
/*********************************************************************
882
* We call this->init() and catch its exception, if one is thrown. *
883
*********************************************************************/
884
try
885
{
886
this->init(theTes.dimt);
887
} // try
888
catch (GenericExcep &e)
889
{
890
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
891
} // catch
892
893
/*********************************************************************
894
* If theTes.data[tSeriesIndex] is not NULL, then it holds a *
895
* non-zero time series. NOTE: If in fact theTes.data[tSeriesIndex] *
896
* is NULL, then we have a time series that is entirely zero. *
897
* Therefore, we do not need to do anything else since the above call *
898
* to this->init() sets all the elements to 0. *
899
*********************************************************************/
900
if (theTes.data[tSeriesIndex])
901
{
902
903
/*********************************************************************
904
* We now switch on the possible data types for theTes. *
905
*********************************************************************/
906
switch(theTes.datatype)
907
{
908
case vb_byte:
909
910
/*********************************************************************
911
* The following for loop assigns each element of this instance of *
912
* VB_Vector. *
913
*********************************************************************/
914
for (long i = 0; i < theTes.dimt; i++)
915
{
916
this->theVector->data[i] = *((unsigned char *) (theTes.data[tSeriesIndex] + (i * theTes.datasize)));
917
} // for i
918
919
break;
920
921
case vb_short:
922
923
/*********************************************************************
924
* The following for loop assigns each element of this instance of *
925
* VB_Vector. *
926
*********************************************************************/
927
for (long i = 0; i < theTes.dimt; i++)
928
{
929
this->theVector->data[i] = *((int16 *) (theTes.data[tSeriesIndex] + (i * theTes.datasize)));
930
} // for i
931
932
break;
933
934
case vb_long:
935
936
/*********************************************************************
937
* The following for loop assigns each element of this instance of *
938
* VB_Vector. *
939
*********************************************************************/
940
for (long i = 0; i < theTes.dimt; i++)
941
{
942
this->theVector->data[i] = *((int32 *) (theTes.data[tSeriesIndex] + (i * theTes.datasize)));
943
} // for i
944
945
break;
946
947
case vb_float:
948
949
/*********************************************************************
950
* The following for loop assigns each element of this instance of *
951
* VB_Vector. *
952
*********************************************************************/
953
for (long i = 0; i < theTes.dimt; i++)
954
{
955
this->theVector->data[i] = *((float *) (theTes.data[tSeriesIndex] + (i * theTes.datasize)));
956
} // for i
957
958
break;
959
960
case vb_double:
961
962
/*********************************************************************
963
* The following for loop assigns each element of this instance of *
964
* VB_Vector. *
965
*********************************************************************/
966
for (long i = 0; i < theTes.dimt; i++)
967
{
968
this->theVector->data[i] = *((double *) (theTes.data[tSeriesIndex] + (i * theTes.datasize)));
969
} // for i
970
971
break;
972
973
} // switch
974
975
} // if
976
977
} // VB_Vector::VB_Vector(const Tes& theTes, const unsigned long tSeriesIndex)
978
979
/*********************************************************************
980
* Method to turn off the current GSL error handler. *
981
*********************************************************************/
982
void VB_Vector::turnOffGSLErrorHandler() const
983
{
984
VB_Vector::currentGSLErrorHandler = gsl_set_error_handler_off();
985
} // void VB_Vector::turnOffGSLErrorHandler() const
986
987
/*********************************************************************
988
* Method to restore the GSL error handler to *
989
* VB_Vector::->currentGSLErrorHandler. *
990
*********************************************************************/
991
void VB_Vector::restoreGSLErrorHandler() const
992
{
993
994
/*********************************************************************
995
* If VB_Vector::currentGSLErrorHandler is not null, then the GSL *
996
* error handler is restored to VB_Vector::currentGSLErrorHandler. *
997
*********************************************************************/
998
if (VB_Vector::currentGSLErrorHandler)
999
{
1000
gsl_set_error_handler(VB_Vector::currentGSLErrorHandler);
1001
} // if
1002
1003
} // void VB_Vector::restoreGSLErrorHandler() const
1004
1005
/*********************************************************************
1006
* The destructor. *
1007
*********************************************************************/
1008
VB_Vector::~VB_Vector()
1009
{
1010
1011
/*********************************************************************
1012
* If this->valid is true, then a gsl_vector struct was allocated to *
1013
* this->theVector. Therefore, we delete the memory. *
1014
*********************************************************************/
1015
if (this->valid)
1016
{
1017
gsl_vector_free(this->theVector);
1018
} // if
1019
1020
} // VB_Vector::~VB_Vector()
1021
1022
/*********************************************************************
1023
* This method simply returns the specified element from the vector. *
1024
*********************************************************************/
1025
double VB_Vector::getElement(const size_t index) const
1026
{
1027
1028
/*********************************************************************
1029
* We now check to make sure that index is within *
1030
* [0, this->theVector->size - 1]. The try/catch blocks are used to *
1031
* process any exception, if one is thrown by *
1032
* this->checkVectorRange(). *
1033
*********************************************************************/
1034
try
1035
{
1036
this->checkVectorRange(index, __LINE__, __FILE__, __FUNCTION__);
1037
} // try
1038
catch (GenericExcep &e)
1039
{
1040
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1041
} // catch
1042
1043
/*********************************************************************
1044
* Now returning the specified element value. *
1045
*********************************************************************/
1046
return gsl_vector_get(this->theVector, index);
1047
1048
} // double VB_Vector::getElement(const size_t index) const
1049
1050
/*********************************************************************
1051
* This method simply returns a pointer to the specified element. *
1052
* NOTE: This method makes it possible to change the value of the *
1053
* specified element by using the return pointer, i.e., a different *
1054
* value can be assigned to the pointer. *
1055
*********************************************************************/
1056
double *VB_Vector::getElementPtr(size_t index) const
1057
{
1058
1059
/*********************************************************************
1060
* We now check to make sure that index is within *
1061
* [0, this->theVector->size - 1]. The try/catch blocks are used to *
1062
* process any exception, if one is thrown by *
1063
* this->checkVectorRange(). *
1064
*********************************************************************/
1065
try
1066
{
1067
this->checkVectorRange(index, __LINE__, __FILE__, __FUNCTION__);
1068
} // try
1069
catch (GenericExcep &e)
1070
{
1071
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1072
} // catch
1073
1074
/*********************************************************************
1075
* Calling gsl_vector_ptr() to return the pointer to the specified *
1076
* vector element. *
1077
*********************************************************************/
1078
return gsl_vector_ptr(this->theVector, index);
1079
1080
} // double *VB_Vector::getElementPtr(size_t index) const
1081
1082
/*********************************************************************
1083
* This method simply returns a constant pointer to the specified *
1084
* element. NOTE: This method does not allow the value of the *
1085
* specified element to be changed since it returns a constant *
1086
* pointer. *
1087
*********************************************************************/
1088
const double *VB_Vector::getElementConstPtr(const size_t index) const
1089
{
1090
1091
/*********************************************************************
1092
* We now check to make sure that index is within *
1093
* [0, this->theVector->size - 1]. The try/catch blocks are used to *
1094
* process any exception, if one is thrown by *
1095
* this->checkVectorRange(). *
1096
*********************************************************************/
1097
try
1098
{
1099
this->checkVectorRange(index, __LINE__, __FILE__, __FUNCTION__);
1100
} // try
1101
catch (GenericExcep &e)
1102
{
1103
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1104
} // catch
1105
1106
/*********************************************************************
1107
* Calling gsl_vector_ptr() to return the pointer to the specified *
1108
* vector element. *
1109
*********************************************************************/
1110
return gsl_vector_ptr(this->theVector, index);
1111
1112
} // const double *VB_Vector::getElementConstPtr(const size_t index)
1113
1114
/*********************************************************************
1115
* This method sets the specified element to the specified value. *
1116
*********************************************************************/
1117
void VB_Vector::setElement(const size_t index, const double value)
1118
{
1119
1120
/*********************************************************************
1121
* We now check to make sure that index is within *
1122
* [0, this->theVector->size - 1]. The try/catch blocks are used to *
1123
* process any exception, if one is thrown by *
1124
* this->checkVectorRange(). *
1125
*********************************************************************/
1126
try
1127
{
1128
this->checkVectorRange(index, __LINE__, __FILE__, __FUNCTION__);
1129
} // try
1130
catch (GenericExcep &e)
1131
{
1132
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1133
} // catch
1134
1135
/*********************************************************************
1136
* Calling gsl_vector_set() to set the specified element to the *
1137
* specified value. *
1138
*********************************************************************/
1139
gsl_vector_set(this->theVector, index, value);
1140
1141
} // void VB_Vector::setElement(const size_t index, const double value)
1142
1143
/*********************************************************************
1144
* This method simply calculates the Euclidean inner product of this *
1145
* instance of VB_Vector and the input VB_Vector, V2. *
1146
*********************************************************************/
1147
double VB_Vector::euclideanProduct(const VB_Vector& V2) const
1148
{
1149
1150
/*********************************************************************
1151
* Now returning the Euclidean inner product. *
1152
*********************************************************************/
1153
return this->euclideanProduct(V2.theVector);
1154
1155
} // double VB_Vector::euclideanProduct(const VB_Vector& V2) const
1156
1157
/*********************************************************************
1158
* This method simply calculates the Euclidean inner product of this *
1159
* instance of VB_Vector and the input VB_Vector, V2. *
1160
*********************************************************************/
1161
double VB_Vector::euclideanProduct(const VB_Vector *V2) const
1162
{
1163
1164
/*********************************************************************
1165
* Now returning the Euclidean inner product. *
1166
*********************************************************************/
1167
return this->euclideanProduct(V2->theVector);
1168
1169
} // double VB_Vector::euclideanProduct(const VB_Vector *V2) const
1170
1171
/*********************************************************************
1172
* This method simply calculates the Euclidean inner product of this *
1173
* instance of VB_Vector and the input gsl_vector struct, V2. *
1174
*********************************************************************/
1175
double VB_Vector::euclideanProduct(const gsl_vector& V2) const
1176
{
1177
1178
/*********************************************************************
1179
* Now returning the Euclidean inner product. *
1180
*********************************************************************/
1181
return this->euclideanProduct(&V2);
1182
1183
} // double VB_Vector::euclideanProduct(const gsl_vector& V2) const
1184
1185
/*********************************************************************
1186
* This method simply calculates the Euclidean inner product of this *
1187
* instance of VB_Vector and the input gsl_vector struct, V2. *
1188
*********************************************************************/
1189
double VB_Vector::euclideanProduct(const gsl_vector *V2) const
1190
{
1191
1192
/*********************************************************************
1193
* By default, if gsl_blas_ddot() has an error, the GSL error *
1194
* handler will be invoked, which calls abort(), creating a core dump.*
1195
* To forgo this behavior, the GSL error handler is turned off and *
1196
* the current error handler is stored in *
1197
* VB_Vector::currentGSLErrorHandler. This is accomplished by the *
1198
* call to this->turnOffGSLErrorHandler(). *
1199
*********************************************************************/
1200
this->turnOffGSLErrorHandler();
1201
1202
/*********************************************************************
1203
* The following try/catch blocks are used to ensure that the vector *
1204
* lengths are equal. *
1205
*********************************************************************/
1206
try
1207
{
1208
VB_Vector::checkVectorLengths(this->theVector, V2, __LINE__, __FILE__, __FUNCTION__);
1209
} // try
1210
catch (GenericExcep &e)
1211
{
1212
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1213
} // catch
1214
1215
/*********************************************************************
1216
* result will hold the Euclidean inner product. *
1217
*********************************************************************/
1218
double result = 0.0;
1219
1220
/*********************************************************************
1221
* Calling gsl_blas_ddot() to compute the Eucliden inner product. *
1222
* It's return value will be processed by VB_Vector::checkGSLStatus() *
1223
* and if the return value is non-zero, then *
1224
* VB_Vector::checkGSLStatus() will throw (and catch) an exception *
1225
* with an appropriate error message (as returned by gsl_strerror()). *
1226
*********************************************************************/
1227
VB_Vector::checkGSLStatus(gsl_blas_ddot(this->theVector, V2, &result), __LINE__, __FILE__, __FUNCTION__);
1228
1229
/*********************************************************************
1230
* We now restore the standard GSL error handler by calling *
1231
* this->restoreGSLErrorHandler(). *
1232
*********************************************************************/
1233
this->restoreGSLErrorHandler();
1234
1235
/*********************************************************************
1236
* Returning the inner product. *
1237
*********************************************************************/
1238
return result;
1239
1240
} // double VB_Vector::euclideanProduct(const gsl_vector *V2) const
1241
1242
/*********************************************************************
1243
* This static method calculates the Euclidean inner product of the *
1244
* 2 input gsl_vector structs. *
1245
*********************************************************************/
1246
double VB_Vector::euclideanProduct(const gsl_vector *V1, const gsl_vector *V2)
1247
{
1248
1249
/*********************************************************************
1250
* Checking to see if the 2 input gsl_vector structs have the same *
1251
* length or not. *
1252
*********************************************************************/
1253
try
1254
{
1255
VB_Vector::checkVectorLengths(V1, V2, __LINE__, __FILE__, __FUNCTION__);
1256
} // try
1257
catch (GenericExcep &e)
1258
{
1259
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1260
} // catch
1261
1262
/*********************************************************************
1263
* result will hold the Euclidean inner product. *
1264
*********************************************************************/
1265
double result = 0.0;
1266
1267
/*********************************************************************
1268
* Calling gsl_blas_ddot() to compute the Eucliden inner product. It's*
1269
* return value will be processed by VB_Vector::checkGSLStatus(). If *
1270
* the return value is non-zero, then VB_Vector::checkGSLStatus() will*
1271
* throw (and catch) an exception with an appropriate error message *
1272
* (as returned by gsl_strerror()). *
1273
*********************************************************************/
1274
VB_Vector::checkGSLStatus(gsl_blas_ddot(V1, V2, &result), __LINE__, __FILE__, __FUNCTION__);
1275
1276
/*********************************************************************
1277
* Returning the Euclidean inner product. *
1278
*********************************************************************/
1279
return result;
1280
1281
} // double VB_Vector::euclideanProduct(const gsl_vector *V1, const gsl_vector *V2)
1282
1283
/*********************************************************************
1284
* This static method calculates the Euclidean inner product of the *
1285
* 2 input gsl_vector structs. *
1286
*********************************************************************/
1287
double VB_Vector::euclideanProduct(const gsl_vector& V1, const gsl_vector& V2)
1288
{
1289
1290
/*********************************************************************
1291
* Calling VB_Vector::euclideanProduct(const gsl_vector *, *
1292
* const gsl_vector *) to return the Euclidean inner product. *
1293
*********************************************************************/
1294
return VB_Vector::euclideanProduct(&V1, &V2);
1295
1296
} // double VB_Vector::euclideanProduct(const gsl_vector& V1, const gsl_vector& V2)
1297
1298
/*********************************************************************
1299
* This method computes the distance between this instance of *
1300
* VB_Vector and the input VB_Vector. *
1301
*********************************************************************/
1302
double VB_Vector::euclideanDistance(const VB_Vector& V2) const
1303
{
1304
1305
/*********************************************************************
1306
* Calling this->euclideanDistance(const gsl_vector *) const. *
1307
*********************************************************************/
1308
return this->euclideanDistance(V2.theVector);
1309
1310
} // double VB_Vector::euclideanDistance(const VB_Vector& V2) const
1311
1312
/*********************************************************************
1313
* This method computes the distance between this instance of *
1314
* VB_Vector and the input VB_Vector. *
1315
*********************************************************************/
1316
double VB_Vector::euclideanDistance(const VB_Vector *V2) const
1317
{
1318
1319
/*********************************************************************
1320
* Calling this->euclideanDistance(const gsl_vector *) const. *
1321
*********************************************************************/
1322
return this->euclideanDistance(V2->theVector);
1323
1324
} // double VB_Vector::euclideanDistance(const VB_Vector *V2) const
1325
1326
/*********************************************************************
1327
* This method computes the distance between this instance of *
1328
* VB_Vector and the input gsl_vector struct. *
1329
*********************************************************************/
1330
double VB_Vector::euclideanDistance(const gsl_vector& V2) const
1331
{
1332
1333
/*********************************************************************
1334
* Calling this->euclideanDistance(const gsl_vector *) const. *
1335
*********************************************************************/
1336
return this->euclideanDistance(&V2);
1337
1338
} // double VB_Vector::euclideanDistance(const gsl_vector& V2) const
1339
1340
/*********************************************************************
1341
* This method computes the distance between this instance of *
1342
* VB_Vector and the input gsl_vector struct. *
1343
*********************************************************************/
1344
double VB_Vector::euclideanDistance(const gsl_vector *V2) const
1345
{
1346
1347
/*********************************************************************
1348
* Checking to see if the 2 vectors have the same lengths or not. *
1349
*********************************************************************/
1350
try
1351
{
1352
VB_Vector::checkVectorLengths(this->theVector, V2, __LINE__,
1353
__FILE__, __FUNCTION__);
1354
} // try
1355
catch (GenericExcep &e)
1356
{
1357
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1358
} // catch
1359
1360
/*********************************************************************
1361
* By default, if the memory allocation fails, *
1362
* then gsl_vector_calloc() will invoke the GSL error handler *
1363
* which calls abort(), creating a core dump. To forgo this behavior, *
1364
* the GSL error handler is turned off and the current error handler *
1365
* is stored in VB_Vector::currentGSLErrorHandler. This is *
1366
* accomplished by the call to this->turnOffGSLErrorHandler(). *
1367
*********************************************************************/
1368
this->turnOffGSLErrorHandler();
1369
1370
/*********************************************************************
1371
* tempVec will hold a copy of this instance of VB_Vector. tempVec is *
1372
* allocated sufficient memory and then this->theVector is copied to *
1373
* it. *
1374
*********************************************************************/
1375
gsl_vector *tempVec = gsl_vector_calloc(this->theVector->size);
1376
1377
/*********************************************************************
1378
* If tempVec is null, then we were unsuccessful in allocating memory *
1379
* for it. *
1380
*********************************************************************/
1381
try
1382
{
1383
VB_Vector::vectorNull(tempVec);
1384
} // try
1385
catch (GenericExcep &e)
1386
{
1387
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1388
} // catch
1389
1390
/*********************************************************************
1391
* Now copying the values from this->theVector to tempVec. The *
1392
* try/catch blocks are used to process the exception, if one is *
1393
* thrown. *
1394
*********************************************************************/
1395
try
1396
{
1397
this->GSLVectorMemcpy(tempVec, this->theVector);
1398
} // try
1399
catch (GenericExcep &e)
1400
{
1401
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1402
} // catch
1403
1404
/*********************************************************************
1405
* gsl_vector_sub() is called to compute the vector difference and its*
1406
* return value is passed to VB_Vector::checkGSLStatus(). If the *
1407
* return value is non-zero, then an error occurred with the GSL *
1408
* function. VB_Vector::checkGSLStatus() will then process the error *
1409
* by writing out an appropriate error message and then causing this *
1410
* program to exit. *
1411
*********************************************************************/
1412
VB_Vector::checkGSLStatus(gsl_vector_sub(tempVec, V2),
1413
__LINE__, __FILE__, __FUNCTION__);
1414
1415
/*********************************************************************
1416
* Now the Euclidean norm of tempVec is stored in distance, which is *
1417
* equal to the distance between this instance of VB_Vector and the *
1418
* input VB_Vector. *
1419
*********************************************************************/
1420
double distance = gsl_blas_dnrm2(tempVec);
1421
1422
/*********************************************************************
1423
* Now freeing the memory previously allocated to tempVec> *
1424
*********************************************************************/
1425
gsl_vector_free(tempVec);
1426
1427
/*********************************************************************
1428
* We now restore the standard GSL error handler by calling *
1429
* this->restoreGSLErrorHandler(). *
1430
*********************************************************************/
1431
this->restoreGSLErrorHandler();
1432
1433
/*********************************************************************
1434
* Returning the Euclidean distance. *
1435
*********************************************************************/
1436
return distance;
1437
1438
} // double VB_Vector::euclideanDistance(const gsl_vector *V2) const
1439
1440
/*********************************************************************
1441
* This overloaded operator returns the vector sum of this instance *
1442
* of VB_Vector and the input gsl_vector. *
1443
*********************************************************************/
1444
VB_Vector VB_Vector::operator+(const gsl_vector *V2) const
1445
{
1446
1447
/*********************************************************************
1448
* If the vector lengths are not equal, then a GenericExcep is thrown.*
1449
* The try/catch blocks are used to process the exception, which will *
1450
* have an appropriate error message. *
1451
*********************************************************************/
1452
try
1453
{
1454
VB_Vector::checkVectorLengths(this->theVector, V2, __LINE__, __FILE__, __FUNCTION__);
1455
} // try
1456
catch (GenericExcep &e)
1457
{
1458
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1459
} // catch
1460
1461
/*********************************************************************
1462
* vectorSum is instantiated to hold the vector sum. *
1463
*********************************************************************/
1464
VB_Vector vectorSum = VB_Vector(this->theVector->size);
1465
1466
/*********************************************************************
1467
* Now copying the values from this->theVector to vectorSum.theVector.*
1468
* The try/catch blocks are used to process the exception, if one is *
1469
* thrown. *
1470
*********************************************************************/
1471
try
1472
{
1473
this->GSLVectorMemcpy(vectorSum.theVector, this->theVector);
1474
} // try
1475
catch (GenericExcep &e)
1476
{
1477
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1478
} // catch
1479
1480
/*********************************************************************
1481
* gsl_vector_add() is called to compute the vector sum and its *
1482
* return value is passed to VB_Vector::checkGSLStatus(). If the *
1483
* return value is non-zero, then an error occurred with the GSL *
1484
* function. VB_Vector::checkGSLStatus() will then process the error *
1485
* by writing out an appropriate error message and then causing this *
1486
* program to exit. *
1487
*********************************************************************/
1488
VB_Vector::checkGSLStatus(gsl_vector_add(vectorSum.theVector, V2),
1489
__LINE__, __FILE__, __FUNCTION__);
1490
1491
/*********************************************************************
1492
* Now returning the vector sum. *
1493
*********************************************************************/
1494
return vectorSum;
1495
1496
} // VB_Vector VB_Vector::operator+(const gsl_vector *V2) const
1497
1498
/*********************************************************************
1499
* This overloaded operator returns the vector sum of this instance *
1500
* of VB_Vector and the input VB_Vector. *
1501
*********************************************************************/
1502
VB_Vector VB_Vector::operator+(const VB_Vector& V2) const
1503
{
1504
1505
/*********************************************************************
1506
* Returning ((*this) + V2.theVector). *
1507
*********************************************************************/
1508
return ((*this) + V2.theVector);
1509
1510
} // VB_Vector VB_Vector::operator+(const VB_Vector& V2) const
1511
1512
/*********************************************************************
1513
* This overloaded operator returns the vector difference of this *
1514
* instance of VB_Vector and the input gsl_vector. *
1515
*********************************************************************/
1516
VB_Vector VB_Vector::operator-(const gsl_vector *V2) const
1517
{
1518
1519
/*********************************************************************
1520
* If the vector lengths are not equal, then a GenericExcep is thrown.*
1521
* The try/catch blocks are used to process the exception, which will *
1522
* have an appropriate error message. *
1523
*********************************************************************/
1524
try
1525
{
1526
VB_Vector::checkVectorLengths(this->theVector, V2, __LINE__, __FILE__, __FUNCTION__);
1527
} // try
1528
catch (GenericExcep &e)
1529
{
1530
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1531
} // catch
1532
1533
/*********************************************************************
1534
* vectorDiff is instantiated to hold the vector difference. *
1535
*********************************************************************/
1536
VB_Vector vectorDiff = VB_Vector(this->theVector->size);
1537
1538
/*********************************************************************
1539
* Now copying the values from this->theVector to *
1540
* vectorDiff.theVector. The try/catch blocks are used to process the *
1541
* exception, if one is thrown. *
1542
*********************************************************************/
1543
try
1544
{
1545
this->GSLVectorMemcpy(vectorDiff.theVector, this->theVector);
1546
} // try
1547
catch (GenericExcep &e)
1548
{
1549
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1550
} // catch
1551
1552
/*********************************************************************
1553
* gsl_vector_sub() is called to compute the vector difference and its*
1554
* return value is passed to VB_Vector::checkGSLStatus(). If the *
1555
* return value is non-zero, then an error occurred with the GSL *
1556
* function. VB_Vector::checkGSLStatus() will then process the error *
1557
* by writing out an appropriate error message and then causing this *
1558
* program to exit. *
1559
*********************************************************************/
1560
VB_Vector::checkGSLStatus(gsl_vector_sub(vectorDiff.theVector, V2),
1561
__LINE__, __FILE__, __FUNCTION__);
1562
1563
/*********************************************************************
1564
* Now returning the vector difference. *
1565
*********************************************************************/
1566
return vectorDiff;
1567
1568
} // VB_Vector VB_Vector::operator-(const gsl_vector *V2) const
1569
1570
/*********************************************************************
1571
* This overloaded operator returns the vector difference of this *
1572
* instance of VB_Vector and the input VB_Vector. *
1573
*********************************************************************/
1574
VB_Vector VB_Vector::operator-(const VB_Vector& V2) const
1575
{
1576
1577
/*********************************************************************
1578
* Returning ((*this) - V2.theVector). *
1579
*********************************************************************/
1580
return ((*this) - V2.theVector);
1581
1582
} // VB_Vector VB_Vector::operator-(const VB_Vector& V2) const
1583
1584
/*********************************************************************
1585
* This overloaded operator returns the vector difference of this *
1586
* instance of VB_Vector and the input VB_Vector. *
1587
*********************************************************************/
1588
VB_Vector VB_Vector::operator-(const VB_Vector *V2) const
1589
{
1590
1591
/*********************************************************************
1592
* Returning the difference of (*this) and V2->theVector. *
1593
*********************************************************************/
1594
return ((*this) - V2->theVector);
1595
1596
} // VB_Vector VB_Vector::operator-(const VB_Vector *V2) const
1597
1598
/*********************************************************************
1599
* This overloaded operator returns the vector difference of this *
1600
* instance of VB_Vector and the input gsl_vector struct. *
1601
*********************************************************************/
1602
VB_Vector VB_Vector::operator-(const gsl_vector& V2) const
1603
{
1604
1605
/*********************************************************************
1606
* Returning the difference of (*this) and (&V2). *
1607
*********************************************************************/
1608
return ((*this) - (&V2));
1609
1610
} // VB_Vector VB_Vector::operator-(const gsl_vector& V2) const
1611
1612
/*********************************************************************
1613
* This is the overloaded [] operator. NOTE: A reference is returned, *
1614
* allowing assignment. *
1615
*********************************************************************/
1616
double& VB_Vector::operator[](const size_t index) const
1617
{
1618
1619
/*********************************************************************
1620
* Returning the specified reference. *
1621
*********************************************************************/
1622
return this->theVector->data[index];
1623
1624
} // double& VB_Vector::operator[](const size_t index) const
1625
1626
/*********************************************************************
1627
* Overloaded "()" operator. Behaves just like the overloaded "[]" *
1628
* operator for this class except index is range checked. *
1629
*********************************************************************/
1630
double& VB_Vector::operator()(const size_t index) const
1631
{
1632
1633
/*********************************************************************
1634
* A try/catch block is used to determine if index is out of range or *
1635
* not. *
1636
*********************************************************************/
1637
try
1638
{
1639
1640
/*********************************************************************
1641
* We now check to make sure that index is within *
1642
* [0, this->theVector->size - 1]. *
1643
*********************************************************************/
1644
this->checkVectorRange(index, __LINE__, __FILE__, __FUNCTION__);
1645
1646
} // try
1647
catch (GenericExcep &e)
1648
{
1649
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1650
} // catch
1651
1652
/*********************************************************************
1653
* Returning the specified reference. *
1654
*********************************************************************/
1655
return this->theVector->data[index];
1656
1657
} // double& VB_Vector::operator()(const size_t index) const
1658
1659
/*********************************************************************
1660
* This overloaded operator checks for equality between this instance *
1661
* of VB_Vector and the input VB_Vector. *
1662
*********************************************************************/
1663
bool VB_Vector::operator==(const VB_Vector *V2) const
1664
{
1665
1666
/*********************************************************************
1667
* Returning the value of the call to operator==(gsl_vector *). *
1668
*********************************************************************/
1669
return ((*this) == V2->theVector);
1670
1671
} // bool VB_Vector::operator==(const VB_Vector *V2) const
1672
1673
/*********************************************************************
1674
* This overloaded operator checks for equality between this instance *
1675
* of VB_Vector and the input VB_Vector. *
1676
*********************************************************************/
1677
bool VB_Vector::operator==(const VB_Vector& V2) const
1678
{
1679
1680
/*********************************************************************
1681
* Returning the value of the call to operator==(gsl_vector *). *
1682
*********************************************************************/
1683
return ((*this) == V2.theVector);
1684
1685
} // bool VB_Vector::operator==(const VB_Vector& V2) const
1686
1687
/*********************************************************************
1688
* This overloaded operator checks for equality between this instance *
1689
* of VB_Vector and the input gsl_vector struct. *
1690
*********************************************************************/
1691
bool VB_Vector::operator==(const gsl_vector *V2) const
1692
{
1693
1694
/*********************************************************************
1695
* If both this->theVector and V2 are NULL, then true is returned. *
1696
*********************************************************************/
1697
if (this->theVector==NULL && V2==NULL)
1698
return true;
1699
1700
/*********************************************************************
1701
* If one of {this->theVector, V2} is NULL and the other is not, then *
1702
* false is returned. *
1703
*********************************************************************/
1704
if ( (this->theVector && V2==NULL) || (this->theVector==NULL && V2))
1705
return false;
1706
1707
/*********************************************************************
1708
* If this->theVector->size and V2->size are not equal, then false is *
1709
* returned. *
1710
*********************************************************************/
1711
if (this->theVector->size != V2->size)
1712
{
1713
return false;
1714
} // if
1715
1716
/*********************************************************************
1717
* memcmp() returns 0 if the 2 memory areas have the same content. *
1718
*********************************************************************/
1719
return (memcmp(this->theVector, V2, this->theVector->size * sizeof(double)) == 0);
1720
1721
} // bool VB_Vector::operator==(const gsl_vector *V2) const
1722
1723
/*********************************************************************
1724
* This overloaded operator checks for equality between this instance *
1725
* of VB_Vector and the input gsl_vector struct. *
1726
*********************************************************************/
1727
bool VB_Vector::operator==(const gsl_vector& V2) const
1728
{
1729
1730
/*********************************************************************
1731
* Returning the value of the call to operator==(gsl_vector *). *
1732
*********************************************************************/
1733
return ((*this) == (&V2));
1734
1735
} // bool VB_Vector::operator==(const gsl_vector& V2) const
1736
1737
/*********************************************************************
1738
* This overloaded operator checks for inequality between this *
1739
* instance of VB_Vector and the input VB_Vector. *
1740
*********************************************************************/
1741
bool VB_Vector::operator!=(const VB_Vector *V2) const
1742
{
1743
1744
/*********************************************************************
1745
* Returning the negation of the equality check. *
1746
*********************************************************************/
1747
return !((*this) == V2);
1748
1749
} // bool VB_Vector::operator!=(const VB_Vector *V2) const
1750
1751
/*********************************************************************
1752
* This overloaded operator checks for inequality between this *
1753
* instance of VB_Vector and the input VB_Vector. *
1754
*********************************************************************/
1755
bool VB_Vector::operator!=(const VB_Vector& V2) const
1756
{
1757
1758
/*********************************************************************
1759
* Returning the negation of the equality check. *
1760
*********************************************************************/
1761
return !((*this) == V2);
1762
1763
} // bool VB_Vector::operator!=(const VB_Vector& V2) const
1764
1765
/*********************************************************************
1766
* This overloaded operator checks for inequality between this *
1767
* instance of VB_Vector and the input gsl_vector struct. *
1768
*********************************************************************/
1769
bool VB_Vector::operator!=(const gsl_vector& V2) const
1770
{
1771
1772
/*********************************************************************
1773
* Returning the negation of the equality check. *
1774
*********************************************************************/
1775
return !((*this) == V2);
1776
1777
} // bool VB_Vector::operator!=(const gsl_vector& V2) const
1778
1779
/*********************************************************************
1780
* This overloaded operator checks for inequality between this *
1781
* instance of VB_Vector and the input gsl_vector struct. *
1782
*********************************************************************/
1783
bool VB_Vector::operator!=(const gsl_vector *V2) const
1784
{
1785
1786
/*********************************************************************
1787
* Returning the negation of the equality check. *
1788
*********************************************************************/
1789
return !((*this) == V2);
1790
1791
} // bool VB_Vector::operator!=(const gsl_vector *V2) const
1792
1793
const VB_Vector&
1794
VB_Vector::operator=(const VB_Vector& V2)
1795
{
1796
if (this==(&V2))
1797
return *this;
1798
if (V2.getLength()==0) {
1799
clear();
1800
return *this;
1801
}
1802
this->init(this->valid, V2.dataType,V2.fileFormat);
1803
this->init(V2.getLength());
1804
if (!(this->theVector))
1805
return *this;
1806
this->fileName=V2.fileName;
1807
this->GSLVectorMemcpy(this->theVector, V2.theVector);
1808
return *this;
1809
}
1810
1811
const VB_Vector &
1812
VB_Vector::zero()
1813
{
1814
gsl_vector_set_zero(this->theVector);
1815
return *this;
1816
}
1817
1818
/*********************************************************************
1819
* This static method returns a VB_Vector object that represents a *
1820
* basis vector. *
1821
*********************************************************************/
1822
VB_Vector VB_Vector::getBasisVector(const size_t length, const size_t index)
1823
{
1824
1825
/*********************************************************************
1826
* Instantiating a VB_Vector whose length will be length. NOTE: All *
1827
* the elements will be zero. *
1828
*********************************************************************/
1829
VB_Vector basis(length);
1830
1831
/*********************************************************************
1832
* Now setting the appropriate element to 1. *
1833
*********************************************************************/
1834
basis.setElement(index, 1.0);
1835
1836
/*********************************************************************
1837
* Now returning the basis VB_Vector. *
1838
*********************************************************************/
1839
return basis;
1840
1841
} // VB_Vector VB_Vector::getBasisVector(const size_t length, const size_t index)
1842
1843
/*********************************************************************
1844
* This method scales the elements of this instance of VB_Vector by *
1845
* alpha and returns the scaled vector. *
1846
*********************************************************************/
1847
VB_Vector VB_Vector::vectorScale(const double alpha) const
1848
{
1849
1850
/*********************************************************************
1851
* By default, if the memory allocation fails, *
1852
* then gsl_vector_calloc() will invoke the GSL error handler *
1853
* which calls abort(), creating a core dump. To forgo this behavior, *
1854
* the GSL error handler is turned off and the current error handler *
1855
* is stored in VB_Vector::currentGSLErrorHandler. This is *
1856
* accomplished by the call to this->turnOffGSLErrorHandler(). *
1857
*********************************************************************/
1858
this->turnOffGSLErrorHandler();
1859
1860
/*********************************************************************
1861
* tempVec will hold a copy of this instance of VB_Vector. tempVec is *
1862
* allocated sufficient memory and then this->theVector is copied to *
1863
* it. *
1864
*********************************************************************/
1865
gsl_vector *tempVec = gsl_vector_calloc(this->getLength());
1866
1867
/*********************************************************************
1868
* We now restore the standard GSL error handler by calling *
1869
* this->restoreGSLErrorHandler(). *
1870
*********************************************************************/
1871
this->restoreGSLErrorHandler();
1872
1873
/*********************************************************************
1874
* If tempVec is null, then we were unsuccessful in allocating memory *
1875
* for it. *
1876
*********************************************************************/
1877
try
1878
{
1879
VB_Vector::vectorNull(tempVec);
1880
} // try
1881
catch (GenericExcep &e)
1882
{
1883
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1884
} // catch
1885
1886
/*********************************************************************
1887
* Now copying the values from this->theVector to tempVec. The *
1888
* try/catch blocks are used to process the exception, if one is *
1889
* thrown. *
1890
*********************************************************************/
1891
try
1892
{
1893
this->GSLVectorMemcpy(tempVec, this->theVector);
1894
} // try
1895
catch (GenericExcep &e)
1896
{
1897
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
1898
} // catch
1899
1900
/*********************************************************************
1901
* gsl_vector_scale() is called to do the scaling and its *
1902
* return value is passed to VB_Vector::checkGSLStatus(). If the *
1903
* return value is non-zero, then an error occurred with the GSL *
1904
* function. VB_Vector::checkGSLStatus() will then process the error *
1905
* by writing out an appropriate error message and then causing this *
1906
* program to exit. *
1907
*********************************************************************/
1908
VB_Vector::checkGSLStatus(gsl_vector_scale(tempVec, alpha),
1909
__LINE__, __FILE__, __FUNCTION__);
1910
1911
/*********************************************************************
1912
* Now returning the scaled vector, tempVec, as a VB_Vector object. *
1913
*********************************************************************/
1914
return VB_Vector(tempVec);
1915
1916
} // VB_Vector VB_Vector::vectorScale(const double alpha) const
1917
1918
/*********************************************************************
1919
* This method adds the constant alpha to all elements of the *
1920
* vector. *
1921
*********************************************************************/
1922
void VB_Vector::addConstant(const double alpha)
1923
{
1924
1925
/*********************************************************************
1926
* gsl_vector_add_constant() is called to add the constant and its *
1927
* return value is passed to VB_Vector::checkGSLStatus(). If the *
1928
* return value is non-zero, then an error occurred with the GSL *
1929
* function. VB_Vector::checkGSLStatus() will then process the error *
1930
* by writing out an appropriate error message and then causing this *
1931
* program to exit. *
1932
*********************************************************************/
1933
VB_Vector::checkGSLStatus(gsl_vector_add_constant(this->theVector, alpha),
1934
__LINE__, __FILE__, __FUNCTION__);
1935
1936
} // void VB_Vector::addConstant(const double alpha)
1937
1938
/*********************************************************************
1939
* This overloaded operator returns the vector sum of this instance *
1940
* of VB_Vector and the input VB_Vector. *
1941
*********************************************************************/
1942
VB_Vector VB_Vector::operator+(const VB_Vector *V2) const
1943
{
1944
1945
/*********************************************************************
1946
* Returning the sum of (*this) and V2->theVector. *
1947
*********************************************************************/
1948
return ((*this) + V2->theVector);
1949
1950
} // VB_Vector VB_Vector::operator+(const VB_Vector *V2) const
1951
1952
/*********************************************************************
1953
* This overloaded operator returns the vector sum of this instance *
1954
* of VB_Vector and the input gsl_vector struct. *
1955
*********************************************************************/
1956
VB_Vector VB_Vector::operator+(const gsl_vector& V2) const
1957
{
1958
1959
/*********************************************************************
1960
* Returning the sum of (*this) and (&V2). *
1961
*********************************************************************/
1962
return ((*this) + (&V2));
1963
1964
} // VB_Vector VB_Vector::operator+(const gsl_vector& V2) const
1965
1966
/*********************************************************************
1967
* This overloaded "+" operator is a friend function. This overloaded *
1968
* operator is implemented to preserve the commutative property for *
1969
* addition between VB_Vector and gsl_vector. *
1970
*********************************************************************/
1971
VB_Vector operator+(const gsl_vector *V1, const VB_Vector& V2)
1972
{
1973
1974
/*********************************************************************
1975
* Returning the sum of V2 and V1. *
1976
*********************************************************************/
1977
return (V2 + V1);
1978
1979
} // VB_Vector operator+(const gsl_vector *V1, const VB_Vector& V2)
1980
1981
/*********************************************************************
1982
* This overloaded "+" operator is a friend function. This overloaded *
1983
* operator is implemented to preserve the commutative property for *
1984
* addition between VB_Vector and gsl_vector. *
1985
*********************************************************************/
1986
VB_Vector operator+(const gsl_vector& V1, const VB_Vector *V2)
1987
{
1988
1989
/*********************************************************************
1990
* Returning the sum of V2 and V1. *
1991
*********************************************************************/
1992
return ((*V2) + (&V1));
1993
1994
} // VB_Vector operator+(const gsl_vector& V1, const VB_Vector *V2)
1995
1996
/*********************************************************************
1997
* This overloaded "+" operator is a friend function. This overloaded *
1998
* operator is implemented to preserve the commutative property for *
1999
* addition between VB_Vector and gsl_vector. *
2000
*********************************************************************/
2001
VB_Vector operator+(const gsl_vector& V1, const VB_Vector& V2)
2002
{
2003
2004
/*********************************************************************
2005
* Returning the sum of V2 and V1. *
2006
*********************************************************************/
2007
return (V2 + (&V1));
2008
2009
} // VB_Vector operator+(const gsl_vector& V1, const VB_Vector& V2)
2010
2011
/*********************************************************************
2012
* This overloaded operator is implemented as a friend function so *
2013
* that commutativity can be preserved for "==". *
2014
*********************************************************************/
2015
bool operator==(const gsl_vector *V1, const VB_Vector& V2)
2016
{
2017
2018
/*********************************************************************
2019
* If V2.theVector->size and V1->size are not equal, then false is *
2020
* returned. *
2021
*********************************************************************/
2022
if (V2.getLength() != V1->size)
2023
{
2024
return false;
2025
} // if
2026
2027
/*********************************************************************
2028
* memcmp() returns 0 if the 2 memory areas are equal. *
2029
*********************************************************************/
2030
return (memcmp(V2.theVector, V1, V2.getLength() * sizeof(double)) == 0);
2031
2032
} // bool operator==(const gsl_vector *V1, const VB_Vector& V2)
2033
2034
/*********************************************************************
2035
* This overloaded operator is implemented as a friend function so *
2036
* that commutativity can be preserved for "==". *
2037
*********************************************************************/
2038
bool operator==(const gsl_vector& V1, const VB_Vector& V2)
2039
{
2040
2041
/*********************************************************************
2042
* memcmp() returns 0 if the 2 memory areas are equal. *
2043
*********************************************************************/
2044
return (memcmp(V2.theVector, &V1, V2.getLength() * sizeof(double)) == 0);
2045
2046
} // bool operator==(const gsl_vector& V1, const VB_Vector& V2)
2047
2048
/*********************************************************************
2049
* This overloaded "!=" is implemented as a friend function to *
2050
* preserve commutativity. *
2051
*********************************************************************/
2052
bool operator!=(const gsl_vector *V1, const VB_Vector& V2)
2053
{
2054
2055
/*********************************************************************
2056
* Returning the negation of the equality check. *
2057
*********************************************************************/
2058
return !(V2 == V1);
2059
2060
} // bool operator!=(const gsl_vector *V1, const VB_Vector& V2)
2061
2062
/*********************************************************************
2063
* This overloaded "!=" is implemented as a friend function to *
2064
* preserve commutativity. *
2065
*********************************************************************/
2066
bool operator!=(const gsl_vector& V1, const VB_Vector& V2)
2067
{
2068
2069
/*********************************************************************
2070
* Returning the negation of the equality check. *
2071
*********************************************************************/
2072
return !(V2 == V1);
2073
2074
} // bool operator!=(const gsl_vector& V1, const VB_Vector& V2)
2075
2076
/*********************************************************************
2077
* This overloaded "<<" operator is implemented as a friend function. *
2078
*********************************************************************/
2079
ostream& operator<<(ostream& outStream, const VB_Vector& V)
2080
{
2081
2082
/*********************************************************************
2083
* Now calling operator<<(ostream&, const VB_Vector *), whose return *
2084
* value, in turn, is returned by this function. *
2085
*********************************************************************/
2086
return (outStream << (&V));
2087
2088
} // ostream& operator<<(ostream& outStream, const VB_Vector& V)
2089
2090
/*********************************************************************
2091
* This overloaded "<<" operator is implemented as a friend function. *
2092
*********************************************************************/
2093
ostream& operator<<(ostream& outStream, const VB_Vector *V)
2094
{
2095
2096
/*********************************************************************
2097
* Printing out the file name, valid flag, data type, and file type *
2098
* data members. *
2099
*********************************************************************/
2100
outStream << "Vector File Name = [" << V->fileName << "]" << endl;
2101
outStream << "Vector Valid = [" << V->valid << "]" << endl;
2102
outStream << "Vector Data Type = [" << DataTypeName(V->dataType) << "]" << endl;
2103
// DYK replaced the following
2104
// outStream << "Vector File Type = [" << (V->fileFormat ? V->fileFormat->getName() : "NONE") << "]" << endl;
2105
outStream << "Vector File Type = [" << V->fileFormat.getName() << "]" << endl;
2106
2107
/*********************************************************************
2108
* If the input VB_Vector does not have a null theVector data member, *
2109
* then the vector elements are printed out. *
2110
*********************************************************************/
2111
if (V->theVector)
2112
{
2113
2114
/*********************************************************************
2115
* Writing out the fields from the data member theVector, a *
2116
* gsl_vector struct. *
2117
*********************************************************************/
2118
outStream << "gsl_vector stride = [" << V->theVector->stride << "]" << endl;
2119
outStream << "gsl_vector owner = [" << V->theVector->owner << "]" << endl;
2120
outStream << "Vector Length = [" << V->theVector->size << "]" << endl;
2121
2122
/*********************************************************************
2123
* The following for loop is used to write out the vector elements. *
2124
*********************************************************************/
2125
for (size_t i = 0; i < V->theVector->size; i++)
2126
{
2127
outStream << "element[" << i << "] = [" << V->theVector->data[i] << "]" << endl;
2128
} // for i
2129
2130
} // if
2131
2132
/*********************************************************************
2133
* If program flow ends up here, then V->getTheVector() returned a *
2134
* null gsl_vector. *
2135
*********************************************************************/
2136
else
2137
{
2138
outStream << "NULL gsl_vector." << endl;
2139
}
2140
2141
/*********************************************************************
2142
* Now returning outStream. *
2143
*********************************************************************/
2144
return outStream;
2145
2146
} // ostream& operator<<(ostream& outStream, const VB_Vector *V)
2147
2148
/*********************************************************************
2149
* This is a static method to compute the Euclidean distance between *
2150
* the two input gsl_vector structs. *
2151
*********************************************************************/
2152
double VB_Vector::euclideanDistance(const gsl_vector *V1, const gsl_vector *V2)
2153
{
2154
2155
/*********************************************************************
2156
* tempVec is instantiated from V1 as a VB_Vector object. *
2157
*********************************************************************/
2158
VB_Vector tempVec(V1);
2159
2160
/*********************************************************************
2161
* Now calling tempVec.euclideanDistance() to return the Euclidean *
2162
* distance between V1 and V2. *
2163
*********************************************************************/
2164
return tempVec.euclideanDistance(V2);
2165
2166
} // double VB_Vector::euclideanDistance(const gsl_vector *V1, const gsl_vector *V2)
2167
2168
/*********************************************************************
2169
* This is a static method to compute the Euclidean distance between *
2170
* the two input gsl_vector structs. *
2171
*********************************************************************/
2172
double VB_Vector::euclideanDistance(const gsl_vector& V1, const gsl_vector& V2)
2173
{
2174
2175
/*********************************************************************
2176
* tempVec is instantiated from V1 as a VB_Vector object. *
2177
*********************************************************************/
2178
VB_Vector tempVec(V1);
2179
2180
/*********************************************************************
2181
* Now calling tempVec.euclideanDistance() to return the Euclidean *
2182
* distance between V1 and V2. *
2183
*********************************************************************/
2184
return tempVec.euclideanDistance(V2);
2185
2186
} // double VB_Vector::euclideanDistance(const gsl_vector& V1, const gsl_vector& V2)
2187
2188
/*********************************************************************
2189
* This method adds V to this instance of VB_Vector in place (this *
2190
* instance of VB_Vector gets overwritten by the vector sum). *
2191
*********************************************************************/
2192
void VB_Vector::addInPlace(const VB_Vector& V)
2193
{
2194
2195
/*********************************************************************
2196
* Calling this->addInPlace(gsl_vector *) to carry out the addition. *
2197
*********************************************************************/
2198
this->addInPlace(V.theVector);
2199
2200
} // void VB_Vector::addInPlace(const VB_Vector& V)
2201
2202
/*********************************************************************
2203
* This method adds V to this instance of VB_Vector in place (this *
2204
* instance of VB_Vector gets overwritten by the vector sum). *
2205
*********************************************************************/
2206
void VB_Vector::addInPlace(const VB_Vector *V)
2207
{
2208
2209
/*********************************************************************
2210
* Calling this->addInPlace(gsl_vector *) to carry out the addition. *
2211
*********************************************************************/
2212
this->addInPlace(V->theVector);
2213
2214
} // void VB_Vector::addInPlace(const VB_Vector *V)
2215
2216
/*********************************************************************
2217
* This method adds V to this instance of VB_Vector in place (this *
2218
* instance of VB_Vector gets overwritten by the vector sum). *
2219
*********************************************************************/
2220
void VB_Vector::addInPlace(const gsl_vector& V)
2221
{
2222
2223
/*********************************************************************
2224
* Calling this->addInPlace(gsl_vector *) to carry out the addition. *
2225
*********************************************************************/
2226
this->addInPlace(&V);
2227
2228
} // void VB_Vector::addInPlace(const gsl_vector& V)
2229
2230
/*********************************************************************
2231
* This method adds V to this instance of VB_Vector in place (this *
2232
* instance of VB_Vector gets overwritten by the vector sum). *
2233
*********************************************************************/
2234
void VB_Vector::addInPlace(const gsl_vector *V)
2235
{
2236
2237
/*********************************************************************
2238
* gsl_vector_add() is called to compute the vector sum and its *
2239
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2240
* return value is non-zero, then a error occurred with the GSL *
2241
* function. VB_Vector::checkGSLStatus() will then process the error *
2242
* by writing out an appropriate error message and then causing this *
2243
* program to exit. *
2244
*********************************************************************/
2245
VB_Vector::checkGSLStatus(gsl_vector_add(this->theVector, V),
2246
__LINE__, __FILE__, __FUNCTION__);
2247
2248
} // void VB_Vector::addInPlace(const gsl_vector *V)
2249
2250
/*********************************************************************
2251
* This method subtracts V from this instance of VB_Vector in place *
2252
* (this instance of VB_Vector gets overwritten by the vector *
2253
* difference). *
2254
*********************************************************************/
2255
void VB_Vector::subtractInPlace(const VB_Vector& V)
2256
{
2257
2258
/*********************************************************************
2259
* Calling this->subtractInPlace(gsl_vector *) to carry out the *
2260
* subtraction. *
2261
*********************************************************************/
2262
this->subtractInPlace(V.theVector);
2263
2264
} // void VB_Vector::subtractInPlace(const VB_Vector& V)
2265
2266
/*********************************************************************
2267
* This method subtracts V from this instance of VB_Vector in place *
2268
* (this instance of VB_Vector gets overwritten by the vector *
2269
* difference). *
2270
*********************************************************************/
2271
void VB_Vector::subtractInPlace(const VB_Vector *V)
2272
{
2273
2274
/*********************************************************************
2275
* Calling this->subtractInPlace(gsl_vector *) to carry out the *
2276
* subtraction. *
2277
*********************************************************************/
2278
this->subtractInPlace(V->theVector);
2279
2280
} // void VB_Vector::subtractInPlace(const VB_Vector *V)
2281
2282
/*********************************************************************
2283
* This method subtracts V from this instance of VB_Vector in place *
2284
* (this instance of VB_Vector gets overwritten by the vector *
2285
* difference). *
2286
*********************************************************************/
2287
void VB_Vector::subtractInPlace(const gsl_vector& V)
2288
{
2289
2290
/*********************************************************************
2291
* Calling this->subtractInPlace(gsl_vector *) to carry out the *
2292
* subtraction. *
2293
*********************************************************************/
2294
this->subtractInPlace(&V);
2295
2296
} // void VB_Vector::subtractInPlace(const gsl_vector& V)
2297
2298
/*********************************************************************
2299
* This method subtracts V from this instance of VB_Vector in place *
2300
* (this instance of VB_Vector gets overwritten by the vector *
2301
* difference). *
2302
*********************************************************************/
2303
void VB_Vector::subtractInPlace(const gsl_vector *V)
2304
{
2305
2306
/*********************************************************************
2307
* gsl_vector_sub() is called to compute the difference and its *
2308
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2309
* return value is non-zero, then a error occurred with the GSL *
2310
* function. VB_Vector::checkGSLStatus() will then process the error *
2311
* by writing out an appropriate error message and then causing this *
2312
* program to exit. *
2313
*********************************************************************/
2314
VB_Vector::checkGSLStatus(gsl_vector_sub(this->theVector, V),
2315
__LINE__, __FILE__, __FUNCTION__);
2316
2317
} // void VB_Vector::subtractInPlace(const gsl_vector *V)
2318
2319
/*********************************************************************
2320
* This method scales the elements of this instance of VB_Vector by *
2321
* alpha in place. *
2322
*********************************************************************/
2323
void VB_Vector::scaleInPlace(const double alpha)
2324
{
2325
2326
/*********************************************************************
2327
* gsl_vector_scale() is called to do the scaling and its *
2328
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2329
* return value is non-zero, then a error occurred with the GSL *
2330
* function. VB_Vector::checkGSLStatus() will then process the error *
2331
* by writing out an appropriate error message and then causing this *
2332
* program to exit. *
2333
*********************************************************************/
2334
VB_Vector::checkGSLStatus(gsl_vector_scale(this->theVector, alpha),
2335
__LINE__, __FILE__, __FUNCTION__);
2336
2337
} // void VB_Vector::scaleInPlace(const double alpha)
2338
2339
/*********************************************************************
2340
* This overloaded "*" operator returns the Euclidean inner product of*
2341
* this instance of VB_Vector with the input vector. *
2342
*********************************************************************/
2343
double VB_Vector::operator*(const VB_Vector *V2) const
2344
{
2345
2346
/*********************************************************************
2347
* Now returning the Euclidean inner product. *
2348
*********************************************************************/
2349
return this->euclideanProduct(V2->theVector);
2350
2351
} // double VB_Vector::operator*(const VB_Vector *V2) const
2352
2353
/*********************************************************************
2354
* This overloaded "*" operator returns the Euclidean inner product of*
2355
* this instance of VB_Vector with the input vector. *
2356
*********************************************************************/
2357
double VB_Vector::operator*(const VB_Vector& V2) const
2358
{
2359
2360
/*********************************************************************
2361
* Now returning the Euclidean inner product. *
2362
*********************************************************************/
2363
return this->euclideanProduct(V2.theVector);
2364
2365
} // double VB_Vector::operator*(const VB_Vector& V2) const
2366
2367
/*********************************************************************
2368
* This overloaded "*" operator returns the Euclidean inner product of*
2369
* this instance of VB_Vector with the input vector. *
2370
*********************************************************************/
2371
double VB_Vector::operator*(const gsl_vector& V2) const
2372
{
2373
2374
/*********************************************************************
2375
* Now returning the Euclidean inner product. *
2376
*********************************************************************/
2377
return this->euclideanProduct(&V2);
2378
2379
} // double VB_Vector::operator*(const gsl_vector& V2) const
2380
2381
/*********************************************************************
2382
* This overloaded "*" operator returns the Euclidean inner product of*
2383
* this instance of VB_Vector with the input vector. *
2384
*********************************************************************/
2385
double VB_Vector::operator*(const gsl_vector *V2) const
2386
{
2387
2388
/*********************************************************************
2389
* Now returning the Euclidean inner product. *
2390
*********************************************************************/
2391
return this->euclideanProduct(V2);
2392
2393
} // double VB_Vector::operator*(const gsl_vector& V2) const
2394
2395
/*********************************************************************
2396
* Overloaded "*" operator (a friend function) that returns the *
2397
* Euclidean inner product of the 2 input vectors. *
2398
*********************************************************************/
2399
double operator*(const gsl_vector *V1, const VB_Vector& V2)
2400
{
2401
2402
/*********************************************************************
2403
* Now calling V2.euclideanProduct(V1) to compute the Euclidean inner *
2404
* product. *
2405
*********************************************************************/
2406
return V2.euclideanProduct(V1);
2407
2408
} // double operator*(const gsl_vector *V1, const VB_Vector& V2)
2409
2410
/*********************************************************************
2411
* Overloaded "*" operator (a friend function) that returns the *
2412
* Euclidean inner product of the 2 input vectors. *
2413
*********************************************************************/
2414
double operator*(const gsl_vector& V1, const VB_Vector& V2)
2415
{
2416
2417
/*********************************************************************
2418
* Now calling V2.euclideanProduct(V1) to compute the Euclidean inner *
2419
* product. *
2420
*********************************************************************/
2421
return V2.euclideanProduct(&V1);
2422
2423
} // double operator*(const gsl_vector& V1, const VB_Vector& V2)
2424
2425
/*********************************************************************
2426
* Overloaded "*" operator (a friend function) that returns the *
2427
* Euclidean inner product of the 2 input vectors. *
2428
*********************************************************************/
2429
double operator*(const gsl_vector& V1, const VB_Vector *V2)
2430
{
2431
2432
/*********************************************************************
2433
* Now calling V2.euclideanProduct(V1) to compute the Euclidean inner *
2434
* product. *
2435
*********************************************************************/
2436
return V2->euclideanProduct(&V1);
2437
2438
} // double operator*(const gsl_vector& V1, const VB_Vector *V2)
2439
2440
/*********************************************************************
2441
* This overloaded operator adds the input scalar to all elements of *
2442
* this instance of VB_Vector, in place. *
2443
*********************************************************************/
2444
VB_Vector& VB_Vector::operator+=(const double alpha)
2445
{
2446
2447
/*********************************************************************
2448
* gsl_vector_add_constant() is called to add the constant and its *
2449
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2450
* return value is non-zero, then a error occurred with the GSL *
2451
* function. VB_Vector::checkGSLStatus() will then process the error *
2452
* by writing out an appropriate error message and then causing this *
2453
* program to exit. *
2454
*********************************************************************/
2455
VB_Vector::checkGSLStatus(gsl_vector_add_constant(this->theVector, alpha),
2456
__LINE__, __FILE__, __FUNCTION__);
2457
2458
/*********************************************************************
2459
* Now returning a reference to this instance of VB_Vector. *
2460
*********************************************************************/
2461
return (*this);
2462
2463
} // VB_Vector& VB_Vector::operator+=(const double alpha)
2464
2465
/*********************************************************************
2466
* This overloaded operator subtracts the input scalar to all *
2467
* elements of this instance of VB_Vector, in place. *
2468
*********************************************************************/
2469
VB_Vector& VB_Vector::operator-=(const double alpha)
2470
{
2471
2472
/*********************************************************************
2473
* gsl_vector_add_constant() is called to subtract the constant and *
2474
* its return value is passed to VB_Vector::checkGSLStatus(). If the *
2475
* return value is non-zero, then a error occurred with the GSL *
2476
* function. VB_Vector::checkGSLStatus() will then process the error *
2477
* by writing out an appropriate error message and then causing this *
2478
* program to exit. *
2479
*********************************************************************/
2480
VB_Vector::checkGSLStatus(gsl_vector_add_constant(this->theVector, -1.0 * alpha),
2481
__LINE__, __FILE__, __FUNCTION__);
2482
2483
/*********************************************************************
2484
* Now returning a reference to this instance of VB_Vector. *
2485
*********************************************************************/
2486
return (*this);
2487
2488
} // VB_Vector& VB_Vector::operator-=(const double alpha)
2489
2490
/*********************************************************************
2491
* This overloaded operator multiplies each element of this instance *
2492
* of VB_Vector by the input scalar, in place. *
2493
*********************************************************************/
2494
VB_Vector& VB_Vector::operator*=(const double alpha)
2495
{
2496
2497
/*********************************************************************
2498
* gsl_vector_scale() is called to do the scaling and its *
2499
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2500
* return value is non-zero, then a error occurred with the GSL *
2501
* function. VB_Vector::checkGSLStatus() will then process the error *
2502
* by writing out an appropriate error message and then causing this *
2503
* program to exit. *
2504
*********************************************************************/
2505
VB_Vector::checkGSLStatus(gsl_vector_scale(this->theVector, alpha),
2506
__LINE__, __FILE__, __FUNCTION__);
2507
2508
/*********************************************************************
2509
* Now returning a reference to this instance of VB_Vector. *
2510
*********************************************************************/
2511
return (*this);
2512
2513
} // VB_Vector& VB_Vector::operator*=(const double alpha)
2514
2515
/*********************************************************************
2516
* This overloaded operator multiplies each element of this instance *
2517
* of VB_Vector by the reciprocal of the input scalar, in place. *
2518
*********************************************************************/
2519
VB_Vector& VB_Vector::operator/=(const double alpha)
2520
{
2521
2522
/*********************************************************************
2523
* If alpha is 0.0, then an exception is thrown. *
2524
*********************************************************************/
2525
if (alpha == 0.0)
2526
{
2527
2528
/*********************************************************************
2529
* Calling VB_Vector::createException() to throw and catch a *
2530
* GenericExcep so that the appropriate error messaqge can be printed.*
2531
*********************************************************************/
2532
VB_Vector::createException(string("Can not divide by a zero scalar value."),
2533
__LINE__, __FILE__, __FUNCTION__);
2534
2535
} // if
2536
2537
/*********************************************************************
2538
* Calling gsl_vector_scale() to carry out the scaling by the *
2539
* reciprocal. status will hold the return value of *
2540
* gsl_vector_scale(). *
2541
*********************************************************************/
2542
double reciprocal = 1.0 / alpha;
2543
2544
/*********************************************************************
2545
* gsl_vector_scale() is called to do the needed scaling and its *
2546
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2547
* return value is non-zero, then a error occurred with the GSL *
2548
* function. VB_Vector::checkGSLStatus() will then process the error *
2549
* by writing out an appropriate error message and then causing this *
2550
* program to exit. *
2551
*********************************************************************/
2552
VB_Vector::checkGSLStatus(gsl_vector_scale(this->theVector, reciprocal),
2553
__LINE__, __FILE__, __FUNCTION__);
2554
2555
/*********************************************************************
2556
* VB_Vector::checkFiniteness() is called to ensure that all elements *
2557
* of this->theVector are finite, i.e., no Inf's nor Nan's. *
2558
*********************************************************************/
2559
VB_Vector::checkFiniteness(this->theVector, __LINE__, __FILE__, __FUNCTION__);
2560
2561
2562
/*********************************************************************
2563
* Now returning a reference to this instance of VB_Vector. *
2564
*********************************************************************/
2565
return (*this);
2566
2567
} // VB_Vector& VB_Vector::operator/=(const double alpha)
2568
2569
/*********************************************************************
2570
* This overloaded operator adds the input vector to this instance of *
2571
* VB_Vector, in place. *
2572
*********************************************************************/
2573
VB_Vector& VB_Vector::operator+=(const VB_Vector& V)
2574
{
2575
2576
/*********************************************************************
2577
* gsl_vector_add() is called to compute the vector sum and its *
2578
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2579
* return value is non-zero, then a error occurred with the GSL *
2580
* function. VB_Vector::checkGSLStatus() will then process the error *
2581
* by writing out an appropriate error message and then causing this *
2582
* program to exit. *
2583
*********************************************************************/
2584
VB_Vector::checkGSLStatus(gsl_vector_add(this->theVector, V.theVector),
2585
__LINE__, __FILE__, __FUNCTION__);
2586
2587
/*********************************************************************
2588
* Now returning a reference to this instance of VB_Vector. *
2589
*********************************************************************/
2590
return (*this);
2591
2592
} // VB_Vector& VB_Vector::operator+=(const VB_Vector& V)
2593
2594
/*********************************************************************
2595
* This overloaded operator subtracts the input vector from this *
2596
* instance of VB_Vector, in place. *
2597
*********************************************************************/
2598
VB_Vector& VB_Vector::operator-=(const VB_Vector& V)
2599
{
2600
2601
/*********************************************************************
2602
* gsl_vector_sub() is called to compute the vector difference and its*
2603
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2604
* return value is non-zero, then a error occurred with the GSL *
2605
* function. VB_Vector::checkGSLStatus() will then process the error *
2606
* by writing out an appropriate error message and then causing this *
2607
* program to exit. *
2608
*********************************************************************/
2609
VB_Vector::checkGSLStatus(gsl_vector_sub(this->theVector, V.theVector),
2610
__LINE__, __FILE__, __FUNCTION__);
2611
2612
/*********************************************************************
2613
* Now returning a reference to this instance of VB_Vector. *
2614
*********************************************************************/
2615
return (*this);
2616
2617
} // VB_Vector& VB_Vector::operator-=(const VB_Vector& V)
2618
2619
/*********************************************************************
2620
* This overloaded operator multiplies each element of this instance *
2621
* of VB_Vector by the corresponding element of the input vector, in *
2622
* place. *
2623
*********************************************************************/
2624
VB_Vector& VB_Vector::operator*=(const VB_Vector& V)
2625
{
2626
2627
/*********************************************************************
2628
* gsl_vector_mul() is called to compute the element by element *
2629
* product and its return value is passed to *
2630
* VB_Vector::checkGSLStatus(). If the return value is non-zero, then *
2631
* a error occurred with the GSL function. VB_Vector::checkGSLStatus()*
2632
* will then process the error by writing out an appropriate error *
2633
* message and then causing this program to exit. *
2634
*********************************************************************/
2635
VB_Vector::checkGSLStatus(gsl_vector_mul(this->theVector, V.theVector),
2636
__LINE__, __FILE__, __FUNCTION__);
2637
2638
/*********************************************************************
2639
* Now returning a reference to this instance of VB_Vector. *
2640
*********************************************************************/
2641
return (*this);
2642
2643
} // VB_Vector& VB_Vector::operator*=(const VB_Vector& V)
2644
2645
/*********************************************************************
2646
* This overloaded operator divides each element of this instance *
2647
* of VB_Vector by the corresponding element of the input vector, in *
2648
* place. *
2649
*********************************************************************/
2650
VB_Vector& VB_Vector::operator/=(const VB_Vector& V)
2651
{
2652
2653
/*********************************************************************
2654
* gsl_vector_div() is called to do the element by element division *
2655
* and its return value is passed to VB_Vector::checkGSLStatus(). If *
2656
* the return value is non-zero, then a error occurred with the GSL *
2657
* function. VB_Vector::checkGSLStatus() will then process the error *
2658
* by writing out an appropriate error message and then causing this *
2659
* program to exit. *
2660
*********************************************************************/
2661
VB_Vector::checkGSLStatus(gsl_vector_div(this->theVector, V.theVector),
2662
__LINE__, __FILE__, __FUNCTION__);
2663
2664
/*********************************************************************
2665
* VB_Vector::checkFiniteness() is called to ensure that all elements *
2666
* of this->theVector are finite, i.e., no Inf's nor Nan's. *
2667
*********************************************************************/
2668
VB_Vector::checkFiniteness(this->theVector, __LINE__, __FILE__, __FUNCTION__);
2669
2670
/*********************************************************************
2671
* Now returning a reference to this instance of VB_Vector. *
2672
*********************************************************************/
2673
return (*this);
2674
2675
} // VB_Vector& VB_Vector::operator/=(const VB_Vector& V)
2676
2677
/*********************************************************************
2678
* This overloaded operator adds the input vector to this instance of *
2679
* VB_Vector, in place. *
2680
*********************************************************************/
2681
VB_Vector& VB_Vector::operator+=(const VB_Vector *V)
2682
{
2683
2684
/*********************************************************************
2685
* gsl_vector_add() is called to compute the vector sum and its *
2686
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2687
* return value is non-zero, then a error occurred with the GSL *
2688
* function. VB_Vector::checkGSLStatus() will then process the error *
2689
* by writing out an appropriate error message and then causing this *
2690
* program to exit. *
2691
*********************************************************************/
2692
VB_Vector::checkGSLStatus(gsl_vector_add(this->theVector, V->theVector),
2693
__LINE__, __FILE__, __FUNCTION__);
2694
2695
/*********************************************************************
2696
* Now returning a reference to this instance of VB_Vector. *
2697
*********************************************************************/
2698
return (*this);
2699
2700
} // VB_Vector& VB_Vector::operator+=(const VB_Vector *V)
2701
2702
/*********************************************************************
2703
* This overloaded operator subtracts the input vector from this *
2704
* instance of VB_Vector, in place. *
2705
*********************************************************************/
2706
VB_Vector& VB_Vector::operator-=(const VB_Vector *V)
2707
{
2708
2709
/*********************************************************************
2710
* gsl_vector_sub() is called to compute the vector difference and its*
2711
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2712
* return value is non-zero, then a error occurred with the GSL *
2713
* function. VB_Vector::checkGSLStatus() will then process the error *
2714
* by writing out an appropriate error message and then causing this *
2715
* program to exit. *
2716
*********************************************************************/
2717
VB_Vector::checkGSLStatus(gsl_vector_sub(this->theVector, V->theVector),
2718
__LINE__, __FILE__, __FUNCTION__);
2719
2720
/*********************************************************************
2721
* Now returning a reference to this instance of VB_Vector. *
2722
*********************************************************************/
2723
return (*this);
2724
2725
} // VB_Vector& VB_Vector::operator-=(const VB_Vector *V)
2726
2727
/*********************************************************************
2728
* This overloaded operator multiplies each element of this instance *
2729
* of VB_Vector by the corresponding element of the input vector, in *
2730
* place. *
2731
*********************************************************************/
2732
VB_Vector& VB_Vector::operator*=(const VB_Vector *V)
2733
{
2734
2735
/*********************************************************************
2736
* gsl_vector_mul() is called to compute the element by element *
2737
* product and its return value is passed to *
2738
* VB_Vector::checkGSLStatus(). If the return value is non-zero, then *
2739
* a error occurred with the GSL function. VB_Vector::checkGSLStatus()*
2740
* will then process the error by writing out an appropriate error *
2741
* message and then causing this program to exit. *
2742
*********************************************************************/
2743
VB_Vector::checkGSLStatus(gsl_vector_mul(this->theVector, V->theVector),
2744
__LINE__, __FILE__, __FUNCTION__);
2745
2746
/*********************************************************************
2747
* Now returning a reference to this instance of VB_Vector. *
2748
*********************************************************************/
2749
return (*this);
2750
2751
} // VB_Vector& VB_Vector::operator*=(const VB_Vector *V)
2752
2753
/*********************************************************************
2754
* This overloaded operator divides each element of this instance *
2755
* of VB_Vector by the corresponding element of the input vector, in *
2756
* place. *
2757
*********************************************************************/
2758
VB_Vector& VB_Vector::operator/=(const VB_Vector *V)
2759
{
2760
2761
/*********************************************************************
2762
* gsl_vector_div() is called to do the element by element division *
2763
* and its return value is passed to VB_Vector::checkGSLStatus(). If *
2764
* the return value is non-zero, then a error occurred with the GSL *
2765
* function. VB_Vector::checkGSLStatus() will then process the error *
2766
* by writing out an appropriate error message and then causing this *
2767
* program to exit. *
2768
*********************************************************************/
2769
VB_Vector::checkGSLStatus(gsl_vector_div(this->theVector, V->theVector),
2770
__LINE__, __FILE__, __FUNCTION__);
2771
2772
/*********************************************************************
2773
* VB_Vector::checkFiniteness() is called to ensure that all elements *
2774
* of this->theVector are finite, i.e., no Inf's nor Nan's. *
2775
*********************************************************************/
2776
VB_Vector::checkFiniteness(this->theVector, __LINE__, __FILE__, __FUNCTION__);
2777
2778
/*********************************************************************
2779
* Now returning a reference to this instance of VB_Vector. *
2780
*********************************************************************/
2781
return (*this);
2782
2783
} // VB_Vector& VB_Vector::operator/=(const VB_Vector *V)
2784
2785
/*********************************************************************
2786
* This overloaded operator adds the input vector to this instance of *
2787
* VB_Vector, in place. *
2788
*********************************************************************/
2789
VB_Vector& VB_Vector::operator+=(const gsl_vector *V)
2790
{
2791
2792
/*********************************************************************
2793
* gsl_vector_add() is called to compute the vector sum and its *
2794
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2795
* return value is non-zero, then a error occurred with the GSL *
2796
* function. VB_Vector::checkGSLStatus() will then process the error *
2797
* by writing out an appropriate error message and then causing this *
2798
* program to exit. *
2799
*********************************************************************/
2800
VB_Vector::checkGSLStatus(gsl_vector_add(this->theVector, V),
2801
__LINE__, __FILE__, __FUNCTION__);
2802
2803
/*********************************************************************
2804
* Now returning a reference to this instance of VB_Vector. *
2805
*********************************************************************/
2806
return (*this);
2807
2808
} // VB_Vector& VB_Vector::operator+=(const gsl_vector *V)
2809
2810
/*********************************************************************
2811
* This overloaded operator subtracts the input vector from this *
2812
* instance of VB_Vector, in place. *
2813
*********************************************************************/
2814
VB_Vector& VB_Vector::operator-=(const gsl_vector *V)
2815
{
2816
2817
/*********************************************************************
2818
* gsl_vector_sub() is called to compute the vector difference and its*
2819
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2820
* return value is non-zero, then a error occurred with the GSL *
2821
* function. VB_Vector::checkGSLStatus() will then process the error *
2822
* by writing out an appropriate error message and then causing this *
2823
* program to exit. *
2824
*********************************************************************/
2825
VB_Vector::checkGSLStatus(gsl_vector_sub(this->theVector, V),
2826
__LINE__, __FILE__, __FUNCTION__);
2827
2828
/*********************************************************************
2829
* Now returning a reference to this instance of VB_Vector. *
2830
*********************************************************************/
2831
return (*this);
2832
2833
} // VB_Vector& VB_Vector::operator-=(const gsl_vector *V)
2834
2835
/*********************************************************************
2836
* This overloaded operator multiplies each element of this instance *
2837
* of VB_Vector by the corresponding element of the input vector, in *
2838
* place. *
2839
*********************************************************************/
2840
VB_Vector& VB_Vector::operator*=(const gsl_vector *V)
2841
{
2842
2843
/*********************************************************************
2844
* gsl_vector_mul() is called to compute the element by element *
2845
* product and its return value is passed to *
2846
* VB_Vector::checkGSLStatus(). If the return value is non-zero, then *
2847
* a error occurred with the GSL function. VB_Vector::checkGSLStatus()*
2848
* will then process the error by writing out an appropriate error *
2849
* message and then causing this program to exit. *
2850
*********************************************************************/
2851
VB_Vector::checkGSLStatus(gsl_vector_mul(this->theVector, V),
2852
__LINE__, __FILE__, __FUNCTION__);
2853
2854
/*********************************************************************
2855
* Now returning a reference to this instance of VB_Vector. *
2856
*********************************************************************/
2857
return (*this);
2858
2859
} // VB_Vector& VB_Vector::operator*=(const gsl_vector *V)
2860
2861
/*********************************************************************
2862
* This overloaded operator divides each element of this instance *
2863
* of VB_Vector by the corresponding element of the input vector, in *
2864
* place. *
2865
*********************************************************************/
2866
VB_Vector& VB_Vector::operator/=(const gsl_vector *V)
2867
{
2868
2869
/*********************************************************************
2870
* gsl_vector_div() is called to do the element by element division *
2871
* and its return value is passed to VB_Vector::checkGSLStatus(). If *
2872
* the return value is non-zero, then a error occurred with the GSL *
2873
* function. VB_Vector::checkGSLStatus() will then process the error *
2874
* by writing out an appropriate error message and then causing this *
2875
* program to exit. *
2876
*********************************************************************/
2877
VB_Vector::checkGSLStatus(gsl_vector_div(this->theVector, V),
2878
__LINE__, __FILE__, __FUNCTION__);
2879
2880
/*********************************************************************
2881
* VB_Vector::checkFiniteness() is called to ensure that all elements *
2882
* of this->theVector are finite, i.e., no Inf's nor Nan's. *
2883
*********************************************************************/
2884
VB_Vector::checkFiniteness(this->theVector, __LINE__, __FILE__, __FUNCTION__);
2885
2886
/*********************************************************************
2887
* Now returning a reference to this instance of VB_Vector. *
2888
*********************************************************************/
2889
return (*this);
2890
2891
} // VB_Vector& VB_Vector::operator/=(const gsl_vector *V)
2892
2893
/*********************************************************************
2894
* This overloaded operator adds the input vector to this instance of *
2895
* VB_Vector, in place. *
2896
*********************************************************************/
2897
VB_Vector& VB_Vector::operator+=(const gsl_vector& V)
2898
{
2899
2900
/*********************************************************************
2901
* gsl_vector_add() is called to compute the vector sum and its *
2902
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2903
* return value is non-zero, then a error occurred with the GSL *
2904
* function. VB_Vector::checkGSLStatus() will then process the error *
2905
* by writing out an appropriate error message and then causing this *
2906
* program to exit. *
2907
*********************************************************************/
2908
VB_Vector::checkGSLStatus(gsl_vector_add(this->theVector, (&V)),
2909
__LINE__, __FILE__, __FUNCTION__);
2910
2911
/*********************************************************************
2912
* Now returning a reference to this instance of VB_Vector. *
2913
*********************************************************************/
2914
return (*this);
2915
2916
} // VB_Vector& VB_Vector::operator+=(const gsl_vector& V)
2917
2918
/*********************************************************************
2919
* This overloaded operator subtracts the input vector from this *
2920
* instance of VB_Vector, in place. *
2921
*********************************************************************/
2922
VB_Vector& VB_Vector::operator-=(const gsl_vector& V)
2923
{
2924
2925
/*********************************************************************
2926
* gsl_vector_sub() is called to compute the vector difference and its*
2927
* return value is passed to VB_Vector::checkGSLStatus(). If the *
2928
* return value is non-zero, then a error occurred with the GSL *
2929
* function. VB_Vector::checkGSLStatus() will then process the error *
2930
* by writing out an appropriate error message and then causing this *
2931
* program to exit. *
2932
*********************************************************************/
2933
VB_Vector::checkGSLStatus(gsl_vector_sub(this->theVector, (&V)),
2934
__LINE__, __FILE__, __FUNCTION__);
2935
2936
/*********************************************************************
2937
* Now returning a reference to this instance of VB_Vector. *
2938
*********************************************************************/
2939
return (*this);
2940
2941
} // VB_Vector& VB_Vector::operator-=(const gsl_vector& V)
2942
2943
/*********************************************************************
2944
* This overloaded operator multiplies each element of this instance *
2945
* of VB_Vector by the corresponding element of the input vector, in *
2946
* place. *
2947
*********************************************************************/
2948
VB_Vector& VB_Vector::operator*=(const gsl_vector& V)
2949
{
2950
2951
/*********************************************************************
2952
* gsl_vector_mul() is called to compute the element by element *
2953
* product and its return value is passed to *
2954
* VB_Vector::checkGSLStatus(). If the return value is non-zero, then *
2955
* a error occurred with the GSL function. VB_Vector::checkGSLStatus()*
2956
* will then process the error by writing out an appropriate error *
2957
* message and then causing this program to exit. *
2958
*********************************************************************/
2959
VB_Vector::checkGSLStatus(gsl_vector_mul(this->theVector, (&V)),
2960
__LINE__, __FILE__, __FUNCTION__);
2961
2962
/*********************************************************************
2963
* Now returning a reference to this instance of VB_Vector. *
2964
*********************************************************************/
2965
return (*this);
2966
2967
} // VB_Vector& VB_Vector::operator*=(const gsl_vector& V)
2968
2969
/*********************************************************************
2970
* This overloaded operator divides each element of this instance *
2971
* of VB_Vector by the corresponding element of the input vector, in *
2972
* place. *
2973
*********************************************************************/
2974
VB_Vector& VB_Vector::operator/=(const gsl_vector& V)
2975
{
2976
2977
/*********************************************************************
2978
* gsl_vector_div() is called to do the element by element division *
2979
* and its return value is passed to VB_Vector::checkGSLStatus(). If *
2980
* the return value is non-zero, then a error occurred with the GSL *
2981
* function. VB_Vector::checkGSLStatus() will then process the error *
2982
* by writing out an appropriate error message and then causing this *
2983
* program to exit. *
2984
*********************************************************************/
2985
VB_Vector::checkGSLStatus(gsl_vector_div(this->theVector, (&V)),
2986
__LINE__, __FILE__, __FUNCTION__);
2987
2988
/*********************************************************************
2989
* VB_Vector::checkFiniteness() is called to ensure that all elements *
2990
* of this->theVector are finite, i.e., no Inf's nor Nan's. *
2991
*********************************************************************/
2992
VB_Vector::checkFiniteness(this->theVector, __LINE__, __FILE__, __FUNCTION__);
2993
2994
/*********************************************************************
2995
* Now returning a reference to this instance of VB_Vector. *
2996
*********************************************************************/
2997
return (*this);
2998
2999
} // VB_Vector& VB_Vector::operator/=(const gsl_vector& V)
3000
3001
/*********************************************************************
3002
* This instance method populates the input array with the element *
3003
* values from this instance of VB_Vector. It is assumes that the *
3004
* input array, theArray[], has been properly sized, i.e., sufficient *
3005
* memory has been allocated to theArray[], before this method is *
3006
* called. *
3007
*********************************************************************/
3008
void VB_Vector::fetchData(double *theArray, const size_t length) const
3009
{
3010
3011
/*********************************************************************
3012
* VB_Vector::checkVectorLengths() is called to see if *
3013
* this->theVector->size equals the input variable length. If they are*
3014
* not equal, then VB_Vector::checkVectorLengths() will display an *
3015
* appropriate error and then cause this program to exit. NOTE: Even *
3016
* this->theVector->size equals length, this is no guarantee that *
3017
* theArray[] has been properly sized. *
3018
*********************************************************************/
3019
VB_Vector::checkVectorLengths(this->getLength(), length, __LINE__,
3020
__FILE__, __FUNCTION__);
3021
3022
/*********************************************************************
3023
* Now copying the elements of this instance of VB_Vector to *
3024
* theArray[]. *
3025
*********************************************************************/
3026
memcpy(theArray, this->theVector->data, length * sizeof(double));
3027
3028
} // void VB_Vector::fetchData(double *theArray, const size_t length) const
3029
3030
/*********************************************************************
3031
* This instance method transforms this instance of VB_Vector into a *
3032
* unit vector. *
3033
*********************************************************************/
3034
void VB_Vector::unitVector()
3035
{
3036
3037
/*********************************************************************
3038
* The Euclidean norm is stored in norm. *
3039
*********************************************************************/
3040
double norm = this->euclideanNorm();
3041
3042
/*********************************************************************
3043
* If norm is positive, then we go ahead and make the unit vector. *
3044
*********************************************************************/
3045
if (norm > 0.0)
3046
{
3047
3048
/*********************************************************************
3049
* we now save the reciprocal of norm in norm. *
3050
*********************************************************************/
3051
norm = 1.0 / norm;
3052
3053
/*********************************************************************
3054
* Now calling this->scaleInPlace() to turn this instance of VB_Vector*
3055
* into a unit vector. *
3056
*********************************************************************/
3057
this->scaleInPlace(norm);
3058
3059
} // if
3060
3061
/*********************************************************************
3062
* If program flow ends up here, then we have a zero vector. *
3063
* Therefore, an exception is thrown with an appropriate error *
3064
* message. *
3065
*********************************************************************/
3066
else
3067
{
3068
3069
/*********************************************************************
3070
* Now throwing the exception. *
3071
*********************************************************************/
3072
VB_Vector::zeroVector(norm, __LINE__, __FILE__, __FUNCTION__);
3073
3074
} // else
3075
3076
} // void VB_Vector::unitVector()
3077
3078
/*********************************************************************
3079
* This method resizes the vector to the desired new length and all *
3080
* of the vector elements will be set to zero. *
3081
*********************************************************************/
3082
void VB_Vector::resize(size_t newLength)
3083
{
3084
3085
/*********************************************************************
3086
* If this->theVector is NULL (no gsl_vector was allocated) or if the *
3087
* current size of this instance of VB_Vector is not equal to the *
3088
* new desired length, then a resizing is carried out. *
3089
*********************************************************************/
3090
if (this->theVector==NULL || (this->theVector->size != newLength))
3091
{
3092
3093
/*********************************************************************
3094
* Calling this->init() to do the resizing. All vector elements will *
3095
* be set to zero. try/catch blocks are used to process any exception *
3096
* thrown by this->init(). *
3097
*********************************************************************/
3098
try
3099
{
3100
this->init(newLength);
3101
} // try
3102
catch (GenericExcep &e)
3103
{
3104
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3105
} // catch
3106
} // if
3107
3108
/*********************************************************************
3109
* If program flow ends up here, then this instance of VB_Vector is *
3110
* already of the required size. Therefore, we simply overwrite *
3111
* this->theVector->data with zeroes. *
3112
*********************************************************************/
3113
else if (this->theVector->size == newLength)
3114
{
3115
memset(this->theVector->data, 0, sizeof(double) * newLength);
3116
} // else if
3117
3118
} // void VB_Vector::resize(size_t newLength)
3119
3120
/*********************************************************************
3121
* This instance method computes and returns the vector sum. *
3122
*********************************************************************/
3123
double VB_Vector::getVectorSum() const
3124
{
3125
double sum = 0.0;
3126
3127
for (unsigned long i = 0; i < this->getLength(); i++)
3128
{
3129
sum += this->theVector->data[i];
3130
} // for i
3131
3132
return sum;
3133
} // double VB_Vector::getVectorSum() const
3134
3135
/*********************************************************************
3136
* This instance method computes and returns the vector mean. *
3137
*********************************************************************/
3138
double VB_Vector::getVectorMean() const
3139
{
3140
3141
/*********************************************************************
3142
* Calling the template function getMean() to return the mean. *
3143
*********************************************************************/
3144
return getMean(this->theVector->data, this->getLength());
3145
3146
} // double VB_Vector::getVectorMean() const
3147
3148
/*********************************************************************
3149
* This static method returns the mean of the input gsl_vector. *
3150
*********************************************************************/
3151
double VB_Vector::getVectorMean(const gsl_vector *V)
3152
{
3153
3154
/*********************************************************************
3155
* Calling the template function getMean() to return the mean. *
3156
*********************************************************************/
3157
return getMean(V->data, V->size);
3158
3159
} // double VB_Vector::getVectorMean(const gsl_vector *V)
3160
3161
/*********************************************************************
3162
* This static method returns the mean of the input gsl_vector. *
3163
*********************************************************************/
3164
double VB_Vector::getVectorMean(const gsl_vector& V)
3165
{
3166
3167
/*********************************************************************
3168
* Returning the mean. *
3169
*********************************************************************/
3170
return VB_Vector::getVectorMean(&V);
3171
3172
} // double VB_Vector::getVectorMean(const gsl_vector& V)
3173
3174
/*********************************************************************
3175
* This static method returns the mean of the input VB_Vector. *
3176
*********************************************************************/
3177
double VB_Vector::getVectorMean(const VB_Vector& V)
3178
{
3179
3180
/*********************************************************************
3181
* Returning the mean. *
3182
*********************************************************************/
3183
return VB_Vector::getVectorMean(V.theVector);
3184
3185
} // double VB_Vector::getVectorMean(const VB_Vector& V)
3186
3187
/*********************************************************************
3188
* This static method returns the mean of the input VB_Vector. *
3189
*********************************************************************/
3190
double VB_Vector::getVectorMean(const VB_Vector *V)
3191
{
3192
3193
/*********************************************************************
3194
* Returning the mean. *
3195
*********************************************************************/
3196
return VB_Vector::getVectorMean(V->theVector);
3197
3198
} // double VB_Vector::getVectorMean(const VB_Vector *V)
3199
3200
/*********************************************************************
3201
* This method convolves this instance of VB_Vector with the input *
3202
* VB_Vector, in place. This method implements the MATLAB style *
3203
* convolution function. *
3204
*********************************************************************/
3205
void VB_Vector::convolve(const VB_Vector& v)
3206
{
3207
this->convolve(v.theVector);
3208
} // void VB_Vector::convolve(const VB_Vector& v)
3209
3210
/*********************************************************************
3211
* This method convolves this instance of VB_Vector with the input *
3212
* VB_Vector, in place. This method implements the MATLAB style *
3213
* convolution function. *
3214
*********************************************************************/
3215
void VB_Vector::convolve(const VB_Vector *v)
3216
{
3217
this->convolve(*v);
3218
} // void VB_Vector::convolve(const VB_Vector *v)
3219
3220
/*********************************************************************
3221
* This method convolves this instance of VB_Vector with the input *
3222
* gsl_vector, in place. *
3223
* *
3224
* Here is the definition of the convolution used: *
3225
* 1. Let the resulting vector of the convolution of 2 vectors A *
3226
* and B be denoted by: (A, B). *
3227
* 2. Let Az denote the vector A concatenated with a zero vector *
3228
* of the same length as A. Similarly for B. *
3229
* 3. Let the length of A be m and the length of B be n. *
3230
* 4. Let [i] represent the i'th element of a vector. *
3231
* 5. Then (A, B)[k] is given by the sum: *
3232
* *
3233
* *
3234
* k *
3235
* ---- *
3236
* \ *
3237
* \ Az[j] * Bz[k - j] *
3238
* / *
3239
* /___ *
3240
* j = 0 *
3241
* *
3242
* for k = 0, 1, ..., m + n - 2. *
3243
* *
3244
* NOTE: Some texts allow the vector resulting from the convolution *
3245
* to have length (m + n) and not length (m + n - 1) as we do here. *
3246
* However, the last element in this vector of length (m + n) is *
3247
* always zero. This method simply omits this last element that is *
3248
* always zero. This method implements the MATLAB style convolution *
3249
* function. *
3250
*********************************************************************/
3251
void VB_Vector::convolve(const gsl_vector *kernel)
3252
{
3253
3254
/*********************************************************************
3255
* Since this method overwrites this instance of VB_Vector with the *
3256
* result of the convolution, we need to save this instance of *
3257
* VB_Vector. *
3258
*********************************************************************/
3259
VB_Vector tempVec(this);
3260
3261
/*********************************************************************
3262
* We now resize this instance of VB_Vector to have the required size *
3263
* of the result of the convolution. try/catch blocks are used in *
3264
* case this->init() throws an exception. *
3265
*********************************************************************/
3266
try
3267
{
3268
this->init(this->getLength() + kernel->size - 1);
3269
} // try
3270
catch (GenericExcep &e)
3271
{
3272
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3273
} // catch
3274
3275
/*********************************************************************
3276
* The following for loop indexes the vector resulting from the *
3277
* convolution. *
3278
*********************************************************************/
3279
for (size_t k = 0; k < this->getLength(); k++)
3280
{
3281
3282
/*********************************************************************
3283
* The following for loop is used to compute the k'th element of the *
3284
* convoluted vector. *
3285
*********************************************************************/
3286
for (size_t j = 0; j <= k; j++)
3287
{
3288
3289
/*********************************************************************
3290
* Recall that the definition of vector convolution used in this *
3291
* method involves concatenating each vector with a zero vector *
3292
* whose length is equal to the length of the vector itself. In order *
3293
* to keep things as efficient as possible, the concatenation of the *
3294
* zero vector was skipped, but we "assume" that each vector was *
3295
* concatenated with a zero vector of appropriate length. However, *
3296
* since neither of the 2 vectors was actually concatenated with a *
3297
* zero vector, we need to check to see if we are actually accessing *
3298
* a valid element in each vector or not. Specifically, when *
3299
* j >= tempVec.theVector->size, tempVec[j] is not a valid in *
3300
* tempVec. Had we actually appeneded a zero vector to tempVec, *
3301
* tempVec[j] (for j >= tempVec.theVector->size and *
3302
* j < (2 * tempVec.theVector->size)) would be a valid element of *
3303
* tempVec and it would equal zero. Similarly for the input vector *
3304
* kernel when (k - j) >= kernel->size and *
3305
* (k - j) < (2 * kernel->size). Therefore, the following if block is *
3306
* is used to ensure that only valid elements of tempVec and kernel *
3307
* are accessed. Doing things this way allows us to skip the step of *
3308
* appending a zero vector to each of the 2 vectors whose convolution *
3309
* we calculate here. *
3310
*********************************************************************/
3311
if ( (j < tempVec.getLength()) && ((k - j) < kernel->size) )
3312
{
3313
(*this)[k] += tempVec[j] * kernel->data[k - j];
3314
} // if
3315
} // for j
3316
} // for k
3317
3318
} // void VB_Vector::convolve(const gsl_vector *kernel)
3319
3320
/*********************************************************************
3321
* This method convolves this instance of VB_Vector with the input *
3322
* gsl_vector, in place. This method implements the MATLAB style *
3323
* convolution function. *
3324
*********************************************************************/
3325
void VB_Vector::convolve(const gsl_vector& v)
3326
{
3327
this->convolve(&v);
3328
} // void VB_Vector::convolve(const gsl_vector& v)
3329
3330
/*********************************************************************
3331
* This method convolves this instance of VB_Vector with the input *
3332
* VB_Vector and returns the convolved vector. This method implements *
3333
* the MATLAB style convolution function. *
3334
*********************************************************************/
3335
VB_Vector VB_Vector::convolve2(const VB_Vector& v) const
3336
{
3337
3338
/*********************************************************************
3339
* A copy of this instance of VB_Vector is saved to vTemp. Then vTemp *
3340
* is convolved with the input VB_Vector. Finally, vTemp is returned. *
3341
*********************************************************************/
3342
VB_Vector vTemp(this); // NOTE: vTemp.fileName is set to this->fileName.
3343
vTemp.convolve(v.theVector);
3344
return vTemp;
3345
3346
} // VB_Vector VB_Vector::convolve2(const VB_Vector& v) const
3347
3348
/*********************************************************************
3349
* This method convolves this instance of VB_Vector with the input *
3350
* VB_Vector and returns the convolved vector. This method implements *
3351
* the MATLAB style convolution function. *
3352
*********************************************************************/
3353
VB_Vector VB_Vector::convolve2(const VB_Vector *v) const
3354
{
3355
3356
/*********************************************************************
3357
* A copy of this instance of VB_Vector is saved to vTemp. Then vTemp *
3358
* is convolved with the input VB_Vector. Finally, vTemp is returned. *
3359
*********************************************************************/
3360
VB_Vector vTemp(this); // NOTE: vTemp.fileName is set to this->fileName.
3361
vTemp.convolve(v->theVector);
3362
return vTemp;
3363
3364
} // VB_Vector VB_Vector::convolve2(const VB_Vector *v) const
3365
3366
/*********************************************************************
3367
* This method convolves this instance of VB_Vector with the input *
3368
* gsl_vector and returns the convolved vector. This method implements*
3369
* the MATLAB style convolution function. *
3370
*********************************************************************/
3371
VB_Vector VB_Vector::convolve2(const gsl_vector& v) const
3372
{
3373
3374
/*********************************************************************
3375
* A copy of this instance of VB_Vector is saved to vTemp. Then vTemp *
3376
* is convolved with the input gsl_vector. Finally, vTemp is returned.*
3377
*********************************************************************/
3378
VB_Vector vTemp(this); // NOTE: vTemp.fileName is set to this->fileName.
3379
vTemp.convolve(&v);
3380
return vTemp;
3381
3382
} // VB_Vector VB_Vector::convolve2(const gsl_vector& v) const
3383
3384
/*********************************************************************
3385
* This method convolves this instance of VB_Vector with the input *
3386
* gsl_vector and returns the convolved vector. This method implements*
3387
* the MATLAB style convolution function. *
3388
*********************************************************************/
3389
VB_Vector VB_Vector::convolve2(const gsl_vector *v) const
3390
{
3391
3392
/*********************************************************************
3393
* A copy of this instance of VB_Vector is saved to vTemp. Then vTemp *
3394
* is convolved with the input gsl_vector. Finally, vTemp is returned.*
3395
*********************************************************************/
3396
VB_Vector v2(this); // NOTE: vTemp.fileName is set to this->fileName.
3397
v2.convolve(v);
3398
return v2;
3399
3400
} // VB_Vector VB_Vector::convolve2(const gsl_vector *v) const
3401
3402
/*********************************************************************
3403
* This is a static method that convolves the 2 input gsl_vector *
3404
* variables and returns a VB_Vector. This method implements the *
3405
* MATLAB style convolution function. *
3406
*********************************************************************/
3407
VB_Vector VB_Vector::convolve(const gsl_vector *v1, const gsl_vector *v2)
3408
{
3409
3410
/*********************************************************************
3411
* vTemp is instantiated from v1 and then convolved with v2. Finally, *
3412
* vTemp is returned. *
3413
*********************************************************************/
3414
VB_Vector vTemp(v1); // NOTE: vTemp.fileName is empty.
3415
vTemp.convolve(v2);
3416
return vTemp;
3417
3418
} // VB_Vector VB_Vector::convolve(const gsl_vector *v1, const gsl_vector *v2)
3419
3420
/*********************************************************************
3421
* This is a static method that convolves the 2 input gsl_vector *
3422
* variables and returns a VB_Vector. This method implements the *
3423
* MATLAB style convolution function. *
3424
*********************************************************************/
3425
VB_Vector VB_Vector::convolve(const gsl_vector& v1, const gsl_vector& v2)
3426
{
3427
3428
/*********************************************************************
3429
* Returning the convolution. *
3430
*********************************************************************/
3431
return VB_Vector::convolve(&v1, &v2);
3432
3433
} // VB_Vector VB_Vector::convolve(const gsl_vector& v1, const gsl_vector& v2)
3434
3435
/*********************************************************************
3436
* This is a static method that convolves the 2 input VB_Vector *
3437
* variables and returns a VB_Vector. This method implements the *
3438
* MATLAB style convolution function. *
3439
*********************************************************************/
3440
VB_Vector VB_Vector::convolve(const VB_Vector& v1, const VB_Vector& v2)
3441
{
3442
3443
/*********************************************************************
3444
* Returning the convolution. *
3445
*********************************************************************/
3446
return VB_Vector::convolve(v1.theVector, v2.theVector);
3447
3448
} // VB_Vector VB_Vector::convolve(const VB_Vector& v1, const VB_Vector& v2)
3449
3450
/*********************************************************************
3451
* This is a static method that convolves the 2 input VB_Vector *
3452
* variables and returns a VB_Vector. This method implements the *
3453
* MATLAB style convolution function. *
3454
*********************************************************************/
3455
VB_Vector VB_Vector::convolve(const VB_Vector *v1, const VB_Vector *v2)
3456
{
3457
3458
/*********************************************************************
3459
* Returning the convolution. *
3460
*********************************************************************/
3461
return VB_Vector::convolve(v1->theVector, v2->theVector);
3462
3463
} // VB_Vector VB_Vector::convolve(const VB_Vector *v1, const VB_Vector *v2)
3464
3465
/*********************************************************************
3466
* This instance method computes the variance of the vector. *
3467
*********************************************************************/
3468
double VB_Vector::getVariance() const
3469
{
3470
3471
/*********************************************************************
3472
* variance will hold the variance of the vector. It's initialized to *
3473
* zero. *
3474
*********************************************************************/
3475
double variance = 0.0;
3476
3477
/*********************************************************************
3478
* mean is assigned the vector mean. *
3479
*********************************************************************/
3480
double mean = this->getVectorMean();
3481
3482
/*********************************************************************
3483
* This for loop is used to compute the sum of the square of the *
3484
* difference between each vector element value and the mean value. *
3485
*********************************************************************/
3486
for (size_t i = 0; i < this->theVector->size; i++)
3487
{
3488
variance += ((*this)[i] - mean) * ((*this)[i] - mean);
3489
} // for
3490
3491
/*********************************************************************
3492
* Now returning the variance. *
3493
*********************************************************************/
3494
return variance / (double ) (this->theVector->size - 1);
3495
3496
} // double VB_Vector::getVariance() const
3497
3498
/*********************************************************************
3499
* This instance method makes this instance of VB_Vector have unit *
3500
* variance. *
3501
*********************************************************************/
3502
void VB_Vector::unitVariance()
3503
{
3504
3505
/*********************************************************************
3506
* stdDev is assigned the standard deviation of the vector. *
3507
*********************************************************************/
3508
double stdDev = sqrt(this->getVariance());
3509
3510
/*********************************************************************
3511
* Variance is computed by using the following formula: *
3512
* *
3513
* n - 1 *
3514
* ---- *
3515
* \ _ *
3516
* sigma^2 = \ x(i) - x *
3517
* / ---------------- *
3518
* /___ (n - 1) *
3519
* i = 0 *
3520
* *
3521
* where n = the vector size. *
3522
* *
3523
* If the standard deviation is not 0, then the following for loop *
3524
* will normalize this->theVector to have unit variance. *
3525
*********************************************************************/
3526
if (stdDev != 0.0)
3527
{
3528
(*this) /= stdDev;
3529
} // if
3530
3531
} // void VB_Vector::unitVariance()
3532
3533
/*********************************************************************
3534
* This instance method returns a VB_Vector with unit variance derived*
3535
* from this instance of VB_Vector. *
3536
*********************************************************************/
3537
VB_Vector VB_Vector::unitVarianceVector() const
3538
{
3539
3540
/*********************************************************************
3541
* This instance of VB_Vector is assigned to vTemp. Then vTemp is *
3542
* made to have unit variance, after which it is returned. *
3543
*********************************************************************/
3544
VB_Vector vTemp(this); // NOTE: vTemp.fileName is set to this->fileName.
3545
vTemp.unitVariance();
3546
return vTemp;
3547
3548
} // VB_Vector VB_Vector::unitVarianceVector() const
3549
3550
/*********************************************************************
3551
* This instance method appends the input VB_Vector to this instance *
3552
* of VB_Vector, in place. *
3553
*********************************************************************/
3554
void VB_Vector::concatenate(const VB_Vector& V)
3555
{
3556
3557
/*********************************************************************
3558
* Calling this->concatenate(const gsl_vector *) to carry out the *
3559
* concatenation operation. *
3560
*********************************************************************/
3561
this->concatenate(V.theVector);
3562
3563
} // void VB_Vector::concatenate(const VB_Vector& V)
3564
3565
/*********************************************************************
3566
* This instance method appends the input VB_Vector to this instance *
3567
* of VB_Vector, in place. *
3568
*********************************************************************/
3569
void VB_Vector::concatenate(const VB_Vector *V)
3570
{
3571
3572
/*********************************************************************
3573
* Calling this->concatenate(const gsl_vector *) to carry out the *
3574
* concatenation operation. *
3575
*********************************************************************/
3576
this->concatenate(V->theVector);
3577
3578
} // void VB_Vector::concatenate(const VB_Vector *V)
3579
3580
/*********************************************************************
3581
* This instance method appends the input gsl_vector to this instance *
3582
* of VB_Vector, in place. *
3583
*********************************************************************/
3584
void VB_Vector::concatenate(const gsl_vector& V)
3585
{
3586
3587
/*********************************************************************
3588
* Calling this->concatenate(const gsl_vector *) to carry out the *
3589
* concatenation operation. *
3590
*********************************************************************/
3591
this->concatenate(&V);
3592
3593
} // void VB_Vector::concatenate(const gsl_vector& V)
3594
3595
/*********************************************************************
3596
* This instance method appends the input gsl_vector to this instance *
3597
* of VB_Vector, in place. *
3598
*********************************************************************/
3599
void VB_Vector::concatenate(const gsl_vector *V)
3600
{
3601
3602
3603
/*********************************************************************
3604
* If both the vector to be concatenated have positive length: *
3605
*********************************************************************/
3606
if ( (this->theVector != NULL) && (V != NULL) )
3607
{
3608
3609
/*********************************************************************
3610
* By default, if the memory allocation fails, *
3611
* then gsl_vector_calloc() will invoke the GSL error handler *
3612
* which calls abort(), creating a core dump. To forgo this behavior, *
3613
* the GSL error handler is turned off and the current error handler *
3614
* is stored in VB_Vector::currentGSLErrorHandler. This is *
3615
* accomplished by the call to this->turnOffGSLErrorHandler(). *
3616
*********************************************************************/
3617
this->turnOffGSLErrorHandler();
3618
3619
/*********************************************************************
3620
* It's possible that the gsl_vector data member of this instance of *
3621
* VB_Vector has the same address as V, i.e., this instance of *
3622
* VB_Vector is being concatenated with itself, which we want to *
3623
* allow. Therefore, two temprary gsl_vector structs are allocated. *
3624
* this->theVector will be copied to tempVec and V will be copied to *
3625
* tempVec2. *
3626
*********************************************************************/
3627
gsl_vector *tempVec = gsl_vector_calloc(this->getLength());
3628
gsl_vector *tempVec2 = gsl_vector_calloc(V->size);
3629
3630
/*********************************************************************
3631
* We now restore the standard GSL error handler by calling *
3632
* this->restoreGSLErrorHandler(). *
3633
*********************************************************************/
3634
this->restoreGSLErrorHandler();
3635
3636
/*********************************************************************
3637
* Ensuring that tempVec is non-null. If tempVec is null, then the *
3638
* try/catch will handle the exception thrown by *
3639
* VB_Vector::vectorNull(), which will have an appropriate error *
3640
* message, as given by gsl_strerror(). *
3641
*********************************************************************/
3642
try
3643
{
3644
VB_Vector::vectorNull(tempVec);
3645
} // try
3646
catch (GenericExcep &e)
3647
{
3648
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3649
} // catch
3650
3651
/*********************************************************************
3652
* Ensuring that tempVec2 is non-null. If tempVec2 is null, then the *
3653
* try/catch will handle the exception thrown by *
3654
* VB_Vector::vectorNull(), which will have an appropriate error *
3655
* message, as given by gsl_strerror(). *
3656
*********************************************************************/
3657
try
3658
{
3659
VB_Vector::vectorNull(tempVec2);
3660
} // try
3661
catch (GenericExcep &e)
3662
{
3663
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3664
} // catch
3665
3666
/*********************************************************************
3667
* Copying this->theVector to tempVec. If VB_Vector::GSLVectorMemcpy()*
3668
* throws an exception (which will have an appropriate error message, *
3669
* as determined by gsl_strerror()) then it is be caught here. *
3670
*********************************************************************/
3671
try
3672
{
3673
VB_Vector::GSLVectorMemcpy(tempVec, this->theVector);
3674
} // try
3675
catch (GenericExcep &e)
3676
{
3677
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3678
} // catch
3679
3680
/*********************************************************************
3681
* Copying V to tempVec. If VB_Vector::GSLVectorMemcpy() throws an *
3682
* exception (which will have an appropriate error message, as *
3683
* determined by gsl_strerror()) then it is be caught here. *
3684
*********************************************************************/
3685
try
3686
{
3687
VB_Vector::GSLVectorMemcpy(tempVec2, V);
3688
} // try
3689
catch (GenericExcep &e)
3690
{
3691
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3692
} // catch
3693
3694
/*********************************************************************
3695
* Resizing this->theVector to (this->theVector->size + V->size). If *
3696
* an exception occurrs, which will have an appropriate error message *
3697
* as determined by gsl_strerror(), then it is caught here. *
3698
*********************************************************************/
3699
try
3700
{
3701
this->init(this->getLength() + V->size);
3702
} // try
3703
catch (GenericExcep &e)
3704
{
3705
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3706
} // catch
3707
3708
/*********************************************************************
3709
* Now assigning the values from tempVec to this instance of *
3710
* VB_Vector. *
3711
*********************************************************************/
3712
memcpy(this->theVector->data, tempVec->data, sizeof(double) * tempVec->size);
3713
3714
/*********************************************************************
3715
* Now assigning the values from tempVec2 to this instance of *
3716
* VB_Vector. *
3717
*********************************************************************/
3718
memcpy(this->theVector->data + tempVec->size, tempVec2->data, sizeof(double) * tempVec2->size);
3719
3720
/*********************************************************************
3721
* Freeing the memory allocated to tempVec and tempVec2. *
3722
*********************************************************************/
3723
gsl_vector_free(tempVec);
3724
gsl_vector_free(tempVec2);
3725
3726
} // if
3727
3728
/*********************************************************************
3729
* If V is NULL, then we do nothing since there's no point in *
3730
* concatenating an empty vector to this instance of VB_Vector. *
3731
*********************************************************************/
3732
else if (V == NULL)
3733
{
3734
}
3735
3736
/*********************************************************************
3737
* If this instance of VB_Vector is empty, i.e., has zero length, or *
3738
* V is empty. *
3739
*********************************************************************/
3740
else if (this->theVector==NULL) // qq Could this->theVector ever be NULL ?
3741
{
3742
3743
/*********************************************************************
3744
* By default, if the memory allocation fails, *
3745
* then gsl_vector_calloc() will invoke the GSL error handler *
3746
* which calls abort(), creating a core dump. To forgo this behavior, *
3747
* the GSL error handler is turned off and the current error handler *
3748
* is stored in VB_Vector::currentGSLErrorHandler. This is *
3749
* accomplished by the call to this->turnOffGSLErrorHandler(). *
3750
*********************************************************************/
3751
this->turnOffGSLErrorHandler();
3752
3753
/*********************************************************************
3754
* Allocating sufficient space for this->theVector. *
3755
*********************************************************************/
3756
this->theVector = gsl_vector_calloc(V->size);
3757
3758
/*********************************************************************
3759
* Ensuring that this->theVector is non-null. If this->theVector is *
3760
* NULL, then the try/catch will handle the exception thrown by *
3761
* VB_Vector::vectorNull(), which will have an appropriate error *
3762
* message, as given by gsl_strerror(). *
3763
*********************************************************************/
3764
try
3765
{
3766
VB_Vector::vectorNull(this->theVector);
3767
} // try
3768
catch (GenericExcep &e)
3769
{
3770
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3771
} // catch
3772
3773
/*********************************************************************
3774
* We now restore the standard GSL error handler by calling *
3775
* this->restoreGSLErrorHandler(). *
3776
*********************************************************************/
3777
this->restoreGSLErrorHandler();
3778
3779
/*********************************************************************
3780
* Copying V to tempVec. If VB_Vector::GSLVectorMemcpy() throws an *
3781
* exception (which will have an appropriate error message, as *
3782
* determined by gsl_strerror()) then it is be caught here. *
3783
*********************************************************************/
3784
try
3785
{
3786
VB_Vector::GSLVectorMemcpy(this->theVector, V);
3787
} // try
3788
catch (GenericExcep &e)
3789
{
3790
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
3791
} // catch
3792
3793
/*********************************************************************
3794
* At this point, we have successfully copied V to this->theVector. *
3795
* Therefore, this->valid is set to true. *
3796
*********************************************************************/
3797
this->valid = true;
3798
3799
} // else if
3800
3801
} // void VB_Vector::concatenate(const gsl_vector *V)
3802
3803
/*********************************************************************
3804
* This instance method appends the input VB_Vector to this instance *
3805
* of VB_Vector, not in place, and returns the concatenated VB_Vector.*
3806
*********************************************************************/
3807
VB_Vector VB_Vector::concatenate2(const VB_Vector& V) const
3808
{
3809
3810
/*********************************************************************
3811
* This instance of VB_Vector is saved to vTemp. *
3812
*********************************************************************/
3813
VB_Vector vTemp(this); // NOTE: vTemp.fileName is set to this->fileName.
3814
3815
/*********************************************************************
3816
* Setting vTemp.fileName to the empty string. *
3817
*********************************************************************/
3818
vTemp.fileName = string("");
3819
3820
/*********************************************************************
3821
* Setting this->valid, this->dataType, and this->fileFormat. init() *
3822
* will also assign a new file name to vTemp.fileName. *
3823
*********************************************************************/
3824
vTemp.init(false, vb_double, "ref1");
3825
3826
/*********************************************************************
3827
* The input vector is concatenated to vTemp and then vTemp is *
3828
* returned. *
3829
*********************************************************************/
3830
vTemp.concatenate(V);
3831
return vTemp;
3832
3833
} // VB_Vector VB_Vector::concatenate2(const VB_Vector& V) const
3834
3835
/*********************************************************************
3836
* This instance method appends the input VB_Vector to this instance *
3837
* of VB_Vector, not in place, and returns the concatenated VB_Vector.*
3838
*********************************************************************/
3839
VB_Vector VB_Vector::concatenate2(const VB_Vector *V) const
3840
{
3841
3842
/*********************************************************************
3843
* This instance of VB_Vector is saved to vTemp. *
3844
*********************************************************************/
3845
VB_Vector vTemp(this); // NOTE: vTemp.fileName is set to this->fileName.
3846
3847
/*********************************************************************
3848
* Setting vTemp.fileName to the empty string. *
3849
*********************************************************************/
3850
vTemp.fileName = string("");
3851
3852
/*********************************************************************
3853
* Setting this->valid, this->dataType, and this->fileFormat. init() *
3854
* will also assign a new file name to vTemp.fileName. *
3855
*********************************************************************/
3856
vTemp.init(false, vb_double, "ref1");
3857
3858
/*********************************************************************
3859
* The input vector is concatenated to vTemp and then vTemp is *
3860
* returned. *
3861
*********************************************************************/
3862
vTemp.concatenate(V);
3863
return vTemp;
3864
3865
} // VB_Vector VB_Vector::concatenate2(const VB_Vector *V) const
3866
3867
/*********************************************************************
3868
* This instance method appends the input gsl_vector to this instance *
3869
* of VB_Vector, not in place, and returns the concatenated VB_Vector.*
3870
*********************************************************************/
3871
VB_Vector VB_Vector::concatenate2(const gsl_vector *V) const
3872
{
3873
3874
/*********************************************************************
3875
* This instance of VB_Vector is saved to vTemp. *
3876
*********************************************************************/
3877
VB_Vector vTemp(this); // NOTE: vTemp.fileName is set to this->fileName.
3878
3879
/*********************************************************************
3880
* Setting vTemp.fileName to the empty string. *
3881
*********************************************************************/
3882
vTemp.fileName = string("");
3883
3884
/*********************************************************************
3885
* Setting this->valid, this->dataType, and this->fileFormat. init() *
3886
* will also assign a new file name to vTemp.fileName. *
3887
*********************************************************************/
3888
vTemp.init(false, vb_double, "ref1");
3889
3890
/*********************************************************************
3891
* The input vector is concatenated to vTemp and then vTemp is *
3892
* returned. *
3893
*********************************************************************/
3894
vTemp.concatenate(V);
3895
return vTemp;
3896
3897
} // VB_Vector VB_Vector::concatenate2(const gsl_vector *V) const
3898
3899
/*********************************************************************
3900
* This instance method appends the input gsl_vector to this instance *
3901
* of VB_Vector, not in place, and returns the concatenated VB_Vector.*
3902
*********************************************************************/
3903
VB_Vector VB_Vector::concatenate2(const gsl_vector& V) const
3904
{
3905
3906
/*********************************************************************
3907
* This instance of VB_Vector is saved to vTemp. *
3908
*********************************************************************/
3909
VB_Vector vTemp(this); // NOTE: vTemp.fileName is set to this->fileName.
3910
3911
/*********************************************************************
3912
* Setting vTemp.fileName to the empty string. *
3913
*********************************************************************/
3914
vTemp.fileName = string("");
3915
3916
/*********************************************************************
3917
* Setting this->valid, this->dataType, and this->fileFormat. init() *
3918
* will also assign a new file name to vTemp.fileName. *
3919
*********************************************************************/
3920
vTemp.init(false, vb_double, "ref1");
3921
3922
/*********************************************************************
3923
* The input vector is concatenated to vTemp and then vTemp is *
3924
* returned. *
3925
*********************************************************************/
3926
vTemp.concatenate(V);
3927
return vTemp;
3928
3929
} // VB_Vector VB_Vector::concatenate2(const gsl_vector& V) const
3930
3931
/*********************************************************************
3932
* This is a static method that concatenates the 2 input vectors and *
3933
* returns the concatenated vector. *
3934
*********************************************************************/
3935
VB_Vector VB_Vector::concatenate(const VB_Vector& V1, const VB_Vector& V2)
3936
{
3937
3938
/*********************************************************************
3939
* V1 is saved to vTemp. *
3940
*********************************************************************/
3941
VB_Vector vTemp(V1); // NOTE: vTemp.fileName is set to V1.fileName.
3942
3943
/*********************************************************************
3944
* Setting vTemp.fileName to the empty string. *
3945
*********************************************************************/
3946
vTemp.fileName = string("");
3947
3948
/*********************************************************************
3949
* Setting this->valid, this->dataType, and this->fileFormat. init() *
3950
* will also assign a new file name to vTemp.fileName. *
3951
*********************************************************************/
3952
vTemp.init(false, vb_double, "ref1");
3953
3954
/*********************************************************************
3955
* The input vector is concatenated to vTemp and then vTemp is *
3956
* returned. *
3957
*********************************************************************/
3958
vTemp.concatenate(V2);
3959
return vTemp;
3960
3961
} // VB_Vector VB_Vector::concatenate(const VB_Vector& V1, const VB_Vector& V2)
3962
3963
/*********************************************************************
3964
* This is a static method that concatenates the 2 input vectors and *
3965
* returns the concatenated vector. *
3966
*********************************************************************/
3967
VB_Vector VB_Vector::concatenate(const VB_Vector *V1, const VB_Vector *V2)
3968
{
3969
3970
/*********************************************************************
3971
* V1 is saved to vTemp. *
3972
*********************************************************************/
3973
VB_Vector vTemp(V1); // NOTE: vTemp.fileName is set to V1->fileName.
3974
3975
/*********************************************************************
3976
* Setting vTemp.fileName to the empty string. *
3977
*********************************************************************/
3978
vTemp.fileName = string("");
3979
3980
/*********************************************************************
3981
* Setting this->valid, this->dataType, and this->fileFormat. init() *
3982
* will also assign a new file name to vTemp.fileName. *
3983
*********************************************************************/
3984
vTemp.init(false, vb_double, "ref1");
3985
3986
/*********************************************************************
3987
* The input vector is concatenated to vTemp and then vTemp is *
3988
* returned. *
3989
*********************************************************************/
3990
vTemp.concatenate(V2);
3991
return vTemp;
3992
3993
} // VB_Vector VB_Vector::concatenate(const VB_Vector *V1, const VB_Vector *V2)
3994
3995
/*********************************************************************
3996
* This is a static method that concatenates the 2 input vectors and *
3997
* returns the concatenated vector. *
3998
*********************************************************************/
3999
VB_Vector VB_Vector::concatenate(const gsl_vector& V1, const gsl_vector& V2)
4000
{
4001
4002
/*********************************************************************
4003
* V1 is saved to vTemp. *
4004
*********************************************************************/
4005
VB_Vector vTemp(V1); // NOTE: vTemp.fileName is empty.
4006
4007
/*********************************************************************
4008
* Setting vTemp.fileName to the empty string. *
4009
*********************************************************************/
4010
vTemp.fileName = string("");
4011
4012
/*********************************************************************
4013
* Setting this->valid, this->dataType, and this->fileFormat. init() *
4014
* will also assign a new file name to vTemp.fileName. *
4015
*********************************************************************/
4016
vTemp.init(false, vb_double, "ref1");
4017
4018
/*********************************************************************
4019
* The input vector is concatenated to vTemp and then vTemp is *
4020
* returned. *
4021
*********************************************************************/
4022
vTemp.concatenate(V2);
4023
return vTemp;
4024
4025
} // VB_Vector VB_Vector::concatenate(const gsl_vector& V1, const gsl_vector& V2)
4026
4027
/*********************************************************************
4028
* This is a static method that concatenates the 2 input vectors and *
4029
* returns the concatenated vector. *
4030
*********************************************************************/
4031
VB_Vector VB_Vector::concatenate(const gsl_vector *V1, const gsl_vector *V2)
4032
{
4033
4034
/*********************************************************************
4035
* V1 is saved to vTemp. Then vTemp is concatenated with V1. Finally, *
4036
* vtemp is returned. *
4037
*********************************************************************/
4038
VB_Vector vTemp(V1); // NOTE: vTemp.fileName is empty.
4039
4040
/*********************************************************************
4041
* Setting vTemp.fileName to the empty string. *
4042
*********************************************************************/
4043
vTemp.fileName = string("");
4044
4045
/*********************************************************************
4046
* Setting this->valid, this->dataType, and this->fileFormat. init() *
4047
* will also assign a new file name to vTemp.fileName. *
4048
*********************************************************************/
4049
vTemp.init(false, vb_double, "ref1");
4050
4051
/*********************************************************************
4052
* The input vector is concatenated to vTemp and then vTemp is *
4053
* returned. *
4054
*********************************************************************/
4055
vTemp.concatenate(V2);
4056
return vTemp;
4057
4058
} // VB_Vector VB_Vector::concatenate(const gsl_vector *V1, const gsl_vector *V2)
4059
4060
/*********************************************************************
4061
* This is the overloaded left shift operator. The last i elements *
4062
* from the right are set to zero. *
4063
*********************************************************************/
4064
VB_Vector& VB_Vector::operator<<(size_t i)
4065
{
4066
4067
/*********************************************************************
4068
* If i is zero, i.e., no shift, then (*this) is simply returned. *
4069
*********************************************************************/
4070
if (i == 0)
4071
{
4072
return (*this);
4073
} // if
4074
4075
/*********************************************************************
4076
* If i equals or exceeds the length of this vector, then the vector *
4077
* is simply set to all zeros. *
4078
*********************************************************************/
4079
if (i >= this->getLength())
4080
{
4081
this->init(this->getLength());
4082
return (*this);
4083
} // if
4084
4085
/*********************************************************************
4086
* This for loop translates the elements of this instance of VB_Vector*
4087
* i spots to the left. *
4088
*********************************************************************/
4089
for (size_t j = 0; j < this->getLength() - i; j++)
4090
{
4091
(*this)[j] = (*this)[j + i];
4092
} // for j
4093
4094
/*********************************************************************
4095
* This for loop sets the last i right hand side elements of this *
4096
* instance of VB_Vector to zero. *
4097
*********************************************************************/
4098
for (size_t j = this->getLength() - i; j < this->getLength(); j++)
4099
{
4100
(*this)[j] = 0.0;
4101
} // for j
4102
4103
/*********************************************************************
4104
* Now returning this instance of VB_Vector. *
4105
*********************************************************************/
4106
return (*this);
4107
4108
} // VB_Vector& VB_Vector::operator<<(size_t i)
4109
4110
/*********************************************************************
4111
* This is the overloaded right shift operator. The first i elements *
4112
* from the left are set to zero. *
4113
*********************************************************************/
4114
VB_Vector& VB_Vector::operator>>(size_t i)
4115
{
4116
4117
/*********************************************************************
4118
* If i is zero, i.e., no shift, then (*this) is simply returned. *
4119
*********************************************************************/
4120
if (i == 0)
4121
{
4122
return (*this);
4123
} // if
4124
4125
/*********************************************************************
4126
* If i equals or exceeds the length of this vector, then the vector *
4127
* is simply set to all zeros. *
4128
*********************************************************************/
4129
if (i >= this->getLength())
4130
{
4131
this->init(this->getLength());
4132
return (*this);
4133
} // if
4134
4135
/*********************************************************************
4136
* This for loop translates the elements of this instance of VB_Vector*
4137
* i spots to the right. First we need to make a copy of this *
4138
* instance of VB_Vector since its elements get overwritten from the *
4139
* left. *
4140
*********************************************************************/
4141
VB_Vector tempVec(this);
4142
for (size_t j = i; j < this->getLength(); j++)
4143
{
4144
(*this)[j] = tempVec[j - i];
4145
} // for j
4146
4147
/*********************************************************************
4148
* This for loop sets the first i left hand side elements of this *
4149
* instance of VB_Vector to zero. *
4150
*********************************************************************/
4151
for (size_t j = 0; j < i; j++)
4152
{
4153
(*this)[j] = 0.0;
4154
} // for j
4155
4156
/*********************************************************************
4157
* Now returning this instance of VB_Vector. *
4158
*********************************************************************/
4159
return (*this);
4160
4161
} // VB_Vector& VB_Vector::operator>>(size_t i)
4162
4163
/*********************************************************************
4164
* This overloaed operator scales this instance of VB_Vector by the *
4165
* input alpha and returns the scaled vector. *
4166
*********************************************************************/
4167
VB_Vector VB_Vector::operator*(const double alpha) const
4168
{
4169
4170
/*********************************************************************
4171
* A copy is made of this instance of VB_Vector. *
4172
*********************************************************************/
4173
VB_Vector tempVec(this);
4174
4175
/*********************************************************************
4176
* gsl_vector_scale() is called to do the scaling and its *
4177
* return value is passed to VB_Vector::checkGSLStatus(). If the *
4178
* return value is non-zero, then a error occurred with the GSL *
4179
* function. VB_Vector::checkGSLStatus() will then process the error *
4180
* by writing out an appropriate error message and then causing this *
4181
* program to exit. *
4182
*********************************************************************/
4183
VB_Vector::checkGSLStatus(gsl_vector_scale(tempVec.theVector, alpha),
4184
__LINE__, __FILE__, __FUNCTION__);
4185
4186
/*********************************************************************
4187
* Now returning the scaled vector. *
4188
*********************************************************************/
4189
return tempVec;
4190
4191
} // VB_Vector VB_Vector::operator*(const double alpha) const
4192
4193
/*********************************************************************
4194
* This friend function scales the input VB_Vector by the input *
4195
* double. This function is implemented to preserve commutativity. *
4196
*********************************************************************/
4197
VB_Vector operator*(const double alpha, const VB_Vector& V)
4198
{
4199
4200
/*********************************************************************
4201
* The copy constructor is used to make a copy of V. *
4202
*********************************************************************/
4203
VB_Vector tempVec(V);
4204
4205
/*********************************************************************
4206
* gsl_vector_scale() is called to do the scaling and its *
4207
* return value is passed to VB_Vector::checkGSLStatus(). If the *
4208
* return value is non-zero, then a error occurred with the GSL *
4209
* function. VB_Vector::checkGSLStatus() will then process the error *
4210
* by writing out an appropriate error message and then causing this *
4211
* program to exit. *
4212
*********************************************************************/
4213
VB_Vector::checkGSLStatus(gsl_vector_scale(tempVec.theVector, alpha),
4214
__LINE__, __FILE__, __FUNCTION__);
4215
/*********************************************************************
4216
* tempVec is scaled by alpha and then returned. *
4217
*********************************************************************/
4218
/* tempVec.scaleInPlace(alpha);*/
4219
return tempVec;
4220
4221
} // VB_Vector operator*(const double alpha, const VB_Vector& V)
4222
4223
4224
void
4225
VB_Vector::meanCenter()
4226
{
4227
(*this) -= this->getVectorMean();
4228
}
4229
4230
/*********************************************************************
4231
* This method resets the elements of this instance of VB_Vector from *
4232
* the input array of doubles. The input variable length is assumed *
4233
* to be the length of the array theData. *
4234
*********************************************************************/
4235
void VB_Vector::setData(const double *theData, const size_t length)
4236
{
4237
4238
/*********************************************************************
4239
* We call this->init() to resize this instance of VB_Vector and *
4240
* catch its exception, if one is thrown. *
4241
*********************************************************************/
4242
try
4243
{
4244
this->init(length);
4245
} // try
4246
catch (GenericExcep &e)
4247
{
4248
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
4249
} // catch
4250
4251
memcpy(this->theVector->data, theData, sizeof(double) * length);
4252
4253
} // void VB_Vector::setData(const double *theData, const size_t length)
4254
4255
/*********************************************************************
4256
* This method resets the data values of this instance of VB_Vector *
4257
* from the data values of the input gsl_vector. *
4258
*********************************************************************/
4259
void VB_Vector::setData(const gsl_vector *v)
4260
{
4261
4262
/*********************************************************************
4263
* Calling this->setData(const double *, const size_t) to reset the *
4264
* values. *
4265
*********************************************************************/
4266
this->setData(v->data, v->size);
4267
4268
} // void VB_Vector::setData(const gsl_vector *v)
4269
4270
/*********************************************************************
4271
* This method reset the data values of this instance of VB_Vector *
4272
* from the data values of the input gsl_vector. *
4273
*********************************************************************/
4274
void VB_Vector::setData(const gsl_vector& v)
4275
{
4276
4277
/*********************************************************************
4278
* Calling this->setData(const double *, const size_t) to reset the *
4279
* values. *
4280
*********************************************************************/
4281
this->setData(v.data, v.size);
4282
4283
} // void VB_Vector::setData(const gsl_vector& v)
4284
4285
/*********************************************************************
4286
* This method reset the data values of this instance of VB_Vector *
4287
* from the data values of the input VB_Vector. *
4288
*********************************************************************/
4289
void VB_Vector::setData(const VB_Vector *V)
4290
{
4291
4292
/*********************************************************************
4293
* Calling this->setData(const double *, const size_t) to reset the *
4294
* values. *
4295
*********************************************************************/
4296
this->setData(V->theVector->data, V->getLength());
4297
4298
} // void VB_Vector::setData(const VB_Vector *V)
4299
4300
/*********************************************************************
4301
* This method reset the data values of this instance of VB_Vector *
4302
* from the data values of the input VB_Vector. *
4303
*********************************************************************/
4304
void VB_Vector::setData(const VB_Vector& V)
4305
{
4306
4307
/*********************************************************************
4308
* Calling this->setData(const double *, const size_t) to reset the *
4309
* values. *
4310
*********************************************************************/
4311
this->setData(V.theVector->data, V.getLength());
4312
4313
} // void VB_Vector::setData(const VB_Vector& V)
4314
4315
/*********************************************************************
4316
* This method multiplies this instance of VB_Vector by the *
4317
* orthogonal projector: *
4318
* *
4319
* T -1 T *
4320
* (I - A (A A) A ) *
4321
* *
4322
* where: *
4323
* *
4324
* A is the matrix created from the VB_Vector's found in the *
4325
* argument reference. *
4326
* I is an appropriately dimensioned identity matrix. *
4327
* *
4328
* Basically, if this instance of VB_Vector is denoted by x, then *
4329
* the following computation is carried out: *
4330
* *
4331
* T -1 T *
4332
* x' = (I - A (A A) A )x *
4333
* *
4334
* where x' is this instance of VB_Vector after all the computations *
4335
* have been carried out. *
4336
* *
4337
* To implement the computations in an efficient manner as possible, *
4338
* the following steps are carried out (assume A is (m x n) and that *
4339
* x is (m x 1)): *
4340
* *
4341
* T *
4342
* 1. Let v = A x (v is then (n x 1)) *
4343
* T -1 *
4344
* ==> x' = x - A (A A) v *
4345
* *
4346
* T -1 *
4347
* 2. Let y = (A A) v (y is (n x 1)) *
4348
* This linear system will be solved for y (this saves us the *
4349
* computationally expensive step of computing a matrix inverse). *
4350
* *
4351
* ==> x' = x - Ay *
4352
* *
4353
* NOTE: It is expected that A has linearly independent columns. *
4354
*********************************************************************/
4355
void VB_Vector::orthogonalize(const vector<VB_Vector> reference)
4356
{
4357
4358
/*********************************************************************
4359
* We now see if the input variable reference represents an *
4360
* overdetermined system or not. If it does, then an exception is *
4361
* thrown (by calling VB_Vector::createException()) and then caught. *
4362
* An error message is assembled first. *
4363
*********************************************************************/
4364
if (this->getLength() < reference.size())
4365
{
4366
4367
/*********************************************************************
4368
* errorMsg[] will hold the error message that the index is out of *
4369
* range. *
4370
*********************************************************************/
4371
char errorMsg[OPT_STRING_LENGTH];
4372
memset(errorMsg, 0, OPT_STRING_LENGTH);
4373
4374
/*********************************************************************
4375
* Populating errorMsg[] with the appropriate error message. *
4376
*********************************************************************/
4377
sprintf(errorMsg, "The vector length [%d] is < the number of column vectors [%d] (overdetermined system).",
4378
(int)(this->getLength()), (int)(reference.size()));
4379
4380
/*********************************************************************
4381
* Now calling VB_Vector::createException() to throw and catch the *
4382
* exception. *
4383
*********************************************************************/
4384
VB_Vector::createException(errorMsg, __LINE__, __FILE__, __FUNCTION__);
4385
4386
} // if
4387
4388
/*********************************************************************
4389
* By default, if the memory allocation fails, *
4390
* then gsl_matrix_calloc() (called by this->initMatrix()) will invoke*
4391
* the GSL error handler which calls abort(), creating a core dump. *
4392
* To forgo this behavior, the GSL error handler is turned off and *
4393
* the current error handler is stored in *
4394
* VB_Vector::currentGSLErrorHandler. This is accomplished by the *
4395
* call to this->turnOffGSLErrorHandler(). *
4396
*********************************************************************/
4397
this->turnOffGSLErrorHandler();
4398
4399
/*********************************************************************
4400
* Allocating a gsl_matrix with rows reference[0].getLength() and *
4401
* columns reference.size(). All the elements in the gsl_matrix will *
4402
* be initialized to 0. try/catch blocks are used to process the *
4403
* exception that may be thrown by this->initMatrix(). *
4404
*********************************************************************/
4405
gsl_matrix *A = NULL;
4406
try
4407
{
4408
A = this->initMatrix(reference[0].getLength(), reference.size());
4409
} // try
4410
catch (GenericExcep &e)
4411
{
4412
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
4413
} // catch
4414
4415
/*********************************************************************
4416
* Now allocating a gsl_matrix struct correctly dimensioned to hold *
4417
* T *
4418
* the product of A and A. All its elements will be initialized to 0.*
4419
* try/catch blocks are used to process the exception that may be *
4420
* thrown by this->initMatrix(). *
4421
*********************************************************************/
4422
gsl_matrix *aTa = NULL;
4423
try
4424
{
4425
aTa = this->initMatrix(reference.size(), reference.size());
4426
} // try
4427
catch (GenericExcep &e)
4428
{
4429
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
4430
} // catch
4431
4432
/*********************************************************************
4433
* We now populate A with the VB_Vectors from reference. Each *
4434
* VB_Vector from reference is seen as a column vector of A. If *
4435
* gsl_matrix_set_col() returns a non-zero return value, then an *
4436
* error occurred, which is handled by checkGSLStatus(). *
4437
*********************************************************************/
4438
for (size_t i = 0; i < reference.size(); i++)
4439
{
4440
VB_Vector::checkGSLStatus(gsl_matrix_set_col(A, i, reference[i].getTheVector()), __LINE__, __FILE__, __FUNCTION__);
4441
} // for i
4442
4443
/*********************************************************************
4444
* T *
4445
* v will be used to hold the matrix-vector product A x. Since A is *
4446
* (m x n) and x is (m x 1), v needs to have length n, which is the *
4447
* number of columns in A, found in the field A->size2. *
4448
*********************************************************************/
4449
VB_Vector v(A->size2);
4450
4451
/*********************************************************************
4452
* Now carrying out the matrix-vector multiplication (Ax). The product*
4453
* will be stored in v. The prototype for the function *
4454
* gsl_blas_dgemv() is: *
4455
* *
4456
* int gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA, *
4457
* double alpha, *
4458
* const gsl_matrix * A, *
4459
* const gsl_vector * x, *
4460
* double beta, *
4461
* gsl_vector * y); *
4462
* *
4463
* and the actual operations carried out by gsl_blas_dgemv() are: *
4464
* *
4465
* y := alpha (op(A)) x + beta (y) *
4466
* *
4467
* where op(A) = *
4468
* T *
4469
* A if TransA is CblasTrans *
4470
* A if TransA is CblasNoTrans *
4471
* T *
4472
* In this specific case, we are simply computing v = (A x). If *
4473
* gsl_blas_dgemv() returns a non-zero value, then an error occurred, *
4474
* which gets handled by checkGSLStatus(). *
4475
*********************************************************************/
4476
VB_Vector::checkGSLStatus(gsl_blas_dgemv(CblasTrans, 1.0, A, this->theVector, 0.0, v.theVector),
4477
__LINE__, __FILE__, __FUNCTION__);
4478
4479
/*********************************************************************
4480
* T *
4481
* Now carrying out the matrix multiplication (A * A). The product *
4482
* will be stored in aTa. The prototype for the function *
4483
* gsl_blas_dgemm() is: *
4484
* *
4485
* int gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, *
4486
* CBLAS_TRANSPOSE_t TransB, *
4487
* double alpha, *
4488
* const gsl_matrix * A, *
4489
* const gsl_matrix * B, *
4490
* double beta, *
4491
* gsl_matrix * C); *
4492
* *
4493
* and the actual operations carried out by gsl_blas_dgemm() are: *
4494
* *
4495
* C := alpha * op(A) * op(B) + beta * C *
4496
* *
4497
* where op(A) = *
4498
* T *
4499
* A if TransA is CblasTrans *
4500
* A if TransA is CblasNoTrans *
4501
* Similarly for op(B). *
4502
* T *
4503
* In this specific case, we are simply computing aTa = A A. *
4504
*********************************************************************/
4505
VB_Vector::checkGSLStatus(gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, A, A, 0.0, aTa),
4506
__LINE__, __FILE__, __FUNCTION__);
4507
4508
/*********************************************************************
4509
* As the second step, we use gsl_linalg_cholesky_decomp() to compute *
4510
* Cholesky factorization of aTa (since is is expected that the *
4511
* columns of A are linearly independent, aTa will be positive *
4512
* definite). After gsl_linalg_cholesky_decomp(), aTa will hold the *
4513
* Cholesky factors. checkGSLStatus() is used to process any non-zero *
4514
* value, i.e., an error, returned by gsl_linalg_cholesky_decomp(). *
4515
*********************************************************************/
4516
VB_Vector::checkGSLStatus(gsl_linalg_cholesky_decomp(aTa),
4517
__LINE__, __FILE__, __FUNCTION__);
4518
4519
/*********************************************************************
4520
* T -1 *
4521
* We need to find the vector y such that y = (A A) v. We will do *
4522
* this by solving the linear system: *
4523
* T *
4524
* (A A)y = v *
4525
* This is why we previously computed the Cholesky factorization of *
4526
* aTa. *
4527
* *
4528
* y is initialized to be a vector of length aTa->size1 and all of *
4529
* its elements will be 0. *
4530
*********************************************************************/
4531
VB_Vector y(aTa->size1);
4532
4533
/*********************************************************************
4534
* Now gsl_linalg_cholesky_solve() is used to solve for y. If *
4535
* gsl_linalg_cholesky_solve() returns a non-zero value, then an error*
4536
* occurred, which gets handled by checkGSLStatus(). *
4537
*********************************************************************/
4538
VB_Vector::checkGSLStatus(gsl_linalg_cholesky_solve(aTa, v.theVector, y.theVector),
4539
__LINE__, __FILE__, __FUNCTION__);
4540
4541
/*********************************************************************
4542
* We now call gsl_blas_dgemv() to do the matrix-vector multiplication*
4543
* (Ay). This product will be contained in v.theVector (since we no *
4544
* longer need v, we can safely overwrite it). However, v needs to be *
4545
* resized to have length A->size1 (the number of rows in A). If *
4546
* gsl_blas_dgemv() returns a non-zero value, then an error occurred, *
4547
* which gets handled by checkGSLStatus(). *
4548
*********************************************************************/
4549
v.init(A->size1);
4550
VB_Vector::checkGSLStatus(gsl_blas_dgemv(CblasNoTrans, 1.0, A, y.theVector, 0.0, v.theVector),
4551
__LINE__, __FILE__, __FUNCTION__);
4552
4553
/*********************************************************************
4554
* Finally, we now subtract v from this instance of VB_Vector and *
4555
* store the difference in this instance of VB_Vector. *
4556
*********************************************************************/
4557
(*this) -= v;
4558
4559
/*********************************************************************
4560
* We now restore the standard GSL error handler by calling *
4561
* this->restoreGSLErrorHandler(). *
4562
*********************************************************************/
4563
this->restoreGSLErrorHandler();
4564
4565
/*********************************************************************
4566
* Freeing the previously allocated memory for A and aTa. *
4567
*********************************************************************/
4568
gsl_matrix_free(A);
4569
gsl_matrix_free(aTa);
4570
4571
} // void VB_Vector::orthogonalize(const vector<VB_Vector> reference)
4572
4573
/*********************************************************************
4574
* This method multiplies the input VB_Vector myVec by the orthogonal *
4575
* projector: *
4576
* *
4577
* T -1 T *
4578
* (I - A (A A) A ) *
4579
* *
4580
* wehre A is a matrix derived from the inout reference and returns *
4581
* the resulting vector. *
4582
*********************************************************************/
4583
VB_Vector VB_Vector::orthogonalize(const VB_Vector &myVec, const vector<VB_Vector> reference) const
4584
{
4585
4586
/*********************************************************************
4587
* Instantiating tempVec from myVec. *
4588
*********************************************************************/
4589
VB_Vector tempVec(myVec);
4590
4591
/*********************************************************************
4592
* Multiplying tempVec by the projector and then returning tempVec. *
4593
*********************************************************************/
4594
tempVec.orthogonalize(reference);
4595
return tempVec;
4596
4597
} // VB_Vector VB_Vector::orthogonalize(const VB_Vector &myVec, const vector<VB_Vector> reference)
4598
4599
/*********************************************************************
4600
* This static method writes out the input gsl_matrix struct to *
4601
* stdout. *
4602
*********************************************************************/
4603
void VB_Vector::printMatrix(const gsl_matrix *M)
4604
{
4605
/* gsl_matrix_fprintf(stdout, M, "%g");*/
4606
4607
for (size_t i = 0; i < M->size1; i++)
4608
{
4609
for (size_t j = 0; j < M->size2; j++)
4610
{
4611
if (j == 0)
4612
{
4613
cout << "[ ";
4614
} // if
4615
cout << gsl_matrix_get(M, i, j);
4616
if (j != (M->size2 - 1))
4617
{
4618
cout << ", ";
4619
} // if
4620
else
4621
{
4622
cout << " ]" << endl;
4623
} // else
4624
} // for j
4625
} // for i
4626
cout << endl;
4627
4628
} // void VB_Vector::printMatrix(const gsl_matrix *M)
4629
4630
/*********************************************************************
4631
* This static method writes out the input gsl_matrix struct to *
4632
* stdout. *
4633
*********************************************************************/
4634
void VB_Vector::printMatrix(const gsl_matrix& M)
4635
{
4636
VB_Vector::printMatrix(&M);
4637
} // void VB_Vector::printMatrix(const gsl_matrix& M)
4638
4639
/*********************************************************************
4640
* This method computes the FFT of this instance of VB_Vector and *
4641
* stores the real and imaginary parts into the 2 input VB_Vectors. *
4642
* *
4643
* INPUT VARIABLES: TYPE: DESCRIPTION: *
4644
* ---------------- ----- ------------ *
4645
* realPart VB_Vector * Will hold the real part of the *
4646
* computed FFT. *
4647
* imagPart VB_Vector * Will hold the imaginary part of *
4648
* the computed FFT. *
4649
* *
4650
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
4651
* ----------------- ----- ------------ *
4652
* None. *
4653
*********************************************************************/
4654
void VB_Vector::fft(VB_Vector *realPart, VB_Vector *imagPart) const
4655
{
4656
4657
/*********************************************************************
4658
* Ensuring that realPart and imagPart are properly sized. *
4659
*********************************************************************/
4660
if (this->getLength() != realPart->getLength())
4661
{
4662
realPart->resize(this->theVector->size);
4663
} // if
4664
if (this->getLength() != imagPart->getLength())
4665
{
4666
imagPart->resize(this->theVector->size);
4667
} // if
4668
4669
/*********************************************************************
4670
* halfSize will be (this->theVector->size / 2) if the length of this *
4671
* instance of VB_Vector is even. Otherwise halfSize will be *
4672
* ((this->theVector->size / 2) - 1). *
4673
*********************************************************************/
4674
const unsigned int halfSize = this->theVector->size / 2;
4675
4676
/*********************************************************************
4677
* even will be true if the length of this instance of VB_Vector is *
4678
* even, false otherwise. *
4679
*********************************************************************/
4680
const bool even = ( (halfSize * 2) == this->theVector->size );
4681
4682
/*********************************************************************
4683
* The data from this instance of VB_Vector is copied to FFT[]. *
4684
*********************************************************************/
4685
double FFT[this->theVector->size];
4686
memcpy(FFT, this->theVector->data, this->theVector->size * sizeof(double));
4687
4688
/*********************************************************************
4689
* The following 2 variables are needed by GSL's *
4690
* gsl_fft_real_transform() function. waveTable is a lookup table for *
4691
* sines and cosines while workSpace is a required work space. *
4692
*********************************************************************/
4693
gsl_fft_real_wavetable *waveTable = gsl_fft_real_wavetable_alloc(this->theVector->size);
4694
gsl_fft_real_workspace *workSpace = gsl_fft_real_workspace_alloc(this->theVector->size);
4695
4696
/*********************************************************************
4697
* If waveTable is NULL, then an exception is thrown and then caught *
4698
* by VB_Vector::createException() to inform the user of this error. *
4699
*********************************************************************/
4700
if (!waveTable)
4701
{
4702
VB_Vector::createException("Unable to allocate gsl_fft_real_wavetable.", __LINE__, __FILE__, __FUNCTION__);
4703
} // if
4704
4705
/*********************************************************************
4706
* If workSpace is NULL, then an exception is thrown and then caught *
4707
* by VB_Vector::createException() to inform the user of this error. *
4708
*********************************************************************/
4709
if (!workSpace)
4710
{
4711
VB_Vector::createException("Unable to allocate gsl_fft_real_workspace.", __LINE__, __FILE__, __FUNCTION__);
4712
} // if
4713
4714
/*********************************************************************
4715
* Now computing the FFT. If the call to gsl_fft_real_transform() is *
4716
* successful, then status will be 0. *
4717
*********************************************************************/
4718
int status = gsl_fft_real_transform(FFT, 1, this->theVector->size, waveTable, workSpace);
4719
4720
/*********************************************************************
4721
* If status is non-zero, then VB_Vector::createException() is called *
4722
* to throw and catch an exception to inform the user of the error. *
4723
* The error message will contain the reason for the error, as *
4724
* returned by gsl_strerror(status). *
4725
*********************************************************************/
4726
if (status)
4727
{
4728
VB_Vector::createException(string(gsl_strerror(status) + string(".")), __LINE__, __FILE__, __FUNCTION__);
4729
} // if
4730
4731
/*********************************************************************
4732
* We will need to scale FFT[] by the reciprocal of the length of *
4733
* this instance of VB_Vector to make this implementation of the FFT *
4734
* identical to IDL's. To that end, oneOverSize is used to hold the *
4735
* reciprocal of the length of this instance of VB_Vector. *
4736
*********************************************************************/
4737
const double oneOverSize = 1.0 / this->theVector->size;
4738
4739
/*********************************************************************
4740
* FFT[] now holds the discrete Fourier transform of this instance of *
4741
* VB_Vector. However, the storage scheme is not straight forward and *
4742
* is called "half complex." (It's advantage is that it takes the *
4743
* exact same mount of space as the origianl vector.) Here is a *
4744
* descpription of the the "half complex" storage scheme from p. 200 *
4745
* of the GSL Reference Manual: *
4746
* *
4747
* ... the half-complex transform of a real sequence is stored with *
4748
* frequencies in increasing order, starting at zero, with the real *
4749
* and imaginary parts of each frequency in neighboring locations. *
4750
* When a value is know to be real, its imaginary part is not stored. *
4751
* *
4752
* Here are 2 examples, the first is for n = 5 (odd case), the second *
4753
* for n = 6 (even case): *
4754
* *
4755
* complex[0].real = halfcomplex[0] *
4756
* complex[0].imag = 0 (always 0) *
4757
* complex[1].real = halfcomplex[1] *
4758
* complex[1].imag = halfcomplex[2] *
4759
* complex[2].real = halfcomplex[3] *
4760
* complex[2].imag = halfcomplex[4] *
4761
* complex[3].real = halfcomplex[3] *
4762
* complex[3].imag = -halfcomplex[4] *
4763
* complex[4].real = halfcomplex[1] *
4764
* complex[4].imag = -halfcomplex[2] *
4765
* *
4766
* complex[0].real = halfcomplex[0] *
4767
* complex[0].imag = 0 (always 0) *
4768
* complex[1].real = halfcomplex[1] *
4769
* complex[1].imag = halfcomplex[2] *
4770
* complex[2].real = halfcomplex[3] *
4771
* complex[2].imag = halfcomplex[4] *
4772
* complex[3].real = halfcomplex[5] *
4773
* complex[3].imag = 0 (always 0 in even case) *
4774
* complex[4].real = halfcomplex[3] *
4775
* complex[4].imag = -halfcomplex[4] *
4776
* complex[5].real = halfcomplex[1] *
4777
* complex[5].imag = -halfcomplex[2] *
4778
* *
4779
* We will have to "unpack" this half-complex storage to populate *
4780
* realPart[] and imagPart[] with the real and imaginary parts of *
4781
* FFT[]. We now set realPart[0] and imagPart[0]. *
4782
*********************************************************************/
4783
(*realPart)[0] = FFT[0] * oneOverSize;
4784
(*imagPart)[0] = 0.0;
4785
4786
/*********************************************************************
4787
* The following for loop is used to populate the the reaminder of *
4788
* realPart[] and imagPart[] from FFT[]. This for loop also scales *
4789
* the entries of FFT[] by oneOverSize. *
4790
*********************************************************************/
4791
for (unsigned int i = 1; i < this->theVector->size; i++)
4792
{
4793
if (i < halfSize)
4794
{
4795
(*realPart)[i] = FFT[(2 * i) - 1] * oneOverSize;
4796
(*imagPart)[i] = FFT[(2 * i)] * oneOverSize;
4797
} // if
4798
else if (i == halfSize)
4799
{
4800
if (even)
4801
{
4802
(*realPart)[i] = FFT[this->theVector->size - 1] * oneOverSize;
4803
(*imagPart)[i] = 0.0;
4804
} // if
4805
else
4806
{
4807
(*realPart)[i] = FFT[this->theVector->size - 2] * oneOverSize;
4808
(*imagPart)[i] = FFT[this->theVector->size - 1] * oneOverSize;
4809
} // if
4810
} // else if
4811
else
4812
{
4813
4814
/*********************************************************************
4815
* NOTE: Since the index (this->theVector->size - i) is < i, *
4816
* (*realPart)[this->theVector->size - i] has already been scaled by *
4817
* oneOverSize. Similarly for (*imagPart)[this->theVector->size - i]. *
4818
*********************************************************************/
4819
(*realPart)[i] = (*realPart)[this->theVector->size - i];
4820
(*imagPart)[i] = -1.0 * (*imagPart)[this->theVector->size - i];
4821
4822
} // else
4823
4824
} // for i
4825
4826
/*********************************************************************
4827
* Now freeing waveTable and workSpace. *
4828
*********************************************************************/
4829
gsl_fft_real_wavetable_free(waveTable);
4830
gsl_fft_real_workspace_free(workSpace);
4831
4832
} // void VB_Vector::fft(VB_Vector *realPart, VB_Vector *imagPart) const
4833
4834
/*********************************************************************
4835
* This method computes the FFT of this instance of VB_Vector and *
4836
* stores the real and imaginary parts into the 2 input VB_Vectors. *
4837
* *
4838
* INPUT VARIABLES: TYPE: DESCRIPTION: *
4839
* ---------------- ----- ------------ *
4840
* realPart VB_Vector& Will hold the real part of the *
4841
* computed FFT. *
4842
* imagPart VB_Vector& Will hold the imaginary part of *
4843
* the computed FFT. *
4844
* *
4845
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
4846
* ----------------- ----- ------------ *
4847
* None. *
4848
*********************************************************************/
4849
void VB_Vector::fft(VB_Vector& realPart, VB_Vector& imagPart) const
4850
{
4851
4852
/*********************************************************************
4853
* We now call this->fft(VB_Vector *realPart, VB_Vector *imagPart). *
4854
*********************************************************************/
4855
this->fft(&realPart, &imagPart);
4856
4857
} // void VB_Vector::fft(VB_Vector& realPart, VB_Vector& imagPart) const
4858
4859
/*********************************************************************
4860
* This method computes the inverse FFT of this instance of VB_Vector *
4861
* and stores the real and imaginary parts into the 2 input *
4862
* VB_Vectors. *
4863
* *
4864
* INPUT VARIABLES: TYPE: DESCRIPTION: *
4865
* ---------------- ----- ------------ *
4866
* realPart VB_Vector * Will hold the real part of the *
4867
* computed inverse FFT. *
4868
* imagPart VB_Vector * Will hold the imaginary part of *
4869
* the computed inverse FFT. *
4870
* *
4871
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
4872
* ----------------- ----- ------------ *
4873
* None. *
4874
*********************************************************************/
4875
void VB_Vector::ifft(VB_Vector *realPart, VB_Vector *imagPart) const
4876
{
4877
4878
/*********************************************************************
4879
* Ensuring that realPart and imagPart are appropriately sized. *
4880
*********************************************************************/
4881
if (this->getLength() != realPart->getLength())
4882
{
4883
realPart->resize(this->theVector->size);
4884
} // if
4885
if (this->getLength() != imagPart->getLength())
4886
{
4887
imagPart->resize(this->theVector->size);
4888
} // if
4889
4890
/*********************************************************************
4891
* status will be used to hold the return status values from the GSL *
4892
* functions. *
4893
*********************************************************************/
4894
int status = 0;
4895
4896
/*********************************************************************
4897
* complexVec[] will be a complex representation of the elements of *
4898
* this instance of VB_Vector. Of course, the imaginary parts will *
4899
* all be zero. The real parts will be in the even numbered indexes *
4900
* of complexVec[] and the imaginary parts will be in the odd *
4901
* numbered indexes of complexVec[]. *
4902
*********************************************************************/
4903
double complexVec[this->theVector->size * 2];
4904
4905
/*********************************************************************
4906
* Calling gsl_fft_real_unpack() to populate complexVec[]. *
4907
*********************************************************************/
4908
status = gsl_fft_real_unpack(this->theVector->data, (gsl_complex_packed_array ) complexVec, 1, this->theVector->size);
4909
4910
/*********************************************************************
4911
* If status is non-zero, then VB_Vector::createException() is called *
4912
* to throw and catch an exception to inform the user of the error. *
4913
* The error message will contain the reason for the error, as *
4914
* returned by gsl_strerror(status). *
4915
*********************************************************************/
4916
if (status)
4917
{
4918
VB_Vector::createException(string(gsl_strerror(status) + string(".")), __LINE__, __FILE__, __FUNCTION__);
4919
} // if
4920
4921
/*********************************************************************
4922
* The following 2 variables are needed by GSL's *
4923
* gsl_fft_complex_backward() function. waveTable is a lookup table *
4924
* for sines and cosines while workSpace is a required work space. *
4925
*********************************************************************/
4926
gsl_fft_complex_wavetable *waveTable = gsl_fft_complex_wavetable_alloc(this->theVector->size);
4927
gsl_fft_complex_workspace *workSpace = gsl_fft_complex_workspace_alloc(this->theVector->size);
4928
4929
/*********************************************************************
4930
* If waveTable is NULL, then an exception is thrown and then caught *
4931
* by VB_Vector::createException() to inform the user of this error. *
4932
*********************************************************************/
4933
if (!waveTable)
4934
{
4935
VB_Vector::createException("Unable to allocate gsl_fft_complex_wavetable.", __LINE__, __FILE__, __FUNCTION__);
4936
} // if
4937
4938
/*********************************************************************
4939
* If workSpace is NULL, then an exception is thrown and then caught *
4940
* by VB_Vector::createException() to inform the user of this error. *
4941
*********************************************************************/
4942
if (!workSpace)
4943
{
4944
VB_Vector::createException("Unable to allocate gsl_fft_complex_workspace.", __LINE__, __FILE__, __FUNCTION__);
4945
} // if
4946
4947
/*********************************************************************
4948
* Now computing the inverse FFT. *
4949
*********************************************************************/
4950
status = gsl_fft_complex_backward(complexVec, 1, this->theVector->size, waveTable, workSpace);
4951
4952
/*********************************************************************
4953
* If status is non-zero, then VB_Vector::createException() is called *
4954
* to throw and catch an exception to inform the user of the error. *
4955
* The error message will contain the reason for the error, as *
4956
* returned by gsl_strerror(status). *
4957
*********************************************************************/
4958
if (status)
4959
{
4960
VB_Vector::createException(string(gsl_strerror(status) + string(".")), __LINE__, __FILE__, __FUNCTION__);
4961
} // if
4962
4963
/*********************************************************************
4964
* At this point, complexVec[] contains the inverse FFT of this *
4965
* instance of VB_Vector; its even indexed elements holds the real *
4966
* parts and its odd indexed elements holds the imaginary part. Now *
4967
* we populate realPart and imagPart appropriately with the following *
4968
* for loop. i indexes complexVec[] and j indexes realPart and *
4969
* imagPart. *
4970
*********************************************************************/
4971
for (unsigned int i = 0, j = 0; j < this->theVector->size; j++)
4972
{
4973
4974
/*********************************************************************
4975
* Assigning the real part of the the j'th complex number from *
4976
* complexVec[] to (*realPart)[j]. Note that i is even and i++ causes *
4977
* i to be odd. *
4978
*********************************************************************/
4979
realPart->theVector->data[j] = complexVec[i++];
4980
4981
/*********************************************************************
4982
* Assigning the imaginary part of the the j'th complex number from *
4983
* complexVec[] to (*imagPart)[j]. Note that i is odd and i++ causes *
4984
* i to be even. *
4985
*********************************************************************/
4986
imagPart->theVector->data[j] = complexVec[i++];
4987
4988
} // for i, j
4989
4990
/*********************************************************************
4991
* Now freeing waveTable and workSpace. *
4992
*********************************************************************/
4993
gsl_fft_complex_wavetable_free(waveTable);
4994
gsl_fft_complex_workspace_free(workSpace);
4995
4996
} // void VB_Vector::ifft(VB_Vector *realPart, VB_Vector *imagPart) const
4997
4998
/*********************************************************************
4999
* This method computes the inverse FFT of this instance of VB_Vector *
5000
* and stores the real and imaginary parts into the 2 input *
5001
* VB_Vectors. *
5002
* *
5003
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5004
* ---------------- ----- ------------ *
5005
* realPart VB_Vector& Will hold the real part of the *
5006
* computed inverse FFT. *
5007
* imagPart VB_Vector& Will hold the imaginary part of *
5008
* the computed inverse FFT. *
5009
* *
5010
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5011
* ----------------- ----- ------------ *
5012
* None. *
5013
*********************************************************************/
5014
void VB_Vector::ifft(VB_Vector& realPart, VB_Vector& imagPart) const
5015
{
5016
5017
/*********************************************************************
5018
* We now call this->ifft(VB_Vector *realPart, VB_Vector *imagPart). *
5019
*********************************************************************/
5020
this->ifft(&realPart, &imagPart);
5021
5022
} // void VB_Vector::ifft(VB_Vector& realPart, VB_Vector& imagPart) const
5023
5024
/*********************************************************************
5025
* Instance method computes the power spectrum of this instance of *
5026
* VB_Vector. This method is derived from VoxBo_Fourier.pro's *
5027
* ReturnPS function. The power spectrum will be stored in result. *
5028
*********************************************************************/
5029
void VB_Vector::getPS(VB_Vector &result) const
5030
{
5031
5032
/*********************************************************************
5033
* Since we will need the FFT of this instance of VB_Vector, realPart *
5034
* is needed to hold the real part of the FFT and imagPart is need to *
5035
* hold the imaginary part of the FFT. *
5036
*********************************************************************/
5037
VB_Vector realPart(this->getLength());
5038
VB_Vector imagPart(this->getLength());
5039
5040
/*********************************************************************
5041
* Ensuring that result is of the appropriate size. *
5042
*********************************************************************/
5043
if (this->theVector->size != result.theVector->size)
5044
{
5045
result.resize(this->theVector->size);
5046
} // if
5047
5048
/*********************************************************************
5049
* Geting the FFT of this instance of VB_Vector. *
5050
*********************************************************************/
5051
this->fft(realPart, imagPart);
5052
5053
/*********************************************************************
5054
* The following for loop is used to assign the power spectrum values *
5055
* to result. The power spectrum values are computed from: *
5056
* *
5057
* fft(this) * conjugate(fft(this)) *
5058
* *
5059
* where fft(this) is the FFT of this instance of VB_Vector, "*" *
5060
* denotes element-by-element multiplication of the 2 VB_Vector's, *
5061
* and conjuagte() is the complex comjugate of a VB_Vector. *
5062
*********************************************************************/
5063
for (unsigned long i = 0; i < this->theVector->size; i++)
5064
{
5065
result[i] = (realPart[i] * realPart[i]) + (imagPart[i] * imagPart[i]);
5066
} // for i
5067
5068
} // void VB_Vector::getPS(VB_Vector &result) const
5069
5070
/*********************************************************************
5071
* Instance method to compute the power spectrum of this instance of *
5072
* VB_Vector. This method is derived from VoxBo_Fourier.pro's *
5073
* ReturnPS function. The power spectrum will be stored in result. *
5074
*********************************************************************/
5075
void VB_Vector::getPS(VB_Vector *result) const
5076
{
5077
this->getPS(*result);
5078
} // void VB_Vector::getPS(VB_Vector *result) const
5079
5080
/*********************************************************************
5081
* Instance method to compute the power spectrum of this instance of *
5082
* VB_Vector. This method is derived from VoxBo_Fourier.pro's *
5083
* ReturnPS function. *
5084
*********************************************************************/
5085
void VB_Vector::getPS() throw()
5086
{
5087
VB_Vector tempVec(this->getLength());
5088
this->getPS(tempVec);
5089
(*this) = tempVec;
5090
5091
} // void VB_Vector::getPS() throw()
5092
5093
/*********************************************************************
5094
* This method simply reverses the vector entries. *
5095
*********************************************************************/
5096
void VB_Vector::reverse()
5097
{
5098
5099
/*********************************************************************
5100
* Calling gsl_vector_reverse() to carry out the reversing of the *
5101
* vector elements. *
5102
*********************************************************************/
5103
gsl_vector_reverse(this->theVector);
5104
5105
} // void VB_Vector::reverse()
5106
5107
/*********************************************************************
5108
* This static method carries out Sinc Interpolation on the input *
5109
* VB_Vector object. This method is derived from the IDL function *
5110
* SincInterpo, found in VoxBo_Fourier.pro. *
5111
* *
5112
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5113
* ---------------- ----- ------------ *
5114
* timeSeries const VB_Vector& Input VB_Vector object that will*
5115
* be sinc interpolated (this *
5116
* VB_Vector represents a time *
5117
* series). *
5118
* expFactor int The number of times the original*
5119
* interval is to be expanded. For *
5120
* example, if the interval is 1 *
5121
* second long in the time series, *
5122
* but half second intervals are *
5123
* desired, then expFactor should *
5124
* be 2. *
5125
* newSignal VB_Vector& Reference to the VB_Vector *
5126
* object that will hold the result*
5127
* of the sinc interpolation. The *
5128
* length will be *
5129
* (timeseries.getLength() * *
5130
* expFactor). *
5131
* *
5132
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5133
* ----------------- ----- ------------ *
5134
* N/A *
5135
********************************************************************/
5136
void VB_Vector::sincInterpolation(const VB_Vector& timeSeries,
5137
const unsigned int expFactor, VB_Vector& newSignal)
5138
{
5139
5140
/*********************************************************************
5141
* If the length of timeSeries is less than 2 (mostly like 1, but *
5142
* could be zero), then an error message is printed and then this *
5143
* program exits. *
5144
*********************************************************************/
5145
if (timeSeries.getLength() < 2)
5146
{
5147
ostringstream errorMsg;
5148
errorMsg << "[" << __FUNCTION__ << "]: Need length to be >= 2. VB_Vector length = ["
5149
<< timeSeries.getLength() << "].";
5150
printErrorMsgAndExit(VB_ERROR, errorMsg.str(), 1);
5151
} // if
5152
5153
/*********************************************************************
5154
* length is used to store the length of timeSeries. *
5155
*********************************************************************/
5156
const unsigned long length = timeSeries.getLength();
5157
5158
/*********************************************************************
5159
* Now ensuring that newSignal is properly sized. *
5160
*********************************************************************/
5161
if (newSignal.getLength() != (length * expFactor))
5162
{
5163
newSignal.resize(length * expFactor);
5164
} // if
5165
5166
/*********************************************************************
5167
* We now compute the FFT of timeSeries. realPartFFT is used to hold *
5168
* the real part of timeSeries' FFT. Similarly, imagPartFFT is used *
5169
* to hold the imaginary part of timeSeries' FFT. *
5170
*********************************************************************/
5171
VB_Vector realPartFFT(timeSeries.getLength());
5172
VB_Vector imagPartFFT(timeSeries.getLength());
5173
timeSeries.fft(realPartFFT, imagPartFFT);
5174
5175
/*********************************************************************
5176
* halfLength is the floor of (timeSeries.getLength()/2). *
5177
*********************************************************************/
5178
int halfLength = timeSeries.getLength() / 2;
5179
5180
/*********************************************************************
5181
* phi[] is an array needed to carry out the sinc interpolation. *
5182
* phi[0] will always be zero. The "lower half" values of phi[] will *
5183
* be calculated and the "upper half" values of phi[] depend on its *
5184
* "lower half" values. Specifically, the "lower half" values are *
5185
* computed up to and including phi[halfLength/2]. Then the "lower *
5186
* half" values are reflected about the x-axis, and then about the *
5187
* vertical line y = halfLength. Assume that the x-axis indexes phi[] *
5188
* and the y-axis is phi[i]. When length is 10, the graph of phi[] *
5189
* looks like: *
5190
* *
5191
* *
5192
* | *
5193
* | *
5194
* | *
5195
* | *
5196
* | *
5197
* 1.570 | o *
5198
* | *
5199
* | *
5200
* | *
5201
* | *
5202
* 1.256 | o *
5203
* | *
5204
* | *
5205
* | *
5206
* .942 | o *
5207
* | *
5208
* | *
5209
* .628 | o *
5210
* | *
5211
* .314 | o *
5212
* | *
5213
* -----o---+---+---+---+---+---+---+---+---+---+ *
5214
* | 1 2 3 4 5 6 7 8 9 *
5215
* -.314 | o *
5216
* | *
5217
* -.628 | o *
5218
* | *
5219
* | *
5220
* -.942 | o *
5221
* | *
5222
* | *
5223
* | *
5224
* -1.256 | o *
5225
* | *
5226
* | *
5227
* | *
5228
* *
5229
*********************************************************************/
5230
double phi[length];
5231
memset(phi, 0, sizeof(double) * length);
5232
5233
/*********************************************************************
5234
* The following for loop is used to carry out the remaining work for *
5235
* the sinc interpolation. *
5236
*********************************************************************/
5237
for (unsigned int i = 0; i < expFactor; i++)
5238
{
5239
5240
/*********************************************************************
5241
* Calculating timeShift. *
5242
*********************************************************************/
5243
double timeShift = (double ) i / (double ) expFactor;
5244
5245
/*********************************************************************
5246
* Even number of values. *
5247
*********************************************************************/
5248
if ( (length % 2) == 0 )
5249
{
5250
5251
/*********************************************************************
5252
* Now calculating the values for phi[] (for the case of an even *
5253
* number of values in timeSeries[]. *
5254
*********************************************************************/
5255
for (int j = 1; j <= halfLength; j++)
5256
{
5257
5258
/*********************************************************************
5259
* The following line calculates the "lower half" values of phi[]. *
5260
*********************************************************************/
5261
phi[j] = (timeShift * TWOPI) / ((double ) length / (double ) j);
5262
5263
/*********************************************************************
5264
* When j is not equal to halfLength, then calculate the "upper half" *
5265
* values of phi[]. *
5266
*********************************************************************/
5267
if (j != halfLength)
5268
{
5269
phi[length - j] = -1.0 * phi[j];
5270
} // if
5271
5272
} // for j
5273
5274
} // if
5275
5276
/*********************************************************************
5277
* Odd number of values. *
5278
*********************************************************************/
5279
else
5280
{
5281
5282
/*********************************************************************
5283
* Now calculating the values for phi[] (for the case of an odd *
5284
* number of values in timeSeries[]. *
5285
*********************************************************************/
5286
for (int j = 1; j <= halfLength; j++)
5287
{
5288
5289
/*********************************************************************
5290
* The following line calculates the "lower half" values of phi[]. *
5291
*********************************************************************/
5292
phi[j] = (timeShift * TWOPI) / ((double ) length / (double ) j);
5293
5294
/*********************************************************************
5295
* The following line calculates the "upper half" values of phi[]. *
5296
*********************************************************************/
5297
phi[length - j] = -1.0 * phi[j];
5298
5299
} // for j
5300
5301
} // else
5302
5303
/*********************************************************************
5304
* The elements of shifterReal will hold the real part of the complex *
5305
* product [Complex(cos(phi[j]), sin(phi[j]) * *
5306
* Complex(realPartFFT[j], imagPartFFT[j])], while the elements of *
5307
* shifterImag hold the imaginary part. The following for loop *
5308
* populates both shifterReal and shifterImag. Eventually, we'll need *
5309
* to take the inverse FFT of (shifterReal + i*shifterImag). *
5310
*********************************************************************/
5311
VB_Vector shifterReal(length);
5312
VB_Vector shifterImag(length);
5313
for (unsigned int j = 0; j < length; j++)
5314
{
5315
5316
/*********************************************************************
5317
* Assigning the value of the jth element of shifterReal and the jth *
5318
* element of shifterImag. *
5319
*********************************************************************/
5320
shifterReal[j] = ((cos(phi[j]) * realPartFFT[j]) -
5321
(sin(phi[j]) * imagPartFFT[j]));
5322
shifterImag[j] = ((cos(phi[j]) * imagPartFFT[j]) +
5323
(sin(phi[j]) * realPartFFT[j]));
5324
5325
} // for j
5326
5327
/*********************************************************************
5328
* We must view shifterReal and shifterImag as a complex vector, *
5329
* whose jth element is (shifterReal[j] + i*shifterImag[j]). Since we *
5330
* do not have an acutal complex array, the inverse FFT of shifterReal*
5331
* and shifterImag need to be computed separately. To that end, *
5332
* realPartFFTShifter is used to hold the real part of the inverse *
5333
* FFT of shifterReal, imagPartFFTShifter is used to hold the *
5334
* imaginary part of of the inverse FFT of shifterReal; *
5335
* realPartFFTShifterImag is used to hold the real part of the *
5336
* inverse FFT of shifterImag and imagPartFFTShifter is resued to *
5337
* hold the imaginary part of of the inverse FFT of shifterImag (we *
5338
* are not interested in the imaginary part of the inverse FFT's of *
5339
* shifterReal). Ultimately, we will need just the real part of the *
5340
* inverse FFT of the complex vector (shifterReal + i*shifterImag). *
5341
*********************************************************************/
5342
VB_Vector realPartFFTShifter(shifterReal.getLength());
5343
VB_Vector imagPartFFTShifter(shifterReal.getLength());
5344
VB_Vector realPartFFTShifterImag(shifterImag.getLength());
5345
5346
/*********************************************************************
5347
* Now computing the inverse FFT's. *
5348
*********************************************************************/
5349
shifterReal.ifft(realPartFFTShifter, imagPartFFTShifter);
5350
shifterImag.ifft(realPartFFTShifterImag, imagPartFFTShifter);
5351
5352
/*********************************************************************
5353
* Let (a + i*b) be the complex vector representing the inverse FFT *
5354
* of shifterReal. Similarly, let (c + i*d) be the complex vector *
5355
* representing the inverse FFT of shifterImag. Then the real part of *
5356
* the the inverse FFT of the complex vector (shifterReal + *
5357
* i*shifterImag) is represented by the vector (a - d). So now we *
5358
* place the real part in realPartFFTShifter. *
5359
*********************************************************************/
5360
realPartFFTShifter -= imagPartFFTShifter;
5361
5362
/*********************************************************************
5363
* Now the real part of the inverse FFT of the complex vector *
5364
* (shifterReal + i*shifterImag) is assigned to newSignal *
5365
* appropriately. *
5366
*********************************************************************/
5367
for (unsigned int j = 0; j < length; j++)
5368
{
5369
newSignal[j * expFactor + i] = realPartFFTShifter[j];
5370
} // for j
5371
5372
} // for i
5373
5374
} // void VB_Vector::sincInterpolation(const VB_Vector& timeSeries,
5375
// const unsigned int expFactor, VB_Vector& newSignal)
5376
5377
/*********************************************************************
5378
* This constant instance method carries out Sinc Interpolation on *
5379
* this instance of VB_Vector. This method is derived from the IDL *
5380
* function SincInterpo, found in VoxBo_Fourier.pro. *
5381
* *
5382
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5383
* ---------------- ----- ------------ *
5384
* expFactor int The number of times the original*
5385
* interval is to be expanded. For *
5386
* example, if the interval is 1 *
5387
* second long in the time series, *
5388
* but half second intervals are *
5389
* desired, then expFactor should *
5390
* be 2. *
5391
* newSignal VB_Vector& Reference to the VB_Vector *
5392
* object that will hold the result*
5393
* of the sinc interpolation. The *
5394
* length will be *
5395
* (this->getLength() * expFactor).*
5396
* *
5397
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5398
* ----------------- ----- ------------ *
5399
* N/A *
5400
********************************************************************/
5401
void VB_Vector::sincInterpolation(const unsigned int expFactor,
5402
VB_Vector& newSignal) const
5403
{
5404
5405
/*********************************************************************
5406
* Simply calling the static method to compute the sinc interpolation.*
5407
*********************************************************************/
5408
VB_Vector::sincInterpolation(*this, expFactor, newSignal);
5409
5410
} // void VB_Vector::sincInterpolation(const unsigned int expFactor,
5411
// VB_Vector& newSignal) const
5412
5413
/*********************************************************************
5414
* This non-constant instance method carries out Sinc Interpolation on*
5415
* this instance of VB_Vector. This method is derived from the IDL *
5416
* function SincInterpo, found in VoxBo_Fourier.pro. *
5417
* *
5418
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5419
* ---------------- ----- ------------ *
5420
* expFactor int The number of times the original*
5421
* interval is to be expanded. For *
5422
* example, if the interval is 1 *
5423
* second long in the time series, *
5424
* but half second intervals are *
5425
* desired, then expFactor should *
5426
* be 2. *
5427
* *
5428
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5429
* ----------------- ----- ------------ *
5430
* N/A *
5431
********************************************************************/
5432
void VB_Vector::sincInterpolation(const unsigned int expFactor)
5433
{
5434
5435
/*********************************************************************
5436
* tempVector is a copy of this instance of VB_Vector. newSignal will *
5437
* be used to hold the interpolated signal. *
5438
*********************************************************************/
5439
VB_Vector tempVector(this);
5440
VB_Vector newSignal = VB_Vector();
5441
5442
/*********************************************************************
5443
* Doing the sinc interpolation and then assigning newSignal to this *
5444
* instance of VB_Vector. *
5445
*********************************************************************/
5446
tempVector.sincInterpolation(expFactor, newSignal);
5447
(*this) = newSignal;
5448
5449
} // void VB_Vector::sincInterpolation(const unsigned int expFactor)
5450
5451
/*********************************************************************
5452
* This method is a port of the IDL function PhaseShift, found in *
5453
* VoxBo_Fourier.pro. This method "shifts" a signal in time to *
5454
* provide an output vector that represents the same continuous *
5455
* signal sampled starting either later or earlier. This is *
5456
* accomplished by a simple shift of the phase of the sines that *
5457
* makes up the signal. Recall that a Fourier transform allows for a *
5458
* representation of any signal as the linear combination of *
5459
* sinusoids of different frequencies and phases. Effectively, we will*
5460
* add a constant to the phase of every frequency, shifting the data *
5461
* in time. *
5462
* *
5463
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5464
* ---------------- ----- ------------ *
5465
* tSeries const VB_Vector& The time series to be shifted.*
5466
* *
5467
* timeShift double Number of images to shift to *
5468
* the right, i.e., to lag the *
5469
* signal with respect to the *
5470
* original domain. *
5471
* *
5472
* shiftedSignal VB_Vector& The VB_Vector object that *
5473
* will hold the shifted signal. *
5474
* *
5475
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5476
* ----------------- ----- ------------ *
5477
* shiftedSignal VB_Vector& The shifted signal will be *
5478
* stored here. *
5479
********************************************************************/
5480
void VB_Vector::phaseShift(const VB_Vector& tSeries, const double timeShift, VB_Vector& shiftedSignal)
5481
{
5482
5483
/*********************************************************************
5484
* Making sure that shiftedSignal is appropriately sized. *
5485
*********************************************************************/
5486
if (tSeries.getLength() != shiftedSignal.getLength())
5487
{
5488
shiftedSignal.resize(tSeries.theVector->size);
5489
} // if
5490
5491
/*********************************************************************
5492
* phi[] is an array needed to carry out the phase shift. Now *
5493
* declaring phi[] and setting all of its elements to zero. *
5494
*********************************************************************/
5495
double *phi = new double[tSeries.theVector->size];
5496
memset(phi, 0, sizeof(double) * tSeries.getLength());
5497
5498
/*********************************************************************
5499
* Now calling makePhi() to fill in the values in phi[]. *
5500
*********************************************************************/
5501
VB_Vector::makePhi(phi, tSeries.getLength(), timeShift);
5502
5503
/*********************************************************************
5504
* Now taking the FFT of tSeries[]. realPartFFT is used to hold the *
5505
* real part of the FFT, while imagPartFFT is used to hold the *
5506
* imaginary part. *
5507
*********************************************************************/
5508
VB_Vector realPartFFT = VB_Vector();
5509
VB_Vector imagPartFFT = VB_Vector();
5510
tSeries.fft(realPartFFT, imagPartFFT);
5511
5512
/*********************************************************************
5513
* shifter is the filter by which the signal will be convolved to *
5514
* introduce the phase shift. It is constructed explicitly in the *
5515
* Fourier, i.e., frequency, domain. In the time domain, it may be *
5516
* described as an impulse (delta function) that has been shifted in *
5517
* time by the amount specified by timeShift. shifter is actually *
5518
* instantiated as two VB_Vector objects: shifterReal and shifterImag.*
5519
*********************************************************************/
5520
VB_Vector shifterReal(tSeries.getLength());
5521
VB_Vector shifterImag(tSeries.getLength());
5522
5523
/*********************************************************************
5524
* The following for loop carries out the convolution. *
5525
*********************************************************************/
5526
for (unsigned int i = 0; i < tSeries.getLength(); i++)
5527
{
5528
shifterReal[i] = (cos(phi[i]) * realPartFFT[i]) - (sin(phi[i]) * imagPartFFT[i]);
5529
shifterImag[i] = (cos(phi[i]) * imagPartFFT[i]) + (sin(phi[i]) * realPartFFT[i]);
5530
} // for i
5531
5532
/*********************************************************************
5533
* Now deleting the memory allocated for phi[]. *
5534
*********************************************************************/
5535
deleteArrMem(phi);
5536
5537
/*********************************************************************
5538
* Now we need to get the inverse FFT of shifter. Therefore, we need *
5539
* to take the inverse FFT's of shifterReal and shifterImag. Each of *
5540
* these inverse FFT's will have a real part and an imaginary part. *
5541
* *
5542
* For shifterReal, its inverse FFT parts will be stored in *
5543
* shifterRealIFFTReal and shifterRealIFFTImag. For shifterImag, its *
5544
* inverse FFT parts will be stored in shifterImagIFFTReal and *
5545
* shifterImagIFFTImag. *
5546
*********************************************************************/
5547
VB_Vector shifterRealIFFTReal(tSeries.getLength());
5548
VB_Vector shifterRealIFFTImag(tSeries.getLength());
5549
VB_Vector shifterImagIFFTReal(tSeries.getLength());
5550
VB_Vector shifterImagIFFTImag(tSeries.getLength());
5551
5552
/*********************************************************************
5553
* Now computing the inverse FFT's of shifterReal and shifterImag. *
5554
*********************************************************************/
5555
shifterReal.ifft(shifterRealIFFTReal, shifterRealIFFTImag);
5556
shifterImag.ifft(shifterImagIFFTReal, shifterImagIFFTImag);
5557
5558
/*********************************************************************
5559
* Now assigning the real part of the inverse FFT of shifter to *
5560
* shiftedSignal. *
5561
*********************************************************************/
5562
shiftedSignal = shifterRealIFFTReal - shifterImagIFFTImag;
5563
5564
} // void VB_Vector::phaseShift(const VB_Vector& tSeries,
5565
// const double timeShift, VB_Vector& shiftedSignal)
5566
5567
/*********************************************************************
5568
* This instance method calls the static method *
5569
* VB_Vector::phaseShift() to compute the phase shift, using this *
5570
* instance of VB_Vector as the first argument to the static *
5571
* VB_Vector::phaseShift() method. *
5572
*********************************************************************/
5573
void VB_Vector::phaseShift(const double timeShift,
5574
VB_Vector& shiftedSignal) const
5575
{
5576
VB_Vector::phaseShift(*this, timeShift, shiftedSignal);
5577
} // void VB_Vector::phaseShift(const double timeShift,
5578
// VB_Vector& shiftedSignal) const
5579
5580
/*********************************************************************
5581
* This instance method calls the static method *
5582
* VB_Vector::phaseShift() to compute the phase shift, using this *
5583
* instance of VB_Vector as the first argument to the static *
5584
* VB_Vector::phaseShift() method. Then the resulting phase shifted *
5585
* VB_Vector is assigned to this instance of VB_Vector. *
5586
*********************************************************************/
5587
void VB_Vector::phaseShift(const double timeShift)
5588
{
5589
5590
/*********************************************************************
5591
* tempVector will be used to hold the phase shifted vector. *
5592
*********************************************************************/
5593
VB_Vector tempVector = VB_Vector();
5594
5595
/*********************************************************************
5596
* Calling the static method VB_Vector::phaseShift() to carry out the *
5597
* phase shifting. Then tempVector is assigned to this instance of *
5598
* VB_Vector. *
5599
*********************************************************************/
5600
VB_Vector::phaseShift(*this, timeShift, tempVector);
5601
(*this) = tempVector;
5602
5603
} // void VB_Vector::phaseShift(const double timeShift)
5604
5605
/*********************************************************************
5606
* phi[] is an array needed to carry out the phase shift. phi[0] *
5607
* will always be zero. The "lower half" values of phi[] will be *
5608
* calculated and the "upper half" values of phi[] depend on its *
5609
* "lower half" values. Specifically, the "lower half" values are *
5610
* computed up to and including phi[length/2]. Then the "lower half" *
5611
* values are reflected about the x-axis, and then about the *
5612
* vertical line y = (length/2). Assume that the x-axis indexes phi[] *
5613
* and the y-axis is phi[i]. If length is 10, the graph of phi[] looks*
5614
* like: *
5615
* *
5616
* *
5617
* | *
5618
* | *
5619
* | *
5620
* | *
5621
* | *
5622
* 1.570 | o *
5623
* | *
5624
* | *
5625
* | *
5626
* | *
5627
* 1.256 | o *
5628
* | *
5629
* | *
5630
* | *
5631
* .942 | o *
5632
* | *
5633
* | *
5634
* .628 | o *
5635
* | *
5636
* .314 | o *
5637
* | *
5638
* -----o---+---+---+---+---+---+---+---+---+---+ *
5639
* | 1 2 3 4 5 6 7 8 9 *
5640
* -.314 | o *
5641
* | *
5642
* -.628 | o *
5643
* | *
5644
* | *
5645
* -.942 | o *
5646
* | *
5647
* | *
5648
* | *
5649
* -1.256 | o *
5650
* | *
5651
* | *
5652
* | *
5653
* *
5654
* *
5655
* *
5656
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5657
* ---------------- ----- ------------ *
5658
* phi double * Array used to hold the result. *
5659
* *
5660
* length int Length of phi[]. *
5661
* *
5662
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5663
* ----------------- ----- ------------ *
5664
* phi double * Holds the results. *
5665
********************************************************************/
5666
void VB_Vector::makePhi(double *phi, const int length, const double timeShift)
5667
{
5668
5669
/*********************************************************************
5670
* Setting all elements in phi[] to zero. *
5671
*********************************************************************/
5672
memset(phi, 0, sizeof(double) * length);
5673
5674
/*********************************************************************
5675
* halfLength holds the floor of length / 2. *
5676
*********************************************************************/
5677
int halfLength = length / 2;
5678
5679
/*********************************************************************
5680
* Even number of values. *
5681
*********************************************************************/
5682
if ( (length % 2) == 0 )
5683
{
5684
5685
/*********************************************************************
5686
* Now calculating the values for phi[] (for the case of an even *
5687
* number of values in timeSeries[]. *
5688
*********************************************************************/
5689
for (int j = 1; j <= halfLength; j++)
5690
{
5691
5692
/*********************************************************************
5693
* The following line calculates the "lower half" values of phi[]. *
5694
*********************************************************************/
5695
phi[j] = -1.0 * (timeShift * TWOPI) / ((double ) length / (double ) j);
5696
5697
/*********************************************************************
5698
* When j is not equal to halfLength, then calculate the "upper half" *
5699
* values of phi[]. *
5700
*********************************************************************/
5701
if (j != halfLength)
5702
{
5703
phi[length - j] = -1.0 * phi[j];
5704
} // if
5705
5706
} // for j
5707
5708
} // if
5709
5710
/*********************************************************************
5711
* Odd number of values. *
5712
*********************************************************************/
5713
else
5714
{
5715
5716
/*********************************************************************
5717
* Now calculating the values for phi[] (for the case of an odd *
5718
* number of values in timeSeries[]. *
5719
*********************************************************************/
5720
for (int j = 1; j <= halfLength; j++)
5721
{
5722
5723
/*********************************************************************
5724
* The following line calculates the "lower half" values of phi[]. *
5725
*********************************************************************/
5726
phi[j] = -1.0 * (timeShift * TWOPI) / ((double ) length / (double ) j);
5727
5728
/*********************************************************************
5729
* The following line calculates the "upper half" values of phi[]. *
5730
*********************************************************************/
5731
phi[length - j] = -1.0 * phi[j];
5732
5733
} // for j
5734
5735
} // else
5736
5737
} // void VB_Vector::makePhi(double *phi, const int length, const double timeShift)
5738
5739
/*********************************************************************
5740
* This private static method computes the sine of a complex number *
5741
* that is in Cartesian form. The real part of the sine is assigned *
5742
* to real. The imaginary part of the sine is assigned to imaginary. *
5743
* *
5744
* Here' how the sine is computed: *
5745
* *
5746
* Let z = x + iy, a complex number. Then *
5747
* *
5748
* sin(z) = (sin(x)) * ((e^y + e^(-y)))/2 + *
5749
* i * (cos(x)) * ((e^y - e^(-y)))/2 *
5750
* *
5751
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5752
* ---------------- ----- ------------ *
5753
* x const double The real part of the complex *
5754
* number whose sine is to be *
5755
* computed. *
5756
* y const double The imaghinary part of the *
5757
* complex number whose sine is to *
5758
* be computed. *
5759
* real double& Will hold the real part of the *
5760
* sine. *
5761
* imaginary double& Will hold the imaginary part of *
5762
* the sine. *
5763
* *
5764
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5765
* ----------------- ----- ------------ *
5766
* None. *
5767
* *
5768
* EXCEPTIONS THROWN: *
5769
* ------------------ *
5770
* None. *
5771
*********************************************************************/
5772
void VB_Vector::complexSin(const double x, const double y, double& real, double& imaginary)
5773
{
5774
5775
/*********************************************************************
5776
* Now getting the needed exponentials. *
5777
*********************************************************************/
5778
double eToTheY = exp(y);
5779
double eToTheMinusY = exp(-y);
5780
5781
/*********************************************************************
5782
* Now computing and assigning the real and imaginary parts of the *
5783
* sine function. *
5784
*********************************************************************/
5785
real = sin(x) * ((eToTheY + eToTheMinusY) / 2.0);
5786
imaginary = cos(x) * ((eToTheY - eToTheMinusY) / 2.0);
5787
5788
} // void VB_Vector::complexSin(const double x, const double y, double& real, double& imaginary)
5789
5790
/*********************************************************************
5791
* This private static method computes the cosine of a complex number *
5792
* that is in Cartesian form. The real part of the cosine is assigned *
5793
* to real. The imaginary part of the cosine is assigned to imaginary.*
5794
* *
5795
* Here' how the cosine is computed: *
5796
* *
5797
* Let z = x + iy, a complex number. Then *
5798
* *
5799
* cos(z) = sin(z + Pi/2) *
5800
* *
5801
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5802
* ---------------- ----- ------------ *
5803
* x const double The real part of the complex *
5804
* number whose cosine is to be *
5805
* computed. *
5806
* y const double The imaghinary part of the *
5807
* complex number whose cosine is *
5808
* to be computed. *
5809
* real double& Will hold the real part of the *
5810
* cosine. *
5811
* imaginary double& Will hold the imaginary part of *
5812
* the cosine. *
5813
* *
5814
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5815
* ----------------- ----- ------------ *
5816
* None. *
5817
* *
5818
* EXCEPTIONS THROWN: *
5819
* ------------------ *
5820
* None. *
5821
*********************************************************************/
5822
void VB_Vector::complexCos(const double x, const double y, double& real, double& imaginary)
5823
{
5824
5825
/*********************************************************************
5826
* Now computing the cosine. *
5827
*********************************************************************/
5828
VB_Vector::complexSin(x + (PI / 2.0), y, real, imaginary);
5829
5830
} // void VB_Vector::complexCos(const double x, const double y, double& real, double& imaginary)
5831
5832
5833
/*********************************************************************
5834
* This instance method normalizes the magnitude component of a signal*
5835
* to unity while preserving the phase component. This method is a *
5836
* port of the IDL NormMag function, found in VoxBo_Fourier.pro. *
5837
* *
5838
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5839
* ---------------- ----- ------------ *
5840
* this VB_Vector * The input signal. *
5841
* *
5842
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5843
* ----------------- ----- ------------ *
5844
* None. *
5845
* *
5846
* EXCEPTIONS THROWN: *
5847
* ------------------ *
5848
* None. *
5849
*********************************************************************/
5850
void VB_Vector::normMag()
5851
{
5852
5853
/*********************************************************************
5854
* We now need to compute the FFT of this instance of VB_Vector. *
5855
* realFFT will hold the real part of the FFT, while imagIFFT will *
5856
* hold the imaginary part of the FFT. *
5857
*********************************************************************/
5858
VB_Vector realFFT = VB_Vector();
5859
VB_Vector imagIFFT = VB_Vector();
5860
this->fft(realFFT, imagIFFT);
5861
5862
/*********************************************************************
5863
* magImage[] will hold the modulus of each complex number from the *
5864
* FFT of this instance of VB_Vector. *
5865
*********************************************************************/
5866
double *magImage = new double[this->getLength()];
5867
5868
/*********************************************************************
5869
* theZeros is used to hold the indices of where magImage is zero. *
5870
*********************************************************************/
5871
vector<unsigned long> theZeros;
5872
5873
/*********************************************************************
5874
* The following for loop populates magImage. *
5875
*********************************************************************/
5876
for (unsigned long i = 0; i < this->getLength(); i++)
5877
{
5878
magImage[i] = sqrt((realFFT[i] * realFFT[i]) + (imagIFFT[i] * imagIFFT[i]));
5879
5880
/*********************************************************************
5881
* NOTE: If the modulus is zero, then the phase angle is undefined. *
5882
* To get around this, magImage[i] is set to 1.0, temporarily. Also, *
5883
* the index is added to theZeros. *
5884
*********************************************************************/
5885
if (magImage[i] == 0.0)
5886
{
5887
magImage[i] = 1.0;
5888
theZeros.push_back(i);
5889
} // if
5890
5891
} // for i
5892
5893
/*********************************************************************
5894
* phaseImage[i] will be used to hold the arc-cosine of realFTT[i] *
5895
* divided by magImage[i]. *
5896
*********************************************************************/
5897
double *phaseImage = new double[this->getLength()];
5898
5899
/*********************************************************************
5900
* This for loop populates realFFT[i]. *
5901
*********************************************************************/
5902
for (unsigned long i = 0; i < this->getLength(); i++)
5903
{
5904
phaseImage[i] = acos(realFFT[i] / magImage[i]);
5905
5906
/*********************************************************************
5907
* Since the domain of acos() is limited, we need to check the sign *
5908
* of imagIFFT[i]. *
5909
*********************************************************************/
5910
if (imagIFFT[i] < 0.0)
5911
{
5912
phaseImage[i] = TWOPI - phaseImage[i];
5913
} // if
5914
} // for i
5915
5916
/*********************************************************************
5917
* The following for loop is sued to "restore" the zero values in *
5918
* magImage[] and phaseImage[]. *
5919
*********************************************************************/
5920
for (unsigned long i = 0; i < theZeros.size(); i++)
5921
{
5922
magImage[i] = 0.0;
5923
phaseImage[i] = 0.0;
5924
} // for i
5925
5926
/*********************************************************************
5927
* magImageMax is assigned the maximum value of magImage[]. *
5928
*********************************************************************/
5929
double magImageMax = getMax(magImage, this->getLength());
5930
5931
/*********************************************************************
5932
* realPart is used to hold the real part of the normalized signal *
5933
* and imagPart is used to hold the imaginary part of the normalized *
5934
* signal. *
5935
*********************************************************************/
5936
VB_Vector realPart(this->getLength());
5937
VB_Vector imagPart(this->getLength());
5938
5939
/*********************************************************************
5940
* The following for loop is used to populate both realPart and *
5941
* imagPart. *
5942
*********************************************************************/
5943
for (unsigned long i =0; i < this->getLength(); i++)
5944
{
5945
realPart[i] = (magImage[i] / magImageMax) * cos(phaseImage[i]);
5946
imagPart[i] = (magImage[i] / magImageMax) * sin(phaseImage[i]);
5947
} // for i
5948
5949
/*********************************************************************
5950
* We now need to convert back to the time domain. realPartFFTReal *
5951
* will be used to hold the real part of the inverse FFT of realPart *
5952
* and realPartFFTImag will be used to hold the imaginary part of the *
5953
* inverse FFT of realPart. *
5954
*********************************************************************/
5955
VB_Vector realPartFFTReal = VB_Vector();
5956
VB_Vector realPartFFTImag = VB_Vector();
5957
realPart.ifft(realPartFFTReal, realPartFFTImag);
5958
5959
/*********************************************************************
5960
* imagPartFFTReal will be used to hold the real part of the inverse *
5961
* FFT of imagPart and imagPartFFTImag will be used to hold the *
5962
* imaginary part of the inverse FFT of imagPart. *
5963
*********************************************************************/
5964
VB_Vector imagPartFFTReal = VB_Vector();
5965
VB_Vector imagPartFFTImag = VB_Vector();
5966
imagPart.ifft(imagPartFFTReal, imagPartFFTImag);
5967
5968
/*********************************************************************
5969
* Now setting this instance of VB_Vector to its normalized form. *
5970
*********************************************************************/
5971
(*this) = realPartFFTReal - imagPartFFTImag;
5972
5973
/*********************************************************************
5974
* Deleting previously allocated memory. *
5975
*********************************************************************/
5976
deleteArrMem(phaseImage);
5977
deleteArrMem(magImage);
5978
5979
} // void VB_Vector::normMag()
5980
5981
/*********************************************************************
5982
* This instance method normalizes the magnitude component of a signal*
5983
* to unity while preserving the phase component. This method is a *
5984
* port of the IDL NormMag function, found in VoxBo_Fourier.pro. *
5985
* *
5986
* INPUT VARIABLES: TYPE: DESCRIPTION: *
5987
* ---------------- ----- ------------ *
5988
* normalizedVec VB_Vector& Will hold the normalized vector.*
5989
* *
5990
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
5991
* ----------------- ----- ------------ *
5992
* None. *
5993
* *
5994
* EXCEPTIONS THROWN: *
5995
* ------------------ *
5996
* None. *
5997
*********************************************************************/
5998
void VB_Vector::normMag(VB_Vector& normalizedVec) const
5999
{
6000
6001
/*********************************************************************
6002
* Setting normalizedVec equal to this instance of VB_Vector and then *
6003
* normalizing normalizedVec. *
6004
*********************************************************************/
6005
normalizedVec = (*this);
6006
normalizedVec.normMag();
6007
6008
} // void VB_Vector::normMag(VB_Vector& normalizedVec) const
6009
6010
/*********************************************************************
6011
* This method applies the input function (actually a function *
6012
* pointer) to each element of this instance of VB_Vector. The *
6013
* input function must take a single double argument and return a *
6014
* double. This method is meant to apply sqrt(), cos(), sin(), etc., *
6015
* to each element. *
6016
* *
6017
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6018
* ---------------- ----- ------------ *
6019
* theFunction double (*) (double) A pointer to a function *
6020
* that takes a single *
6021
* double argument and *
6022
* returns a double. *
6023
* *
6024
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6025
* ----------------- ----- ------------ *
6026
* None. *
6027
* *
6028
* EXCEPTIONS THROWN: *
6029
* ------------------ *
6030
* None. *
6031
*********************************************************************/
6032
void VB_Vector::applyFunction(double (*theFunction)(double)) throw()
6033
{
6034
6035
/*********************************************************************
6036
* The following for loop is used to traverse this instance of *
6037
* VB_Vector and apply the input function to each element. *
6038
*********************************************************************/
6039
for (unsigned long i = 0; i < this->getLength(); i++)
6040
{
6041
(*this)[i] = (*theFunction)((*this)[i]);
6042
} // for i
6043
6044
} // void VB_Vector::applyFunction(double (*theFunction)(double)) throw()
6045
6046
/*********************************************************************
6047
* This method applies the input function (actually a function *
6048
* pointer) to each element of this instance of VB_Vector and stores *
6049
* the result in theResult. The input function must take a single *
6050
* double argument and return a double. This method is meant to *
6051
* apply sqrt(), cos(), sin(), etc., to each element. *
6052
* *
6053
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6054
* ---------------- ----- ------------ *
6055
* theFunction double (*) (double) A pointer to a function *
6056
* that takes a single *
6057
* double argument and *
6058
* returns a double. *
6059
* theResult VB_Vector& Will hold the result. *
6060
* *
6061
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6062
* ----------------- ----- ------------ *
6063
* None. *
6064
* *
6065
* EXCEPTIONS THROWN: *
6066
* ------------------ *
6067
* None. *
6068
*********************************************************************/
6069
void VB_Vector::applyFunction(double (*theFunction)(double), VB_Vector& theResult) const throw()
6070
{
6071
6072
/*********************************************************************
6073
* Setting theResult equal to this instance of VB_Vector and then *
6074
* calling VB_Vector::applyFunction(). *
6075
*********************************************************************/
6076
theResult = (*this);
6077
theResult.applyFunction(theFunction);
6078
6079
} // void VB_Vector::applyFunction(double (*theFunction)(double), VB_Vector& theResult) const throw()
6080
6081
/*********************************************************************
6082
* This method applies the input function (actually a function *
6083
* pointer) to each element of this instance of VB_Vector and stores *
6084
* the result in theResult. The input function must take a single *
6085
* double argument and return a double. This method is meant to *
6086
* apply sqrt(), cos(), sin(), etc., to each element. *
6087
* *
6088
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6089
* ---------------- ----- ------------ *
6090
* theFunction double (*) (double) A pointer to a function *
6091
* that takes a single *
6092
* double argument and *
6093
* returns a double. *
6094
* theResult VB_Vector * Will hold the result. *
6095
* *
6096
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6097
* ----------------- ----- ------------ *
6098
* None. *
6099
* *
6100
* EXCEPTIONS THROWN: *
6101
* ------------------ *
6102
* None. *
6103
*********************************************************************/
6104
void VB_Vector::applyFunction(double (*theFunction)(double), VB_Vector *theResult) const throw()
6105
{
6106
6107
/*********************************************************************
6108
* Setting theResult equal to this instance of VB_Vector and then *
6109
* calling VB_Vector::applyFunction(). *
6110
*********************************************************************/
6111
(*theResult) = (*this);
6112
theResult->applyFunction(theFunction);
6113
6114
} // void VB_Vector::applyFunction(double (*theFunction)(double), VB_Vector *theResult) const throw()
6115
6116
/*********************************************************************
6117
* This method carries out an element-by-element multiplication *
6118
* between this instance of VB_Vector and the input VB_Vector. This *
6119
* instance of VB_Vector is overwritten with the new values. *
6120
* *
6121
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6122
* ---------------- ----- ------------ *
6123
* vec const VB_Vector * The input vector. *
6124
* *
6125
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6126
* ----------------- ----- ------------ *
6127
* N/A *
6128
* *
6129
* EXCEPTIONS THROWN: *
6130
* ------------------ *
6131
* None. *
6132
*********************************************************************/
6133
void VB_Vector::elementByElementMult(const VB_Vector *vec) throw()
6134
{
6135
6136
/*********************************************************************
6137
* The following try/catch blocks are used to ensure that the vector *
6138
* lengths are equal. *
6139
*********************************************************************/
6140
try
6141
{
6142
VB_Vector::checkVectorLengths(this->theVector, vec->theVector, __LINE__, __FILE__, __FUNCTION__);
6143
} // try
6144
catch (GenericExcep &e)
6145
{
6146
e.whatNoExit(__LINE__, __FILE__, __FUNCTION__);
6147
} // catch
6148
6149
/*********************************************************************
6150
* The following for loop carries out the element-by-element *
6151
* multiplication. *
6152
*********************************************************************/
6153
for (unsigned long i = 0; i < this->getLength(); i++)
6154
{
6155
(*this)[i] *= (*vec)[i];
6156
} // for i
6157
6158
} // void VB_Vector::elementByElementMult(const VB_Vector *vec) throw()
6159
6160
/*********************************************************************
6161
* This method carries out an element-by-element multiplication *
6162
* between this instance of VB_Vector and the input VB_Vector. This *
6163
* instance of VB_Vector is overwritten with the new values. *
6164
* *
6165
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6166
* ---------------- ----- ------------ *
6167
* vec const VB_Vector& The input vector. *
6168
* *
6169
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6170
* ----------------- ----- ------------ *
6171
* N/A *
6172
* *
6173
* EXCEPTIONS THROWN: *
6174
* ------------------ *
6175
* None. *
6176
*********************************************************************/
6177
void VB_Vector::elementByElementMult(const VB_Vector& vec) throw()
6178
{
6179
6180
/*********************************************************************
6181
* Calling VB_Vector::elementByElementMult(const VB_Vector *vec) to *
6182
* carry out the element-by-element multiplication. *
6183
*********************************************************************/
6184
this->elementByElementMult(&vec);
6185
6186
} // void VB_Vector::elementByElementMult(const VB_Vector& vec) throw()
6187
6188
/*********************************************************************
6189
* This static method multiplies 2 "complex" VB_Vectors: the first *
6190
* complex VB_Vector has its real and imaginary parts passed in; *
6191
* similarly for the second VB_Vector. The real part of the product *
6192
* will be saved to realProd and the imaginary part will be saved to *
6193
* imagProd. *
6194
* *
6195
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6196
* ---------------- ----- ------------ *
6197
* real1 const VB_Vector& The real part of the first *
6198
* VB_Vector. *
6199
* imag1 const VB_Vector& The imaginary part of the *
6200
* first VB_Vector. *
6201
* real2 const VB_Vector& The real part of the second *
6202
* VB_Vector. *
6203
* imag2 const VB_Vector& The imaginary part of the *
6204
* second VB_Vector. *
6205
* realProd VB_Vector& Will hold the real part of *
6206
* the product. *
6207
* imagProd VB_Vector& Will hold the imaginary part *
6208
* of the product. *
6209
* *
6210
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6211
* ----------------- ----- ------------ *
6212
* N/A *
6213
* *
6214
* EXCEPTIONS THROWN: *
6215
* ------------------ *
6216
* GenericExcep *
6217
*********************************************************************/
6218
void VB_Vector::compMult(const VB_Vector& real1, const VB_Vector& imag1,
6219
const VB_Vector& real2, const VB_Vector& imag2, VB_Vector& realProd,
6220
VB_Vector& imagProd) throw (GenericExcep)
6221
{
6222
6223
/*********************************************************************
6224
* First, we need to ensure that all the vector lengths are equal. If *
6225
* any of the sizes do not match, then a GenericExcep is thrown. *
6226
*********************************************************************/
6227
VB_Vector::checkVectorLengths(real1.theVector, imag1.theVector, __LINE__, __FILE__, __FUNCTION__);
6228
VB_Vector::checkVectorLengths(real2.theVector, imag2.theVector, __LINE__, __FILE__, __FUNCTION__);
6229
VB_Vector::checkVectorLengths(real1.theVector, imag2.theVector, __LINE__, __FILE__, __FUNCTION__);
6230
6231
/*********************************************************************
6232
* Ensuring that realProd and imagProd are appropriately sized. *
6233
*********************************************************************/
6234
if (real1.getLength() != realProd.getLength())
6235
{
6236
realProd.resize(real1.theVector->size);
6237
} // if
6238
if (real1.getLength() != imagProd.getLength())
6239
{
6240
imagProd.resize(real1.theVector->size);
6241
} // if
6242
6243
/*********************************************************************
6244
* The following for loop carries out the complex multiplication. The *
6245
* real part of the product is saved to realProd, while the imaginary *
6246
* part is saved to imagProd. *
6247
*********************************************************************/
6248
for (unsigned long i = 0; i < real1.theVector->size; i++)
6249
{
6250
realProd.theVector->data[i] = (real1.theVector->data[i] * real2.theVector->data[i]) - (imag1.theVector->data[i] * imag2.theVector->data[i]);
6251
imagProd.theVector->data[i] = (real1.theVector->data[i] * imag2.theVector->data[i]) + (imag1.theVector->data[i] * real2.theVector->data[i]);
6252
} // for i
6253
6254
} // void VB_Vector::compMult(const VB_Vector& real1, const VB_Vector& imag1,
6255
// const VB_Vector& real2, const VB_Vector& imag2, VB_Vector& realProd,
6256
// VB_Vector& imagProd) throw (GenericExcep)
6257
6258
/*********************************************************************
6259
* This static method computes the FFT of a "complex" VB_Vector. The *
6260
* real part of the VB_Vector is passed in as real and the imaginary *
6261
* part of the VB_Vector is passed in as imag. The real part of the *
6262
* FFT is saved in realFFT while the imaginary part is saved to *
6263
* imagIFFT. *
6264
* *
6265
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6266
* ---------------- ----- ------------ *
6267
* real const VB_Vector& The real part of the input *
6268
* VB_Vector. *
6269
* imag const VB_Vector& The imaginary part of the *
6270
* input VB_Vector. *
6271
* realFFT VB_Vector& Will hold the real part *
6272
* of the FFT. *
6273
* imagIFFT VB_Vector& Will hold the imaginary part *
6274
* of the FFT. *
6275
* *
6276
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6277
* ----------------- ----- ------------ *
6278
* N/A *
6279
* *
6280
* EXCEPTIONS THROWN: *
6281
* ------------------ *
6282
* GenericExcep *
6283
*********************************************************************/
6284
void VB_Vector::complexFFT(const VB_Vector& real,
6285
const VB_Vector& imag, VB_Vector& realFFT, VB_Vector& imagIFFT)
6286
throw (GenericExcep)
6287
{
6288
6289
/*********************************************************************
6290
* Ensuring that the input VB_Vectors real and imag have the same *
6291
* size. If they don't have the same size, then a GenericExcep is *
6292
* thrown. *
6293
*********************************************************************/
6294
VB_Vector::checkVectorLengths(real.theVector, imag.theVector, __LINE__, __FILE__, __FUNCTION__);
6295
6296
/*********************************************************************
6297
* Ensuring that realFFT and imagIFFT are sized appropriately. *
6298
*********************************************************************/
6299
if (real.getLength() != realFFT.getLength())
6300
{
6301
realFFT.resize(real.theVector->size);
6302
} // if
6303
if (real.getLength() != imagIFFT.getLength())
6304
{
6305
imagIFFT.resize(real.theVector->size);
6306
} // if
6307
6308
/*********************************************************************
6309
* realFFTreal and realFFTimag will hold the real and imaginary parts *
6310
* of the FFT of real, respectively. *
6311
*********************************************************************/
6312
VB_Vector realFFTreal(real.getLength());
6313
VB_Vector realFFTimag(real.getLength());
6314
6315
/*********************************************************************
6316
* imagIFFTreal and imagIFFTimag will hold the real and imaginary *
6317
* parts of the FFT of imag, respectively. *
6318
*********************************************************************/
6319
VB_Vector imagIFFTreal(real.getLength());
6320
VB_Vector imagIFFTimag(real.getLength());
6321
6322
/*********************************************************************
6323
* Getting the FFT of real. *
6324
*********************************************************************/
6325
real.fft(realFFTreal, realFFTimag);
6326
6327
/*********************************************************************
6328
* Getting the FFT of imag. *
6329
*********************************************************************/
6330
imag.fft(imagIFFTreal, imagIFFTimag);
6331
6332
/*********************************************************************
6333
* Now realFFT will hold the real part of the FFT of *
6334
* (real + i * imag). *
6335
*********************************************************************/
6336
realFFT = realFFTreal - imagIFFTimag;
6337
6338
/*********************************************************************
6339
* Now imagIFFT will hold the imaginary part of the FFT of *
6340
* (real + i * imag). *
6341
*********************************************************************/
6342
imagIFFT = realFFTimag + imagIFFTreal;
6343
6344
} // void VB_Vector::complexFFT(const VB_Vector& real,
6345
// const VB_Vector& imag, VB_Vector& realFFT, VB_Vector& imagIFFT)
6346
// throw (GenericExcep)
6347
6348
/*********************************************************************
6349
* This static method computes the inverse FFT of a "complex" *
6350
* VB_Vector. The real part of the VB_Vector is passed in as real and *
6351
* the imaginary part of the VB_Vector is passed in as imag. The *
6352
* real part of the inverse FFT is saved in realIFFT while the *
6353
* imaginary part is saved to imagIFFT. *
6354
* *
6355
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6356
* ---------------- ----- ------------ *
6357
* real const VB_Vector& The real part of the input *
6358
* VB_Vector. *
6359
* imag const VB_Vector& The imaginary part of the *
6360
* input VB_Vector. *
6361
* realIFFT VB_Vector& Will hold the real part *
6362
* of the inverse FFT. *
6363
* imagIFFT VB_Vector& Will hold the imaginary part *
6364
* of the inverse FFT. *
6365
* *
6366
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6367
* ----------------- ----- ------------ *
6368
* N/A *
6369
* *
6370
* EXCEPTIONS THROWN: *
6371
* ------------------ *
6372
* GenericExcep *
6373
*********************************************************************/
6374
void VB_Vector::complexIFFT(const VB_Vector& real,
6375
const VB_Vector& imag, VB_Vector& realIFFT, VB_Vector& imagIFFT)
6376
throw (GenericExcep)
6377
{
6378
6379
/*********************************************************************
6380
* Ensuring that the input VB_Vectors real and imag have the same *
6381
* size. If they don't have the same size, then a GenericExcep is *
6382
* thrown. *
6383
*********************************************************************/
6384
VB_Vector::checkVectorLengths(real.theVector, imag.theVector, __LINE__, __FILE__, __FUNCTION__);
6385
6386
/*********************************************************************
6387
* Ensuring that realIFFT and imagIFFT are sized appropriately. *
6388
*********************************************************************/
6389
if (real.getLength() != realIFFT.getLength())
6390
{
6391
realIFFT.resize(real.theVector->size);
6392
} // if
6393
if (real.getLength() != imagIFFT.getLength())
6394
{
6395
imagIFFT.resize(real.theVector->size);
6396
} // if
6397
6398
/*********************************************************************
6399
* realIFFTreal and realIFFTimag will hold the real and imaginary *
6400
* parts of the inverse FFT of real, respectively. *
6401
*********************************************************************/
6402
VB_Vector realIFFTreal(real.getLength());
6403
VB_Vector realIFFTimag(real.getLength());
6404
6405
/*********************************************************************
6406
* imagIFFTreal and imagIFFTimag will hold the real and imaginary *
6407
* parts of the inverse FFT of imag, respectively. *
6408
*********************************************************************/
6409
VB_Vector imagIFFTreal(real.getLength());
6410
VB_Vector imagIFFTimag(real.getLength());
6411
6412
/*********************************************************************
6413
* Getting the inverse FFT of real. *
6414
*********************************************************************/
6415
real.ifft(realIFFTreal, realIFFTimag);
6416
6417
/*********************************************************************
6418
* Getting the inverse FFT of imag. *
6419
*********************************************************************/
6420
imag.ifft(imagIFFTreal, imagIFFTimag);
6421
6422
/*********************************************************************
6423
* Now realIFFT will hold the real part of the inverse FFT of *
6424
* (real + i * imag). *
6425
*********************************************************************/
6426
realIFFT = realIFFTreal - imagIFFTimag;
6427
6428
/*********************************************************************
6429
* Now imagIFFT will hold the imaginary part of the inverse FFT of *
6430
* (real + i * imag). *
6431
*********************************************************************/
6432
imagIFFT = realIFFTimag + imagIFFTreal;
6433
6434
} // void VB_Vector::complexIFFT(const VB_Vector& real,
6435
// const VB_Vector& imag, VB_Vector& realIFFT, VB_Vector& imagIFFT)
6436
// throw (GenericExcep)
6437
6438
/*********************************************************************
6439
* This static method computes the inverse FFT of a "complex" *
6440
* VB_Vector. The real part of the VB_Vector is passed in as real and *
6441
* the imaginary part of the VB_Vector is passed in as imag. The *
6442
* real part of the inverse FFT is saved in realIFFT while the *
6443
* imaginary part is not saved. *
6444
* *
6445
* INPUT VARIABLES: TYPE: DESCRIPTION: *
6446
* ---------------- ----- ------------ *
6447
* real const VB_Vector& The real part of the input *
6448
* VB_Vector. *
6449
* imag const VB_Vector& The imaginary part of the *
6450
* input VB_Vector. *
6451
* realIFFT VB_Vector& Will hold the real part *
6452
* of the inverse FFT. *
6453
* *
6454
* OUTPUT VARIABLES: TYPE: DESCRIPTION: *
6455
* ----------------- ----- ------------ *
6456
* N/A *
6457
* *
6458
* EXCEPTIONS THROWN: *
6459
* ------------------ *
6460
* GenericExcep *
6461
*********************************************************************/
6462
void VB_Vector::complexIFFTReal(const VB_Vector& real,
6463
const VB_Vector& imag, VB_Vector& realIFFT) throw (GenericExcep)
6464
{
6465
6466
/*********************************************************************
6467
* Ensuring that the input VB_Vectors real and imag have the same *
6468
* size. If they don't have the same size, then a GenericExcep is *
6469
* thrown. *
6470
*********************************************************************/
6471
VB_Vector::checkVectorLengths(real.theVector, imag.theVector, __LINE__, __FILE__, __FUNCTION__);
6472
6473
/*********************************************************************
6474
* Ensuring that realIFFT is sized appropriately. *
6475
*********************************************************************/
6476
if (real.getLength() != realIFFT.getLength())
6477
{
6478
realIFFT.resize(real.theVector->size);
6479
} // if
6480
6481
/*********************************************************************
6482
* realIFFTreal and realIFFTimag will hold the real and imaginary *
6483
* parts of the inverse FFT of real, respectively. *
6484
*********************************************************************/
6485
VB_Vector realIFFTreal(real.getLength());
6486
VB_Vector realIFFTimag(real.getLength());
6487
6488
/*********************************************************************
6489
* imagIFFTreal and imagIFFTimag will hold the real and imaginary *
6490
* parts of the inverse FFT of imag, respectively. *
6491
*********************************************************************/
6492
VB_Vector imagIFFTreal(real.getLength());
6493
VB_Vector imagIFFTimag(real.getLength());
6494
6495
/*********************************************************************
6496
* Getting the inverse FFT of real. *
6497
*********************************************************************/
6498
real.ifft(realIFFTreal, realIFFTimag);
6499
6500
/*********************************************************************
6501
* Getting the inverse FFT of imag. *
6502
*********************************************************************/
6503
imag.ifft(imagIFFTreal, imagIFFTimag);
6504
6505
/*********************************************************************
6506
* Now realIFFT will hold the real part of the inverse FFT of *
6507
* (real + i * imag). *
6508
*********************************************************************/
6509
realIFFT = realIFFTreal - imagIFFTimag;
6510
6511
} // void VB_Vector::complexIFFTReal(const VB_Vector& real,
6512
// const VB_Vector& imag, VB_Vector& realIFFT) throw (GenericExcep)
6513
6514
6515
int
6516
VB_Vector::meanNormalize()
6517
{
6518
double avg=this->getVectorMean();
6519
if (fabs(avg) < 1.0) {
6520
if (avg < 0.0) {
6521
*this-=1.0;
6522
avg -= 1.0;
6523
}
6524
else if (avg >= 0.0) {
6525
*this += 1.0;
6526
avg += 1.0;
6527
}
6528
}
6529
*this/=avg;
6530
return 0;
6531
}
6532
6533
int
6534
VB_Vector::removeDrift()
6535
{
6536
double intercept = 0.0, slope = 0.0, cov00 = 0.0, cov01 = 0.0, cov11 = 0.0, chisq = 0.0;
6537
int size = getLength();
6538
double mean = 0;
6539
double x[size], y[size], w[size];
6540
for (int i = 0; i < size; i++) {
6541
x[i] = i;
6542
y[i] = getElement(i);
6543
w[i] = 1.0;
6544
}
6545
gsl_fit_wlinear(x,1,w,1,y,1,size,&intercept,&slope,&cov00,&cov01,&cov11,&chisq);
6546
mean = getVectorMean();
6547
for (int index = 0; index < size; index++)
6548
setElement(index,(getElement(index)-(intercept + slope*index)) + mean);
6549
return 0;
6550
}
6551
6552
double* VB_Vector::begin() const
6553
{
6554
if (this->theVector)
6555
return this->theVector->data;
6556
return NULL;
6557
}
6558
6559
double* VB_Vector::end() const
6560
{
6561
if (this->theVector)
6562
return this->theVector->data + this->theVector->size;
6563
return NULL;
6564
} // double* VB_Vector::end() const
6565
6566
// FIXME : DYK : added the following functions to handle the new i/o system
6567
6568
int
6569
VB_Vector::ReadFile(const string &fname)
6570
{
6571
Vec iovec;
6572
int err;
6573
if ((err=iovec.ReadFile(fname)))
6574
return err;
6575
fileFormat=iovec.fileformat;
6576
init(iovec.size());
6577
memcpy(this->theVector->data, iovec.data, sizeof(double) * this->getLength());
6578
return 0;
6579
}
6580
6581
int
6582
VB_Vector::WriteFile()
6583
{
6584
Vec iovec;
6585
int err;
6586
iovec.resize(getLength());
6587
memcpy(iovec.data, this->theVector->data, sizeof(double) * this->getLength());
6588
iovec.fileformat=fileFormat;
6589
iovec.SetFileName(fileName);
6590
err=iovec.WriteFile();
6591
return err;
6592
}
6593
6594
int
6595
VB_Vector::WriteFile(const string &fname)
6596
{
6597
setFileName(fname);
6598
return WriteFile();
6599
}
6600
6601
// DYK: added the following non-member function to do a t-test
6602
6603
double
6604
ttest(const VB_Vector &v1,const VB_Vector &v2)
6605
{
6606
double s1=v1.getVariance();
6607
double s2=v2.getVariance();
6608
double mean1=v1.getVectorMean();
6609
double mean2=v2.getVectorMean();
6610
return
6611
(mean1-mean2)/sqrt((s1/v1.getLength())+(s2/v2.getLength()));
6612
}
6613
6614
// DYK: added the following non-member function to calculate power
6615
// spectrum
6616
6617
VB_Vector
6618
fftnyquist(VB_Vector &vv)
6619
{
6620
int len=vv.getLength();
6621
VB_Vector fullFFT(len);
6622
vv.getPS(&fullFFT);
6623
6624
int newlen=len/2+1;
6625
VB_Vector halfFFT(newlen);
6626
6627
for (int i=0; i<newlen; i++) {
6628
halfFFT.setElement(i,fullFFT.getElement(i));
6629
}
6630
return halfFFT;
6631
}
6632
6633
6634
// DYK: added to the following better-behaved function for resampling
6635
// vectors
6636
6637
VB_Vector
6638
cspline_resize(VB_Vector vec,double factor)
6639
{
6640
int newsize=(int)((float)vec.size()*factor);
6641
VB_Vector bogus(vec.size());
6642
for (int i=0; i<vec.size(); i++)
6643
bogus.setElement(i,i);
6644
VB_Vector newvector(newsize);
6645
double interval=1.0/factor;
6646
6647
double *xptr=bogus.getTheVector()->data;
6648
double *yptr=vec.getTheVector()->data;
6649
gsl_interp *myinterp=gsl_interp_alloc(gsl_interp_cspline,vec.size());
6650
// gsl_interp *myinterp=gsl_interp_alloc(gsl_interp_linear,vec.size());
6651
gsl_interp_init(myinterp,xptr,yptr,vec.size());
6652
double val,pos=0.0;
6653
for (int i=0; i<newsize; i++) {
6654
val=gsl_interp_eval(myinterp,xptr,yptr,pos,NULL);
6655
newvector.setElement(i,val);
6656
pos+=interval;
6657
}
6658
gsl_interp_free(myinterp);
6659
return newvector;
6660
}
6661
6662
// DYK: added the following useful functions
6663
6664
double
6665
correlation(const VB_Vector &v1,const VB_Vector &v2)
6666
{
6667
if (v1.size()!=v2.size())
6668
return nan("");
6669
double sd1=sqrt(v1.getVariance());
6670
double sd2=sqrt(v2.getVariance());
6671
return (covariance(v1,v2)/(sd1*sd2));
6672
}
6673
6674
double
6675
covariance(const VB_Vector &v1,const VB_Vector &v2)
6676
{
6677
if (v1.size()!=v2.size())
6678
return nan("");
6679
return gsl_stats_covariance(v1.getTheVector()->data,1,
6680
v2.getTheVector()->data,1,
6681
v1.size());
6682
}
6683
6684
// DYK: frequency domain filtering. dfftfilter_frequencies lets you
6685
// specify the frequencies you'd like to remove in Hz
6686
6687
6688
// int
6689
// dfftfilter_frequencies(VB_Vector &signal,VB_Vector filtermask,
6690
// double samplingrate,double low,double high)
6691
// {
6692
// int nfreq=signal.getLength()/2;
6693
// signal.print();
6694
// if (signal.getLength()<4)
6695
// return 101;
6696
// VB_Vector sig_real(signal.getLength());
6697
// VB_Vector sig_imag(signal.getLength());
6698
// signal.fft(sig_real,sig_imag);
6699
6700
// sig_real+=99.0;
6701
6702
// sig_real.print();
6703
// // sig_imag.print();
6704
// //for (int i=0; i<signal.getLength(); i++)
6705
// //printf("freq %d is %.3f\n",i,(double)i/(signal.getLength()*samplingrate));
6706
6707
// for (int i=1; i<2; i++) {
6708
// sig_real(i)=0.0;
6709
// sig_real(16-(i+1))=0.0;
6710
// }
6711
// sig_real.print();
6712
// sig_real[1]=0;
6713
// sig_real[15]=0;
6714
// VB_Vector::complexIFFTReal(sig_real,sig_imag,signal);
6715
// signal.print();
6716
// return 0; // no error!
6717
// }
6718
6719
int
6720
VB_Vector::permute(const VB_Vector &v)
6721
{
6722
if (size()!=v.size())
6723
return 1;
6724
VB_Vector tmp(size());
6725
for (int i=0; i<size(); i++)
6726
tmp[i]=getElement((int)v[i]);
6727
for (int i=0; i<size(); i++)
6728
setElement(i,tmp[i]);
6729
return 0;
6730
}
6731
6732
int
6733
VB_Vector::permute(VBMatrix &m,int col)
6734
{
6735
VB_Vector v;
6736
v=m.GetColumn(col);
6737
return permute(v);
6738
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1