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 
Init()33 void Vector::Init()
34 {
35     dim = size = 0;
36     label = "Unknown";
37     data = NULL;
38 }
39 
~Vector()40 Vector::~Vector()
41 {
42     // printf(" Deleting vector %s ...\n", (const char *) label);
43     if (data != NULL) delete [] data;
44 }
45 
Dimension(int d)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 
Dimension(int d,double value)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 
Negate()90 void Vector::Negate()
91 {
92     for (int i = 0; i < dim; i++)
93         data[i] = -data[i];
94 }
95 
Add(double n)96 void Vector::Add(double n)
97 {
98     for (int i = 0; i< dim; i++)
99         data[i] += n;
100 }
101 
Multiply(double k)102 void Vector::Multiply(double k)
103 {
104     for (int i = 0; i < dim; i++)
105         data[i] *= k;
106 }
107 
Copy(const Vector & v)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 
operator =(const Vector & rhs)117 Vector & Vector::operator = (const Vector & rhs)
118 {
119     Copy(rhs);
120     return *this;
121 }
122 
Add(Vector & v)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 
AddMultiple(double k,Vector & v)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 
Subtract(Vector & v)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 
Zero()158 void Vector::Zero()
159 {
160     for (int i = 0; i < dim; i++)
161         data[i] = 0.0;
162 }
163 
Set(double k)164 void Vector::Set(double k)
165 {
166     for (int i = 0; i < dim; i++)
167         data[i] = k;
168 }
169 
SetMultiple(double k,Vector & v)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 
InnerProduct(Vector & v)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 
Insert(int n,double value)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 
DeleteDimension(int n)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 
Product(Matrix & m,Vector & v)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 
Average() const224 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 
Product() const233 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 
Sum() const243 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 
SumSquares() const253 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 
AveVar(double & ave,double & var) const263 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 
Var() const292 double Vector::Var() const
293 {
294     double mean, var;
295     AveVar(mean, var);
296     return var;
297 }
298 
StandardDeviation() const299 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 
Print(FILE * f,int d)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 
CompareDouble(const double * a,const double * b)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 
Sort()325 void Vector::Sort()
326 {
327     QuickSort(data, dim, sizeof(double), COMPAREFUNC CompareDouble);
328 }
329 
Sort(Vector & freeRider)330 void Vector::Sort(Vector & freeRider)
331 {
332     QuickSort2(data, freeRider.data, dim, sizeof(double),
333                COMPAREFUNC CompareDouble);
334 }
335 
BinarySearch(double element)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 
RemoveDuplicates()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 
operator ==(const Vector & rhs) const358 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 
CountIfGreater(double threshold) const371 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 
CountIfGreaterOrEqual(double treshold) const382 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 
Min() const396 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 
Max() const410 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 
Push(double value)427 void Vector::Push(double value)
428 {
429     Dimension(dim + 1);
430     data[dim - 1] = value;
431 }
432 
Stack(const Vector & v)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 
isAscending()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 
isDescending()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 
VectorFunc()465 VectorFunc::VectorFunc()
466 {
467     f = NULL;
468 }
469 
VectorFunc(double (* func)(Vector &))470 VectorFunc::VectorFunc(double(*func)(Vector &))
471 {
472     f = func;
473 }
474 
Evaluate(Vector & v)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 
Derivative(Vector & x,Vector & d,double h_start)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 
SafeCount() const570 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 
SafeMin() const581 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 
SafeMax() const600 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 
Reverse()619 void Vector::Reverse()
620 {
621     for (int i = 0, j = dim - 1; i < j; i++, j--)
622         Swap(i, j);
623 }
624 
InsertInSortedList(int value)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 
Swap(Vector & rhs)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 
Average(double returnIfNull)663 double Vector::Average(double returnIfNull)
664 {
665     if (Length() == 0)
666         return returnIfNull;
667 
668     return Average();
669 }
670 
Var(double returnIfNull)671 double Vector::Var(double returnIfNull)
672 {
673     if (Length() == 0)
674         return returnIfNull;
675 
676     return Var();
677 }
678 
StandardDeviation(double returnIfNull)679 double Vector::StandardDeviation(double returnIfNull)
680 {
681     if (Length() == 0)
682         return returnIfNull;
683 
684     return StandardDeviation();
685 }
686