libStatGen Software 1
Loading...
Searching...
No Matches
MathVector.cpp
1/*
2 * Copyright (C) 2010 Regents of the University of Michigan
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18#include "MathVector.h"
19#include "MathMatrix.h"
20#include "MathConstant.h"
21#include "Sort.h"
22#include "Error.h"
23
24#ifdef _MSC_VER
25#define _USE_MATH_DEFINES
26#endif
27
28#include <string.h>
29#include <math.h>
30
31int Vector::alloc = 32;
32
33void Vector::Init()
34{
35 dim = size = 0;
36 label = "Unknown";
37 data = NULL;
38}
39
40Vector::~Vector()
41{
42 // printf(" Deleting vector %s ...\n", (const char *) label);
43 if (data != NULL) delete [] data;
44}
45
46void Vector::Dimension(int d)
47{
48 if (d > size)
49 {
50 if (size < 1024)
51 {
52 size = (d + alloc) / alloc * alloc;
53 double * newData = new double [size];
54 if (data != NULL)
55 {
56 for (int i = 0; i < dim; i++)
57 newData[i] = data[i];
58 delete [] data;
59 }
60 data = newData;
61 }
62 else
63 {
64 while (size <= d)
65 size *= 2;
66
67 double * newData = new double [size];
68 if (data != NULL)
69 {
70 for (int i = 0; i < dim; i++)
71 newData[i] = data[i];
72 delete [] data;
73 }
74 data = newData;
75 }
76 }
77 dim = d;
78}
79
80void Vector::Dimension(int d, double value)
81{
82 int original = dim;
83
84 Dimension(d);
85
86 for (int i = original; i < dim; i++)
87 data[i] = value;
88}
89
90void Vector::Negate()
91{
92 for (int i = 0; i < dim; i++)
93 data[i] = -data[i];
94}
95
96void Vector::Add(double n)
97{
98 for (int i = 0; i< dim; i++)
99 data[i] += n;
100}
101
102void Vector::Multiply(double k)
103{
104 for (int i = 0; i < dim; i++)
105 data[i] *= k;
106}
107
108void Vector::Copy(const Vector & v)
109{
110 Dimension(v.dim);
111
112 if (v.data != NULL)
113 for (int i=0; i < dim; i++)
114 data[i] = v.data[i];
115}
116
117Vector & Vector::operator = (const Vector & rhs)
118{
119 Copy(rhs);
120 return *this;
121}
122
123void Vector::Add(Vector & v)
124{
125 if (dim != v.dim)
126 error("Vector::Add - vectors have different dimensions\n"
127 "Vectors - %s [%d] + %s [%d] ",
128 (const char *) label, dim, (const char *) v.label, v.dim);
129
130 for (int i = 0; i < dim; i++)
131 data[i] += v.data[i];
132}
133
134void Vector::AddMultiple(double k, Vector & v)
135{
136 if (dim != v.dim)
137 error("Vector::AddMultiple - vectors are incompatible\n"
138 "Vectors - %s [%d] + %s [%d] ",
139 (const char *) label, dim, (const char *) v.label, v.dim);
140
141 for (int i = 0; i < dim; i++)
142 data[i] += k * v.data[i];
143}
144
145
146void Vector::Subtract(Vector & v)
147{
148 if (dim != v.dim)
149 error("Vector::Subtract - vectors have different dimensions\n"
150 "Vectors - %s [%d] + %s [%d] ",
151 (const char *) label, dim, (const char *) v.label, v.dim);
152
153 for (int i = 0; i < dim; i++)
154 data[i] -= v.data[i];
155}
156
157
158void Vector::Zero()
159{
160 for (int i = 0; i < dim; i++)
161 data[i] = 0.0;
162}
163
164void Vector::Set(double k)
165{
166 for (int i = 0; i < dim; i++)
167 data[i] = k;
168}
169
170void Vector::SetMultiple(double k, Vector & v)
171{
172 Dimension(v.dim);
173
174 for (int i = 0; i < dim; i++)
175 data[i] = k * v[i];
176}
177
178double Vector::InnerProduct(Vector & v)
179{
180 if (dim != v.dim)
181 error("Vector::InnerProduct - vectors have different dimensions\n"
182 "Vectors - %s[%d] * %s[%d] ",
183 (const char *) label, dim, (const char *) v.label, v.dim);
184
185 double sum = 0.0;
186 for (int i = 0; i < dim; i++)
187 sum += data[i] * v.data[i];
188
189 return sum;
190}
191
192void Vector::Insert(int n, double value)
193{
194 Dimension(dim + 1);
195
196 for (int i = dim - 1; i > n; i--)
197 data[i] = data[i - 1];
198 data[n] = value;
199}
200
201void Vector::DeleteDimension(int n)
202{
203 for (int i = n; i < dim - 1; i++)
204 data[i] = data[i + 1];
205 dim --;
206}
207
208void Vector::Product(Matrix & m, Vector & v)
209{
210 if (m.cols != v.dim)
211 error("Vector::Product - Cannot Multiply Matrix by Vector\n"
212 "Vectors - %s [%d, %d] * %s [%d]\n",
213 (const char *) m.label, m.rows, m.cols,
214 (const char *) v.label, v.dim);
215
216 Dimension(m.rows);
217 Zero();
218
219 for (int i = 0; i < m.rows; i++)
220 for (int j = 0; j < m.cols; j++)
221 data[i] += m[i][j] * v[j];
222}
223
224double Vector::Average() const
225{
226 if (dim == 0)
227 error("Average undefined for null vector %s",
228 (const char *) label);
229
230 return Sum() / dim;
231}
232
233double Vector::Product() const
234{
235 double product = 1.0;
236
237 for (int j = 0; j < dim; j++)
238 product *= data[j];
239
240 return product;
241}
242
243double Vector::Sum() const
244{
245 double sum = 0.0;
246
247 for (int j=0; j<dim; j++)
248 sum += data[j];
249
250 return sum;
251}
252
253double Vector::SumSquares() const
254{
255 double sum = 0.0;
256
257 for (int j=0; j<dim; j++)
258 sum += data[j] * data[j];
259
260 return sum;
261}
262
263void Vector::AveVar(double & ave, double & var) const
264{
265 // uses a two pass method to correct for
266 // round-off errors
267
268 if (dim == 0)
269 error("Average and Variance undefined for null vector %s",
270 (const char *) label);
271
272 double s, ep;
273
274 ave = var = ep = 0.0;
275
276 for (int j=0; j<dim; j++)
277 ave += data[j];
278
279 ave /= dim;
280
281 for (int j=0; j<dim; j++)
282 {
283 s = data[j] - ave;
284 ep += s;
285 var += s*s;
286 }
287
288 if (dim > 1)
289 var = (var - ep*ep/dim)/(dim-1);
290}
291
292double Vector::Var() const
293{
294 double mean, var;
295 AveVar(mean, var);
296 return var;
297}
298
299double Vector::StandardDeviation() const
300{
301 double var = Var();
302
303 if (var < 0.0) var = 0.0;
304
305 return sqrt(var);
306}
307
308void Vector::Print(FILE * f, int d)
309{
310 if (d == -1 || d > dim) d = dim;
311
312 fprintf(f, "%.15s : ", (const char *) label);
313 for (int i = 0; i < d; i++)
314 fprintf(f, "%7.3f ", data[i]);
315 fprintf(f, "\n");
316}
317
318int Vector::CompareDouble(const double * a, const double * b)
319{
320 if (*a < *b) return -1;
321 if (*a > *b) return 1;
322 return 0;
323}
324
325void Vector::Sort()
326{
327 QuickSort(data, dim, sizeof(double), COMPAREFUNC CompareDouble);
328}
329
330void Vector::Sort(Vector & freeRider)
331{
332 QuickSort2(data, freeRider.data, dim, sizeof(double),
333 COMPAREFUNC CompareDouble);
334}
335
336int Vector::BinarySearch(double element)
337{
338 void * pointer = ::BinarySearch
339 (&element, data, dim, sizeof(double), COMPAREFUNC CompareDouble);
340
341 if (pointer == NULL)
342 return -1;
343
344 return ((double *) pointer) - data;
345}
346
347void Vector::RemoveDuplicates()
348{
349 int out = 0;
350
351 for (int in = 1; in < Length(); in++)
352 if (data[in] != data[out])
353 data[++out] = data[in];
354
355 Dimension(out + 1);
356}
357
358bool Vector::operator == (const Vector & rhs) const
359{
360 if (rhs.dim != dim) return false;
361
362 for (int i = 0; i < dim; i++)
363 if (data[i] != rhs[i])
364 return false;
365 return true;
366}
367
368// These functions are useful for simulation
369//
370
371int Vector::CountIfGreater(double threshold) const
372{
373 int count = 0;
374
375 for (int i = 0; i < dim; i++)
376 if (data[i] > threshold)
377 count++;
378
379 return count;
380}
381
382int Vector::CountIfGreaterOrEqual(double treshold) const
383{
384 int count = 0;
385
386 for (int i = 0; i < dim; i++)
387 if (data[i] >= treshold)
388 count++;
389
390 return count;
391}
392
393// Min and max functions
394//
395
396double Vector::Min() const
397{
398 if (dim == 0)
399 return 0.0;
400
401 double min = data[0];
402
403 for (int i = 1; i < dim; i++)
404 if (data[i] < min)
405 min = data[i];
406
407 return min;
408}
409
410double Vector::Max() const
411{
412 if (dim == 0)
413 return 0.0;
414
415 double max = data[0];
416
417 for (int i = 1; i < dim; i++)
418 if (data[i] > max)
419 max = data[i];
420
421 return max;
422}
423
424// Push and Pop functions for using vector as a stack
425//
426
427void Vector::Push(double value)
428{
429 Dimension(dim + 1);
430 data[dim - 1] = value;
431}
432
433void Vector::Stack(const Vector & v)
434{
435 int end = dim;
436
437 Dimension(dim + v.dim);
438
439 for (int i = 0; i < v.dim; i++)
440 data[i + end] = v[i];
441}
442
443// Check if all values are in ascending or descending order
444//
445
446bool Vector::isAscending()
447{
448 for (int i = 1; i < dim; i++)
449 if (data[i] < data[i - 1])
450 return false;
451 return true;
452}
453
454bool Vector::isDescending()
455{
456 for (int i = 1; i < dim; i++)
457 if (data[i] > data[i - 1])
458 return false;
459 return true;
460}
461
462// VectorFunc class
463//
464
465VectorFunc::VectorFunc()
466{
467 f = NULL;
468}
469
470VectorFunc::VectorFunc(double(*func)(Vector &))
471{
472 f = func;
473}
474
475double VectorFunc::Evaluate(Vector & v)
476{
477 return f(v);
478}
479
480#ifndef M_SQRT2
481#define M_SQRT2 1.41421356
482#endif
483
484#define MAXROUNDS 10
485#define SQRT_HALF (1.0/M_SQRT2)
486#define TWO (M_SQRT2 * M_SQRT2)
487
488void VectorFunc::Derivative(Vector & x, Vector & d, double h_start)
489{
490 double a[MAXROUNDS][MAXROUNDS];
491
492 // Calculate derivatives along each direction ...
493 for (int k = 0; k < x.dim; k++)
494 {
495 double left, right;
496 double save_x = x[k];
497 double h = h_start;
498
499 // Evaluate function to the left of x along direction k
500 x[k] = save_x - h;
501 left = Evaluate(x);
502
503 // Initialize or update dfmin if appropriate...
504 if (k == 0 || left < dfmin)
505 dfmin = left, dpmin = x;
506
507 // Evaluate function to the right of x along direction k
508 x[k] = save_x + h;
509 right = Evaluate(x);
510
511 // Update dfmin
512 if (right < dfmin)
513 dfmin = left, dpmin = x;
514
515 // Initial crude estimate
516 a[0][0] = (right - left) / (2.0 * h);
517
518 // Initial guess of error is large
519 double err = 1e30;
520
521 // At each round, update Neville tableau with smaller stepsize and higher
522 // order extrapolation ...
523 for (int i = 1; i < MAXROUNDS; i++)
524 {
525 // Decrease h
526 h *= SQRT_HALF;
527
528 // Re-evaluate function and update dfmin as required
529 x[k] = save_x - h;
530 left = Evaluate(x);
531 if (left < dfmin) dfmin = left, dpmin = x;
532 x[k] = save_x + h;
533 right = Evaluate(x);
534 if (right < dfmin) dfmin = right, dpmin = x;
535
536 // Improved estimate of derivative
537 a[0][i] = (right - left) / (2.0 * h);
538
539 // Calculate extrapolations of various orders ...
540 double factor = TWO;
541
542 for (int j = 1; j <= i; j++)
543 {
544 a[j][i] = (a[j-1][i] * factor - a[j-1][i-1])/(factor - 1.0);
545
546 factor *= TWO;
547
548 double error = max(fabs(a[j][i] - a[j-1][i]), fabs(a[j][i] - a[j-1][i-1]));
549
550 // Did we improve solution?
551 if (error < err)
552 {
553 err = error;
554 d[k] = a[j][i];
555 }
556 }
557
558 // Stop if solution is deteriorating ...
559 if (fabs(a[i][i] - a[i-1][i-1]) >= 2.0 * err)
560 {
561 x[k] = save_x;
562 break;
563 }
564 }
565
566 x[k] = save_x;
567 }
568}
569
570int Vector::SafeCount() const
571{
572 int nonMissing = dim;
573
574 for (int i = 0; i < dim; i++)
575 if (data[i] == _NAN_)
576 nonMissing--;
577
578 return nonMissing;
579}
580
581double Vector::SafeMin() const
582{
583 double min = _NAN_;
584 int i;
585
586 for (i = 0; i < dim; i++)
587 if (data[i] != _NAN_)
588 {
589 min = data[i];
590 break;
591 }
592
593 for (; i < dim; i++)
594 if (data[i] != _NAN_ && data[i] < min)
595 min = data[i];
596
597 return min;
598}
599
600double Vector::SafeMax() const
601{
602 double max = _NAN_;
603 int i;
604
605 for (i = 0; i < dim; i++)
606 if (data[i] != _NAN_)
607 {
608 max = data[i];
609 break;
610 }
611
612 for (; i < dim; i++)
613 if (data[i] != _NAN_ && data[i] > max)
614 max = data[i];
615
616 return max;
617}
618
619void Vector::Reverse()
620{
621 for (int i = 0, j = dim - 1; i < j; i++, j--)
622 Swap(i, j);
623}
624
625void Vector::InsertInSortedList(int value)
626{
627 // Skip through large elements
628 int pos = dim - 1;
629
630 while (pos >= 0 && data[pos] > value)
631 pos--;
632
633 // If the value is already in the list, we are done
634 if (pos >= 0 && data[pos] == value)
635 return;
636
637 // Otherwise we need to grow array
638 Dimension(dim + 1);
639
640 // And then shift larger elements to the right
641 pos++;
642 for (int i = dim - 1; i > pos; i--)
643 data[i] = data[i - 1];
644
645 data[pos] = value;
646}
647
648void Vector::Swap(Vector & rhs)
649{
650 double * temp = rhs.data;
651 rhs.data = data;
652 data = temp;
653
654 int swap = rhs.dim;
655 rhs.dim = dim;
656 dim = swap;
657
658 swap = rhs.size;
659 rhs.size = size;
660 size = swap;
661}
662
663double Vector::Average(double returnIfNull)
664{
665 if (Length() == 0)
666 return returnIfNull;
667
668 return Average();
669}
670
671double Vector::Var(double returnIfNull)
672{
673 if (Length() == 0)
674 return returnIfNull;
675
676 return Var();
677}
678
679double Vector::StandardDeviation(double returnIfNull)
680{
681 if (Length() == 0)
682 return returnIfNull;
683
684 return StandardDeviation();
685}