00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "MathVector.h"
00019 #include "MathMatrix.h"
00020 #include "MathConstant.h"
00021 #include "Sort.h"
00022 #include "Error.h"
00023
00024 #ifdef _MSC_VER
00025 #define _USE_MATH_DEFINES
00026 #endif
00027
00028 #include <string.h>
00029 #include <math.h>
00030
00031 int Vector::alloc = 32;
00032
00033 void Vector::Init()
00034 {
00035 dim = size = 0;
00036 label = "Unknown";
00037 data = NULL;
00038 }
00039
00040 Vector::~Vector()
00041 {
00042
00043 if (data != NULL) delete [] data;
00044 }
00045
00046 void Vector::Dimension(int d)
00047 {
00048 if (d > size)
00049 {
00050 if (size < 1024)
00051 {
00052 size = (d + alloc) / alloc * alloc;
00053 double * newData = new double [size];
00054 if (data != NULL)
00055 {
00056 for (int i = 0; i < dim; i++)
00057 newData[i] = data[i];
00058 delete [] data;
00059 }
00060 data = newData;
00061 }
00062 else
00063 {
00064 while (size <= d)
00065 size *= 2;
00066
00067 double * newData = new double [size];
00068 if (data != NULL)
00069 {
00070 for (int i = 0; i < dim; i++)
00071 newData[i] = data[i];
00072 delete [] data;
00073 }
00074 data = newData;
00075 }
00076 }
00077 dim = d;
00078 }
00079
00080 void Vector::Dimension(int d, double value)
00081 {
00082 int original = dim;
00083
00084 Dimension(d);
00085
00086 for (int i = original; i < dim; i++)
00087 data[i] = value;
00088 }
00089
00090 void Vector::Negate()
00091 {
00092 for (int i = 0; i < dim; i++)
00093 data[i] = -data[i];
00094 }
00095
00096 void Vector::Add(double n)
00097 {
00098 for (int i = 0; i< dim; i++)
00099 data[i] += n;
00100 }
00101
00102 void Vector::Multiply(double k)
00103 {
00104 for (int i = 0; i < dim; i++)
00105 data[i] *= k;
00106 }
00107
00108 void Vector::Copy(const Vector & v)
00109 {
00110 Dimension(v.dim);
00111
00112 if (v.data != NULL)
00113 for (int i=0; i < dim; i++)
00114 data[i] = v.data[i];
00115 }
00116
00117 Vector & Vector::operator = (const Vector & rhs)
00118 {
00119 Copy(rhs);
00120 return *this;
00121 }
00122
00123 void Vector::Add(Vector & v)
00124 {
00125 if (dim != v.dim)
00126 error("Vector::Add - vectors have different dimensions\n"
00127 "Vectors - %s [%d] + %s [%d] ",
00128 (const char *) label, dim, (const char *) v.label, v.dim);
00129
00130 for (int i = 0; i < dim; i++)
00131 data[i] += v.data[i];
00132 }
00133
00134 void Vector::AddMultiple(double k, Vector & v)
00135 {
00136 if (dim != v.dim)
00137 error("Vector::AddMultiple - vectors are incompatible\n"
00138 "Vectors - %s [%d] + %s [%d] ",
00139 (const char *) label, dim, (const char *) v.label, v.dim);
00140
00141 for (int i = 0; i < dim; i++)
00142 data[i] += k * v.data[i];
00143 }
00144
00145
00146 void Vector::Subtract(Vector & v)
00147 {
00148 if (dim != v.dim)
00149 error("Vector::Subtract - vectors have different dimensions\n"
00150 "Vectors - %s [%d] + %s [%d] ",
00151 (const char *) label, dim, (const char *) v.label, v.dim);
00152
00153 for (int i = 0; i < dim; i++)
00154 data[i] -= v.data[i];
00155 }
00156
00157
00158 void Vector::Zero()
00159 {
00160 for (int i = 0; i < dim; i++)
00161 data[i] = 0.0;
00162 }
00163
00164 void Vector::Set(double k)
00165 {
00166 for (int i = 0; i < dim; i++)
00167 data[i] = k;
00168 }
00169
00170 void Vector::SetMultiple(double k, Vector & v)
00171 {
00172 Dimension(v.dim);
00173
00174 for (int i = 0; i < dim; i++)
00175 data[i] = k * v[i];
00176 }
00177
00178 double Vector::InnerProduct(Vector & v)
00179 {
00180 if (dim != v.dim)
00181 error("Vector::InnerProduct - vectors have different dimensions\n"
00182 "Vectors - %s[%d] * %s[%d] ",
00183 (const char *) label, dim, (const char *) v.label, v.dim);
00184
00185 double sum = 0.0;
00186 for (int i = 0; i < dim; i++)
00187 sum += data[i] * v.data[i];
00188
00189 return sum;
00190 }
00191
00192 void Vector::Insert(int n, double value)
00193 {
00194 Dimension(dim + 1);
00195
00196 for (int i = dim - 1; i > n; i--)
00197 data[i] = data[i - 1];
00198 data[n] = value;
00199 }
00200
00201 void Vector::DeleteDimension(int n)
00202 {
00203 for (int i = n; i < dim - 1; i++)
00204 data[i] = data[i + 1];
00205 dim --;
00206 }
00207
00208 void Vector::Product(Matrix & m, Vector & v)
00209 {
00210 if (m.cols != v.dim)
00211 error("Vector::Product - Cannot Multiply Matrix by Vector\n"
00212 "Vectors - %s [%d, %d] * %s [%d]\n",
00213 (const char *) m.label, m.rows, m.cols,
00214 (const char *) v.label, v.dim);
00215
00216 Dimension(m.rows);
00217 Zero();
00218
00219 for (int i = 0; i < m.rows; i++)
00220 for (int j = 0; j < m.cols; j++)
00221 data[i] += m[i][j] * v[j];
00222 }
00223
00224 double Vector::Average() const
00225 {
00226 if (dim == 0)
00227 error("Average undefined for null vector %s",
00228 (const char *) label);
00229
00230 return Sum() / dim;
00231 }
00232
00233 double Vector::Product() const
00234 {
00235 double product = 1.0;
00236
00237 for (int j = 0; j < dim; j++)
00238 product *= data[j];
00239
00240 return product;
00241 }
00242
00243 double Vector::Sum() const
00244 {
00245 double sum = 0.0;
00246
00247 for (int j=0; j<dim; j++)
00248 sum += data[j];
00249
00250 return sum;
00251 }
00252
00253 double Vector::SumSquares() const
00254 {
00255 double sum = 0.0;
00256
00257 for (int j=0; j<dim; j++)
00258 sum += data[j] * data[j];
00259
00260 return sum;
00261 }
00262
00263 void Vector::AveVar(double & ave, double & var) const
00264 {
00265
00266
00267
00268 if (dim == 0)
00269 error("Average and Variance undefined for null vector %s",
00270 (const char *) label);
00271
00272 double s, ep;
00273
00274 ave = var = ep = 0.0;
00275
00276 for (int j=0; j<dim; j++)
00277 ave += data[j];
00278
00279 ave /= dim;
00280
00281 for (int j=0; j<dim; j++)
00282 {
00283 s = data[j] - ave;
00284 ep += s;
00285 var += s*s;
00286 }
00287
00288 if (dim > 1)
00289 var = (var - ep*ep/dim)/(dim-1);
00290 }
00291
00292 double Vector::Var() const
00293 {
00294 double mean, var;
00295 AveVar(mean, var);
00296 return var;
00297 }
00298
00299 void Vector::Print(FILE * f, int d)
00300 {
00301 if (d == -1 || d > dim) d = dim;
00302
00303 fprintf(f, "%.15s : ", (const char *) label);
00304 for (int i = 0; i < d; i++)
00305 fprintf(f, "%7.3f ", data[i]);
00306 fprintf(f, "\n");
00307 }
00308
00309 int Vector::CompareDouble(const double * a, const double * b)
00310 {
00311 if (*a < *b) return -1;
00312 if (*a > *b) return 1;
00313 return 0;
00314 }
00315
00316 void Vector::Sort()
00317 {
00318 QuickSort(data, dim, sizeof(double), COMPAREFUNC CompareDouble);
00319 }
00320
00321 void Vector::Sort(Vector & freeRider)
00322 {
00323 QuickSort2(data, freeRider.data, dim, sizeof(double),
00324 COMPAREFUNC CompareDouble);
00325 }
00326
00327 int Vector::BinarySearch(double element)
00328 {
00329 void * pointer = ::BinarySearch
00330 (&element, data, dim, sizeof(double), COMPAREFUNC CompareDouble);
00331
00332 if (pointer == NULL)
00333 return -1;
00334
00335 return ((double *) pointer) - data;
00336 }
00337
00338 void Vector::RemoveDuplicates()
00339 {
00340 int out = 0;
00341
00342 for (int in = 1; in < Length(); in++)
00343 if (data[in] != data[out])
00344 data[++out] = data[in];
00345
00346 Dimension(out + 1);
00347 }
00348
00349 bool Vector::operator == (const Vector & rhs) const
00350 {
00351 if (rhs.dim != dim) return false;
00352
00353 for (int i = 0; i < dim; i++)
00354 if (data[i] != rhs[i])
00355 return false;
00356 return true;
00357 }
00358
00359
00360
00361
00362 int Vector::CountIfGreater(double threshold) const
00363 {
00364 int count = 0;
00365
00366 for (int i = 0; i < dim; i++)
00367 if (data[i] > threshold)
00368 count++;
00369
00370 return count;
00371 }
00372
00373 int Vector::CountIfGreaterOrEqual(double treshold) const
00374 {
00375 int count = 0;
00376
00377 for (int i = 0; i < dim; i++)
00378 if (data[i] >= treshold)
00379 count++;
00380
00381 return count;
00382 }
00383
00384
00385
00386
00387 double Vector::Min() const
00388 {
00389 if (dim == 0)
00390 return 0.0;
00391
00392 double min = data[0];
00393
00394 for (int i = 1; i < dim; i++)
00395 if (data[i] < min)
00396 min = data[i];
00397
00398 return min;
00399 }
00400
00401 double Vector::Max() const
00402 {
00403 if (dim == 0)
00404 return 0.0;
00405
00406 double max = data[0];
00407
00408 for (int i = 1; i < dim; i++)
00409 if (data[i] > max)
00410 max = data[i];
00411
00412 return max;
00413 }
00414
00415
00416
00417
00418 void Vector::Push(double value)
00419 {
00420 Dimension(dim + 1);
00421 data[dim - 1] = value;
00422 }
00423
00424 void Vector::Stack(const Vector & v)
00425 {
00426 int end = dim;
00427
00428 Dimension(dim + v.dim);
00429
00430 for (int i = 0; i < v.dim; i++)
00431 data[i + end] = v[i];
00432 }
00433
00434
00435
00436
00437 bool Vector::isAscending()
00438 {
00439 for (int i = 1; i < dim; i++)
00440 if (data[i] < data[i - 1])
00441 return false;
00442 return true;
00443 }
00444
00445 bool Vector::isDescending()
00446 {
00447 for (int i = 1; i < dim; i++)
00448 if (data[i] > data[i - 1])
00449 return false;
00450 return true;
00451 }
00452
00453
00454
00455
00456 VectorFunc::VectorFunc()
00457 {
00458 f = NULL;
00459 }
00460
00461 VectorFunc::VectorFunc(double(*func)(Vector &))
00462 {
00463 f = func;
00464 }
00465
00466 double VectorFunc::Evaluate(Vector & v)
00467 {
00468 return f(v);
00469 }
00470
00471 #ifndef M_SQRT2
00472 #define M_SQRT2 1.41421356
00473 #endif
00474
00475 #define MAXROUNDS 10
00476 #define SQRT_HALF (1.0/M_SQRT2)
00477 #define TWO (M_SQRT2 * M_SQRT2)
00478
00479 void VectorFunc::Derivative(Vector & x, Vector & d, double h_start)
00480 {
00481 double a[MAXROUNDS][MAXROUNDS];
00482
00483
00484 for (int k = 0; k < x.dim; k++)
00485 {
00486 double left, right;
00487 double save_x = x[k];
00488 double h = h_start;
00489
00490
00491 x[k] = save_x - h;
00492 left = Evaluate(x);
00493
00494
00495 if (k == 0 || left < dfmin)
00496 dfmin = left, dpmin = x;
00497
00498
00499 x[k] = save_x + h;
00500 right = Evaluate(x);
00501
00502
00503 if (right < dfmin)
00504 dfmin = left, dpmin = x;
00505
00506
00507 a[0][0] = (right - left) / (2.0 * h);
00508
00509
00510 double err = 1e30;
00511
00512
00513
00514 for (int i = 1; i < MAXROUNDS; i++)
00515 {
00516
00517 h *= SQRT_HALF;
00518
00519
00520 x[k] = save_x - h;
00521 left = Evaluate(x);
00522 if (left < dfmin) dfmin = left, dpmin = x;
00523 x[k] = save_x + h;
00524 right = Evaluate(x);
00525 if (right < dfmin) dfmin = right, dpmin = x;
00526
00527
00528 a[0][i] = (right - left) / (2.0 * h);
00529
00530
00531 double factor = TWO;
00532
00533 for (int j = 1; j <= i; j++)
00534 {
00535 a[j][i] = (a[j-1][i] * factor - a[j-1][i-1])/(factor - 1.0);
00536
00537 factor *= TWO;
00538
00539 double error = max(fabs(a[j][i] - a[j-1][i]), fabs(a[j][i] - a[j-1][i-1]));
00540
00541
00542 if (error < err)
00543 {
00544 err = error;
00545 d[k] = a[j][i];
00546 }
00547 }
00548
00549
00550 if (fabs(a[i][i] - a[i-1][i-1]) >= 2.0 * err)
00551 {
00552 x[k] = save_x;
00553 break;
00554 }
00555 }
00556
00557 x[k] = save_x;
00558 }
00559 }
00560
00561 int Vector::SafeCount() const
00562 {
00563 int nonMissing = dim;
00564
00565 for (int i = 0; i < dim; i++)
00566 if (data[i] == _NAN_)
00567 nonMissing--;
00568
00569 return nonMissing;
00570 }
00571
00572 double Vector::SafeMin() const
00573 {
00574 double min = _NAN_;
00575 int i;
00576
00577 for (i = 0; i < dim; i++)
00578 if (data[i] != _NAN_)
00579 {
00580 min = data[i];
00581 break;
00582 }
00583
00584 for (; i < dim; i++)
00585 if (data[i] != _NAN_ && data[i] < min)
00586 min = data[i];
00587
00588 return min;
00589 }
00590
00591 double Vector::SafeMax() const
00592 {
00593 double max = _NAN_;
00594 int i;
00595
00596 for (i = 0; i < dim; i++)
00597 if (data[i] != _NAN_)
00598 {
00599 max = data[i];
00600 break;
00601 }
00602
00603 for (; i < dim; i++)
00604 if (data[i] != _NAN_ && data[i] > max)
00605 max = data[i];
00606
00607 return max;
00608 }
00609
00610 void Vector::Reverse()
00611 {
00612 for (int i = 0, j = dim - 1; i < j; i++, j--)
00613 Swap(i, j);
00614 }
00615
00616 void Vector::InsertInSortedList(int value)
00617 {
00618
00619 int pos = dim - 1;
00620
00621 while (pos >= 0 && data[pos] > value)
00622 pos--;
00623
00624
00625 if (pos >= 0 && data[pos] == value)
00626 return;
00627
00628
00629 Dimension(dim + 1);
00630
00631
00632 pos++;
00633 for (int i = dim - 1; i > pos; i--)
00634 data[i] = data[i - 1];
00635
00636 data[pos] = value;
00637 }
00638
00639 void Vector::Swap(Vector & rhs)
00640 {
00641 double * temp = rhs.data;
00642 rhs.data = data;
00643 data = temp;
00644
00645 int swap = rhs.dim;
00646 rhs.dim = dim;
00647 dim = swap;
00648
00649 swap = rhs.size;
00650 rhs.size = size;
00651 size = swap;
00652 }
00653
00654 double Vector::Average(double returnIfNull)
00655 {
00656 if (Length() == 0)
00657 return returnIfNull;
00658
00659 return Average();
00660 }
00661
00662 double Vector::Var(double returnIfNull)
00663 {
00664 if (Length() == 0)
00665 return returnIfNull;
00666
00667 return Var();
00668 }
00669