LinAlg: VecMat.cpp Source File

00001 //===========================================================================
00048 //===========================================================================
00049 
00050 
00051 
00052 #include <LinAlg/VecMat.h>
00053 #include <LinAlg/LinAlg.h>
00054 
00055 
00056 Vector::Vector()
00057 : Array<double>(0u)
00058 {
00059 }
00060 
00061 Vector::Vector(unsigned int len, bool zero)
00062 : Array<double>(len)
00063 {
00064     if (zero) (*(Array<double>*)this) = 0.0;
00065 }
00066 
00067 Vector::Vector(const Array<double>& other)
00068 : Array<double>(other)
00069 {
00070     SIZE_CHECK(ndim() == 1);
00071 }
00072 
00073 Vector::Vector(const Vector& other)
00074 : Array<double>(other)
00075 {
00076 }
00077 
00078 Vector::~Vector()
00079 {
00080 }
00081 
00082 
00083 double Vector::norm2() const
00084 {
00085     double ret = 0.0;
00086     int i, ic = dim(0);
00087     for (i=0; i<ic; i++)
00088     {
00089         double val = operator () (i);
00090         ret += val * val;
00091     }
00092     return ret;
00093 }
00094 
00095 bool Vector::normalize()
00096 {
00097     double n = norm();
00098     if (n == 0.0) return false;
00099     int i, ic = dim(0);
00100     for (i=0; i<ic; i++) operator () (i) /= n;
00101     return true;
00102 }
00103 
00104 ArrayBase* Vector::clone() const
00105 {
00106     return new Vector(*this);
00107 }
00108 
00109 ArrayBase* Vector::cloneEmpty() const
00110 {
00111     return new Vector(0, false);
00112 }
00113 
00114 Vector::Vector(Array<double>& other, bool reference)
00115 : Array<double>(other.dimvec(), other.elemvec(), other.ndim(), other.nelem(), reference)
00116 {
00117 }
00118 
00119 Vector::Vector(unsigned int* d, double* e, unsigned int nd, unsigned int ne, bool reference)
00120 : Array<double>(d, e, nd, ne, reference)
00121 {
00122 }
00123 
00124 void Vector::resize_i(unsigned int* d, unsigned int nd, bool copy)
00125 {
00126     SIZE_CHECK(nd == 1);
00127     Array<double>::resize_i(d, 1, copy);
00128 }
00129 
00130 
00132 
00133 
00134 Matrix::Matrix()
00135 : Array<double>(0u, 0u)
00136 {
00137 }
00138 
00139 Matrix::Matrix(unsigned int rows, unsigned int cols, bool zero)
00140 : Array<double>(rows, cols)
00141 {
00142     if (zero) (*(Array<double>*)this) = 0.0;
00143 }
00144 
00145 Matrix::Matrix(const Array<double>& other)
00146 : Array<double>(other)
00147 {
00148     SIZE_CHECK(ndim() == 2);
00149 }
00150 
00151 Matrix::Matrix(const Matrix& other)
00152 : Array<double>(other)
00153 {
00154     SIZE_CHECK(ndim() == 2);
00155 }
00156 
00157 Matrix::~Matrix()
00158 {
00159 }
00160 
00161 
00162 // static
00163 Matrix Matrix::unitmatrix(unsigned int dim)
00164 {
00165     Matrix ret(dim, dim, true);
00166     unsigned int i;
00167     for (i=0; i<dim; i++) ret(i, i) = 1.0;
00168     return ret;
00169 }
00170 
00171 // static
00172 Matrix Matrix::diagonal(const Array<double>& diag)
00173 {
00174     SIZE_CHECK(diag.ndim() == 1);
00175     unsigned int dim = diag.dim(0);
00176     Matrix ret(dim, dim, true);
00177     unsigned int i;
00178     for (i=0; i<dim; i++) ret(i, i) = 1.0;
00179     return ret;
00180 }
00181 
00182 Matrix Matrix::transpose() const
00183 {
00184     SIZE_CHECK(ndim() == 2);
00185     Matrix ret(*this);
00186 	::transpose(ret);
00187     return ret;
00188 }
00189 
00190 Matrix Matrix::inverse(unsigned maxIterations, double tolerance, bool ignoreThreshold) const
00191 {
00192     SIZE_CHECK(ndim() == 2);
00193     Matrix ret(dim(1), dim(0));
00194     g_inverse(*this, ret, maxIterations, tolerance, ignoreThreshold);
00195     return ret;
00196 }
00197 
00198 Matrix Matrix::inverseCholesky(double thresholdFactor) const
00199 {
00200     SIZE_CHECK(ndim() == 2);
00201     Matrix ret(dim(1), dim(0));
00202     g_inverseCholesky(*this, ret,thresholdFactor);
00203     return ret;
00204 }
00205 
00206 Matrix Matrix::inverseMoorePenrose() const
00207 {
00208     SIZE_CHECK(ndim() == 2);
00209     Matrix ret(dim(1), dim(0));
00210     g_inverseMoorePenrose(*this, ret);
00211     return ret;
00212 }
00213 
00214 Matrix Matrix::inverseSymm() const
00215 {
00216     SIZE_CHECK(ndim() == 2);
00217     Matrix ret(dim(1), dim(0));
00218     invertSymm(ret, *this);
00219     return ret;
00220 }
00221 
00222 Matrix Matrix::inverseSymmPositiveDefinite() const
00223 {
00224     SIZE_CHECK(ndim() == 2);
00225     Matrix ret(dim(1), dim(0));
00226     invertSymmPositiveDefinite(ret, *this);
00227     return ret;
00228 }
00229 
00230 void Matrix::svd(Matrix& U, Matrix& V, Vector& lambda, unsigned int maxIterations, bool ignoreThreshold) const
00231 {
00232 	::svd(*this, U, V, lambda, maxIterations, ignoreThreshold);
00233 }
00234 
00235 void Matrix::eigensymm(Matrix& eigenvectors, Vector& eigenvalues) const
00236 {
00237     SIZE_CHECK(ndim() == 2);
00238     SIZE_CHECK(dim(0) == dim(1));
00239 	::eigensymm(*this, eigenvectors, eigenvalues);
00240 }
00241 
00242 void Matrix::eigenvalues(Vector& real, Vector* imaginary) const
00243 {
00244     SIZE_CHECK(ndim() == 2);
00245     SIZE_CHECK(dim(0) == dim(1));
00246     real.resize(dim(0), false);
00247     if (imaginary != NULL)
00248     {
00249         imaginary->resize(dim(0), false);
00250 		::eigen(*this, real, *imaginary);
00251     }
00252     else
00253     {
00254         Vector im(dim(0));
00255 		::eigen(*this, real, im);
00256     }
00257 }
00258 
00259 double Matrix::detSymm() const
00260 {
00261     SIZE_CHECK(ndim() == 2);
00262     SIZE_CHECK(dim(0) == dim(1));
00263     Matrix a(*this);
00264     Matrix v(dim(0), dim(1));
00265     Vector d(dim(0));
00266     return ::detsymm(a, v, d);
00267 }
00268 
00269 double Matrix::logDetSymm() const
00270 {
00271     SIZE_CHECK(ndim() == 2);
00272     SIZE_CHECK(dim(0) == dim(1));
00273     Matrix a(*this);
00274     Matrix v(dim(0), dim(1));
00275     Vector d(dim(0));
00276     return ::logdetsymm(a, v, d);
00277 }
00278 
00279 double Matrix::trace() const
00280 {
00281     SIZE_CHECK(isSquare());
00282 
00283     double ret = 0.0;
00284     unsigned int i, ic = dim(0);
00285     for (i=0; i<ic; i++) ret += operator () (i, i);
00286     return ret;
00287 }
00288 
00289 Vector Matrix::row(unsigned int r) const
00290 {
00291     SIZE_CHECK(ndim() == 2);
00292     RANGE_CHECK(r < dim(0));
00293 
00294     unsigned int i, ic = dim(1);
00295     Vector ret(ic);
00296     for (i=0; i<ic; i++) ret(i) = operator () (r, i);
00297     return ret;
00298 }
00299 
00300 Vector Matrix::col(unsigned int c) const
00301 {
00302     SIZE_CHECK(ndim() == 2);
00303     RANGE_CHECK(c < dim(1));
00304 
00305     unsigned int i, ic = dim(0);
00306     Vector ret(ic);
00307     for (i=0; i<ic; i++) ret(i) = operator () (i, c);
00308     return ret;
00309 }
00310 
00311 ArrayBase* Matrix::clone() const
00312 {
00313     return new Matrix(*this);
00314 }
00315 
00316 ArrayBase* Matrix::cloneEmpty() const
00317 {
00318     return new Matrix(0, 0);
00319 }
00320 
00321 void Matrix::resize_i(unsigned int* d, unsigned int nd, bool copy)
00322 {
00323     SIZE_CHECK(nd == 2);
00324     Array<double>::resize_i(d, 2, copy);
00325 }
00326 
00327 
00329 
00330 
00331 // comparison
00332 bool operator == (const Vector& v1, const Vector& v2)
00333 {
00334     return ((Array<double>&)(v1) == (Array<double>&)(v2)); 
00335 }
00336 
00337 bool operator != (const Vector& v1, const Vector& v2)
00338 {
00339     return ((Array<double>&)(v1) != (Array<double>&)(v2)); 
00340 }
00341 
00342 bool operator == (const Matrix& m1, const Matrix& m2)
00343 {
00344     return ((Array<double>&)(m1) == (Array<double>&)(m2)); 
00345 }
00346 
00347 bool operator != (const Matrix& m1, const Matrix& m2)
00348 {
00349     return ((Array<double>&)(m1) != (Array<double>&)(m2)); 
00350 }
00351 
00352 // addition and subtraction
00353 Vector operator + (const Vector& v1, const Vector& v2)
00354 {
00355     Vector ret(v1);
00356     ret += v2;
00357     return ret;
00358 }
00359 
00360 Vector operator - (const Vector& v1, const Vector& v2)
00361 {
00362     Vector ret(v1);
00363     ret -= v2;
00364     return ret;
00365 }
00366 
00367 Vector& operator += (Vector& v1, const Vector& v2)
00368 {
00369     unsigned int i, ic = v1.dim(0);
00370     SIZE_CHECK(v2.dim(0) == ic);
00371 
00372     for (i=0; i<ic; i++) v1(i) += v2(i);
00373 
00374     return v1;
00375 }
00376 
00377 Vector& operator -= (Vector& v1, const Vector& v2)
00378 {
00379     unsigned int i, ic = v1.dim(0);
00380     SIZE_CHECK(v2.dim(0) == ic);
00381 
00382     for (i=0; i<ic; i++) v1(i) -= v2(i);
00383 
00384     return v1;
00385 }
00386 
00387 Matrix operator + (const Matrix& m1, const Matrix& m2)
00388 {
00389     Matrix ret(m1);
00390     ret += m2;
00391     return ret;
00392 }
00393 
00394 Matrix operator - (const Matrix& m1, const Matrix& m2)
00395 {
00396     Matrix ret(m1);
00397     ret -= m2;
00398     return ret;
00399 }
00400 
00401 Matrix& operator += (Matrix& m1, const Matrix& m2)
00402 {
00403     unsigned int x, xc = m1.dim(1);
00404     unsigned int y, yc = m1.dim(0);
00405     SIZE_CHECK(m2.dim(0) == yc && m2.dim(1) == xc);
00406 
00407     for (y=0; y<yc; y++) for (x=0; x<xc; x++) m1(y, x) += m2(y, x);
00408 
00409     return m1;
00410 }
00411 
00412 Matrix& operator -= (Matrix& m1, const Matrix& m2)
00413 {
00414     unsigned int x, xc = m1.dim(1);
00415     unsigned int y, yc = m1.dim(0);
00416     SIZE_CHECK(m2.dim(0) == yc && m2.dim(1) == xc);
00417 
00418     for (y=0; y<yc; y++) for (x=0; x<xc; x++) m1(y, x) -= m2(y, x);
00419 
00420     return m1;
00421 }
00422 
00423 
00424 // scalar multiplication
00425 Vector operator * (double s, const Vector& v)
00426 {
00427     int y, yc = v.dim(0);
00428     Vector ret(yc);
00429     for (y=0; y<yc; y++) ret(y) = s * v(y);
00430     return ret;
00431 }
00432 
00433 Matrix operator * (double s, const Matrix& m)
00434 {
00435     int x, xc = m.dim(1);
00436     int y, yc = m.dim(0);
00437     Matrix ret(yc, xc);
00438     for (y=0; y<yc; y++) for (x=0; x<xc; x++) ret(y, x) = s * m(y, x);
00439     return ret;
00440 }
00441 
00442 Vector& operator *= (Vector& v, double s)
00443 {
00444     int y, yc = v.dim(0);
00445     for (y=0; y<yc; y++) v(y) *= s;
00446     return v;
00447 }
00448 
00449 Matrix& operator *= (Matrix& m, double s)
00450 {
00451     int x, xc = m.dim(1);
00452     int y, yc = m.dim(0);
00453     Matrix ret(yc, xc);
00454     for (y=0; y<yc; y++) for (x=0; x<xc; x++) m(y, x) *= s;
00455     return m;
00456 }
00457 
00458 
00459 // standard multiplication
00460 Vector operator * (const Matrix& m, const Vector& v)
00461 {
00462     unsigned int x, xc = m.dim(1);
00463     unsigned int y, yc = m.dim(0);
00464     SIZE_CHECK(v.dim(0) == xc);
00465     Vector ret(yc);
00466     for (y=0; y<yc; y++)
00467     {
00468         double a = 0.0;
00469         for (x=0; x<xc; x++) a += m(y, x) * v(x);
00470         ret(y) = a;
00471     }
00472     return ret;
00473 }
00474 
00475 Vector operator * (const Vector& v, const Matrix& m)
00476 {
00477     unsigned int x, xc = m.dim(1);
00478     unsigned int y, yc = m.dim(0);
00479     SIZE_CHECK(v.dim(0) == yc);
00480     Vector ret(xc);
00481     for (x=0; x<xc; x++)
00482     {
00483         double a = 0.0;
00484         for (y=0; y<yc; y++) a += m(y, x) * v(y);
00485         ret(x) = a;
00486     }
00487     return ret;
00488 }
00489 
00490 Matrix operator * (const Matrix& m1, const Matrix& m2)
00491 {
00492     unsigned int x, xc = m2.dim(1);
00493     unsigned int i, ic = m2.dim(0);
00494     unsigned int y, yc = m1.dim(0);
00495     SIZE_CHECK(m1.dim(1) == ic);
00496     Matrix ret(yc, xc);
00497     for (y=0; y<yc; y++)
00498     {
00499         for (x=0; x<xc; x++)
00500         {
00501             double a = 0.0;
00502             for (i=0; i<ic; i++) a += m1(y, i) * m2(i, x);
00503             ret(y, x) = a;
00504         }
00505     }
00506     return ret;
00507 }
00508 
00509 Matrix& operator *= (Matrix& m1, const Matrix& m2)
00510 {
00511     Matrix tmp(m1.dim(0), m1.dim(1));
00512     tmp = m1 * m2;
00513     m1 = tmp;
00514     return m1;
00515 }
00516 
00517 // inner product
00518 double operator * (const Vector& v1, const Vector& v2)
00519 {
00520     unsigned int i, ic = v1.dim(0);
00521     SIZE_CHECK(v2.dim(0) == ic);
00522     double ret = 0.0;
00523     for (i=0; i<ic; i++) ret += v1(i) * v2(i);
00524     return ret;
00525 }
00526 
00527 // outer product
00528 Matrix operator % (const Vector& v1, const Vector& v2)
00529 {
00530     int y, yc = v1.dim(0);
00531     int x, xc = v2.dim(0);
00532     Matrix ret(yc, xc);
00533     for (y=0; y<yc; y++) for (x=0; x<xc; x++) ret(y, x) = v1(y) * v2(x);
00534     return ret;
00535 }
00536 
00537 
00538 // advanced math
00539 Matrix powerseries(const Matrix& m, const double* coeff, int ncoeff, double acc)
00540 {
00541     SIZE_CHECK(m.isSquare());
00542 
00543     int dim = m.dim(0);
00544     Matrix ret(dim, dim);
00545     Matrix power(Matrix::unitmatrix(dim));
00546 
00547     int i;
00548     for (i=0; i<ncoeff; i++)
00549     {
00550         if (coeff[i] != 0.0)
00551         {
00552             if (acc > 0.0)
00553             {
00554                 double m = 0.0, M = 0.0;
00555                 double mp, mr;
00556                 minmaxElement(power, m, M);
00557                 mp = M; if (m < 0.0 && -m > M) mp = -m;
00558                 minmaxElement(ret, m, M);
00559                 mr = M; if (m < 0.0 && -m > M) mr = -m;
00560                 if (coeff[i] * mp < acc * mr) break;
00561             }
00562             ret += coeff[i] * power;
00563         }
00564         power *= m;
00565     }
00566 
00567     return ret;
00568 }
00569 
00570 Matrix exp(const Matrix& m)
00571 {
00572     double coeff[60] = {
00573         1.0, 1.0, 0.5,
00574         0.1666666666666666667, 0.04166666666666666667, 0.008333333333333333333,
00575         0.001388888888888888889, 0.0001984126984126984127, 2.480158730158730159e-05,
00576         2.755731922398589065e-06, 2.755731922398589065e-07, 2.505210838544171878e-08,
00577         2.087675698786809898e-09, 1.60590438368216146e-10, 1.147074559772972471e-11,
00578         7.647163731819816476e-13, 4.779477332387385297e-14, 2.811457254345520763e-15,
00579         1.561920696858622646e-16, 8.220635246624329717e-18, 4.110317623312164858e-19,
00580         1.957294106339126123e-20, 8.896791392450573287e-22, 3.868170170630684038e-23,
00581         1.611737571096118349e-24, 6.446950284384473396e-26, 2.47959626322479746e-27,
00582         9.183689863795546148e-29, 3.27988923706983791e-30, 1.130996288644771693e-31,
00583         3.769987628815905644e-33, 1.21612504155351795e-34, 3.800390754854743593e-36,
00584         1.151633562077195028e-37, 3.387157535521161847e-39, 9.677592958631890992e-41,
00585         2.688220266286636387e-42, 7.265460179153071315e-44, 1.911963205040281925e-45,
00586         4.902469756513543398e-47, 1.225617439128385849e-48, 2.989310827142404511e-50,
00587         7.117406731291439311e-52, 1.655210867742195189e-53, 3.761842881232261792e-55,
00588         8.359650847182803983e-57, 1.817315401561479127e-58, 3.866628513960593887e-60,
00589         8.055476070751237264e-62, 1.643974708316579034e-63, 3.287949416633158067e-65,
00590         6.44695964045717268e-67, 1.239799930857148592e-68, 2.339245152560657722e-70,
00591         4.331935467704921706e-72, 7.876246304918039466e-74, 1.406472554449649905e-75,
00592         2.467495709560789306e-77, 4.254302947518602253e-79, 7.210682961895936021e-81,
00593     };
00594 
00595     return powerseries(m, coeff, 60, 1e-14);
00596 }