libStatGen Software  1
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 
31 int Vector::alloc = 32;
32 
33 void Vector::Init()
34 {
35  dim = size = 0;
36  label = "Unknown";
37  data = NULL;
38 }
39 
40 Vector::~Vector()
41 {
42  // printf(" Deleting vector %s ...\n", (const char *) label);
43  if (data != NULL) delete [] data;
44 }
45 
46 void 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 
80 void 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 
90 void Vector::Negate()
91 {
92  for (int i = 0; i < dim; i++)
93  data[i] = -data[i];
94 }
95 
96 void Vector::Add(double n)
97 {
98  for (int i = 0; i< dim; i++)
99  data[i] += n;
100 }
101 
102 void Vector::Multiply(double k)
103 {
104  for (int i = 0; i < dim; i++)
105  data[i] *= k;
106 }
107 
108 void 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 
117 Vector & Vector::operator = (const Vector & rhs)
118 {
119  Copy(rhs);
120  return *this;
121 }
122 
123 void 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 
134 void 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 
146 void 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 
158 void Vector::Zero()
159 {
160  for (int i = 0; i < dim; i++)
161  data[i] = 0.0;
162 }
163 
164 void Vector::Set(double k)
165 {
166  for (int i = 0; i < dim; i++)
167  data[i] = k;
168 }
169 
170 void 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 
178 double 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 
192 void 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 
201 void Vector::DeleteDimension(int n)
202 {
203  for (int i = n; i < dim - 1; i++)
204  data[i] = data[i + 1];
205  dim --;
206 }
207 
208 void 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 
224 double 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 
233 double 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 
243 double 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 
253 double 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 
263 void 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 
292 double Vector::Var() const
293 {
294  double mean, var;
295  AveVar(mean, var);
296  return var;
297 }
298 
299 double Vector::StandardDeviation() const
300 {
301  double var = Var();
302 
303  if (var < 0.0) var = 0.0;
304 
305  return sqrt(var);
306 }
307 
308 void 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 
318 int 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 
325 void Vector::Sort()
326 {
327  QuickSort(data, dim, sizeof(double), COMPAREFUNC CompareDouble);
328 }
329 
330 void Vector::Sort(Vector & freeRider)
331 {
332  QuickSort2(data, freeRider.data, dim, sizeof(double),
333  COMPAREFUNC CompareDouble);
334 }
335 
336 int 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 
347 void 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 
358 bool 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 
371 int 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 
382 int 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 
396 double 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 
410 double 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 
427 void Vector::Push(double value)
428 {
429  Dimension(dim + 1);
430  data[dim - 1] = value;
431 }
432 
433 void 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 
446 bool 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 
454 bool 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 
465 VectorFunc::VectorFunc()
466 {
467  f = NULL;
468 }
469 
470 VectorFunc::VectorFunc(double(*func)(Vector &))
471 {
472  f = func;
473 }
474 
475 double 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 
488 void 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 
570 int 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 
581 double 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 
600 double 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 
619 void Vector::Reverse()
620 {
621  for (int i = 0, j = dim - 1; i < j; i++, j--)
622  Swap(i, j);
623 }
624 
625 void 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 
648 void 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 
663 double Vector::Average(double returnIfNull)
664 {
665  if (Length() == 0)
666  return returnIfNull;
667 
668  return Average();
669 }
670 
671 double Vector::Var(double returnIfNull)
672 {
673  if (Length() == 0)
674  return returnIfNull;
675 
676  return Var();
677 }
678 
679 double Vector::StandardDeviation(double returnIfNull)
680 {
681  if (Length() == 0)
682  return returnIfNull;
683 
684  return StandardDeviation();
685 }
Matrix
Definition: MathMatrix.h:77
Vector
Definition: MathVector.h:29