PLearn: TMat_maths_impl.h Source File

00001 // -*- C++ -*- 00002 00003 // PLearn (A C++ Machine Learning Library) 00004 // Copyright (C) 1998 Pascal Vincent 00005 // Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal 00006 // 00007 00008 // Redistribution and use in source and binary forms, with or without 00009 // modification, are permitted provided that the following conditions are met: 00010 // 00011 // 1. Redistributions of source code must retain the above copyright 00012 // notice, this list of conditions and the following disclaimer. 00013 // 00014 // 2. Redistributions in binary form must reproduce the above copyright 00015 // notice, this list of conditions and the following disclaimer in the 00016 // documentation and/or other materials provided with the distribution. 00017 // 00018 // 3. The name of the authors may not be used to endorse or promote 00019 // products derived from this software without specific prior written 00020 // permission. 00021 // 00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR 00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN 00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00032 // 00033 // This file is part of the PLearn library. For more information on the PLearn 00034 // library, go to the PLearn Web site at www.plearn.org 00035 00036 00037 00038 00039 /* ******************************************************* 00040 * $Id: TMat_maths_impl.h,v 1.59 2004/07/22 14:00:31 tihocan Exp $ 00041 * AUTHORS: Pascal Vincent & Yoshua Bengio & Rejean Ducharme 00042 * This file is part of the PLearn library. 00043 ******************************************************* */ 00044 00047 #ifndef TMat_maths_impl_H 00048 #define TMat_maths_impl_H 00049 00050 #include <algorithm> 00051 00052 namespace PLearn { 00053 using namespace std; 00054 00055 template <class T> 00056 T max(const TVec<T>& vec) 00057 { 00058 #ifdef BOUNDCHECK 00059 if(vec.length()==0) 00060 PLERROR("IN max(const NumericVec& vec) TVec has zero length()"); 00061 #endif 00062 T* pv = vec.data(); 00063 T maxval = *pv++; 00064 int n = vec.length(); 00065 while(--n) 00066 { 00067 if(*pv>maxval) 00068 maxval = *pv; 00069 ++pv; 00070 } 00071 return maxval; 00072 } 00073 00075 template <class T> 00076 void softmax(const TVec<T>& x, const TVec<T>& y) 00077 { 00078 int n = x.length(); 00079 if (n>0) 00080 { 00081 T* yp = y.data(); 00082 T* xp = x.data(); 00083 T maxx = max(x); 00084 real s = 0; 00085 for (int i=0;i<n;i++) 00086 s += (yp[i] = safeexp(xp[i]-maxx)); 00087 if (s == 0) PLERROR("trying to divide by 0 in softmax"); 00088 s = 1.0 / s; 00089 for (int i=0;i<n;i++) 00090 yp[i] *= s; 00091 } 00092 } 00093 00094 // returns y = log(sofmax(x)) 00095 template <class T> 00096 void log_softmax(const TVec<T> &x, TVec<T> &y) 00097 { 00098 if (x.length() > 0) { 00099 y << x; 00100 y -= max(x); 00101 y -= logadd(y); 00102 } 00103 } 00104 00106 template <class T> 00107 void exp(const TVec<T>& x, TVec<T>& y) 00108 { 00109 y.resize(x.length()); 00110 int n = x.length(); 00111 T* xp = x.data(); 00112 T* yp = y.data(); 00113 while(n--) 00114 *yp++ = exp(*xp++); 00115 } 00116 00118 template<class T> 00119 T sumsquare(const TVec<T>& x) 00120 { 00121 T* v = x.data(); 00122 T res = square(v[0]); 00123 int l = x.length(); 00124 for(int i=1; i<l; i++) 00125 res += square(v[i]); 00126 return res; 00127 } 00128 00130 template<class T> 00131 T sumabs(const TVec<T>& x) 00132 { 00133 T* v = x.data(); 00134 T res = (T)(fabs((real)v[0])); 00135 int l = x.length(); 00136 for(int i=1; i<l; i++) 00137 res += (T)(fabs((real)v[i])); 00138 return res; 00139 } 00140 00142 template<class T> 00143 void squareElements(const TVec<T>& x) 00144 { 00145 T* ptr = x.data(); 00146 int l = x.length(); 00147 while(l--) 00148 { 00149 *ptr *= *ptr; 00150 ++ptr; 00151 } 00152 } 00153 00155 template<class T> 00156 void squareElements(const TMat<T>& m) 00157 { 00158 if (m.size()==0) return; 00159 if(m.isCompact()) { 00160 typename TMat<T>::compact_iterator it = m.compact_begin(); 00161 typename TMat<T>::compact_iterator itend = m.compact_end(); 00162 for(; it != itend; ++it) 00163 *it = square(*it); 00164 } else { 00165 typename TMat<T>::iterator it = m.begin(); 00166 typename TMat<T>::iterator itend = m.end(); 00167 for(; it != itend; ++it) 00168 *it = square(*it); 00169 } 00170 } 00171 00173 template<class T> 00174 T sumsquare(const TMat<T>& m) 00175 { 00176 if (m.size()==0) return 0; 00177 if(m.isCompact()) 00178 { 00179 typename TMat<T>::compact_iterator it = m.compact_begin(); 00180 typename TMat<T>::compact_iterator itend = m.compact_end(); 00181 T res = square(*it); 00182 ++it; 00183 for(; it!=itend; ++it) 00184 res += square(*it); 00185 return res; 00186 } 00187 else 00188 { 00189 typename TMat<T>::iterator it = m.begin(); 00190 typename TMat<T>::iterator itend = m.end(); 00191 T res = square(*it); 00192 ++it; 00193 for(; it!=itend; ++it) 00194 res += square(*it); 00195 return res; 00196 } 00197 } 00198 00199 00201 template<class T> 00202 T sumabs(const TMat<T>& m) 00203 { 00204 if (m.size()==0) return 0; 00205 if(m.isCompact()) 00206 { 00207 typename TMat<T>::compact_iterator it = m.compact_begin(); 00208 typename TMat<T>::compact_iterator itend = m.compact_end(); 00209 T res = fabs(*it); 00210 ++it; 00211 for(; it!=itend; ++it) 00212 res += fabs(*it); 00213 return res; 00214 } 00215 else 00216 { 00217 typename TMat<T>::iterator it = m.begin(); 00218 typename TMat<T>::iterator itend = m.end(); 00219 T res = fabs(*it); 00220 ++it; 00221 for(; it!=itend; ++it) 00222 res += fabs(*it); 00223 return res; 00224 } 00225 } 00226 00227 00229 template <class T> 00230 inline void multiply(const TVec<T>& source1, T source2, TVec<T>& destination) 00231 { 00232 int n=source1.length(); 00233 if (n!=destination.length()) 00234 destination.resize(n); 00235 T* s1=source1.data(); 00236 T* d=destination.data(); 00237 for (int i=0;i<n;i++) 00238 d[i] = s1[i]*source2; 00239 } 00240 00241 00242 //------- These were previously in Vec_maths 00243 00244 template<class T> 00245 T sum(const TVec<T>& vec, bool ignore_missing=false) 00246 { 00247 double res = 0.0; 00248 T* v = vec.data(); 00249 for(int i=0; i<vec.length(); i++) 00250 { 00251 if (!is_missing(v[i])) res += v[i]; 00252 else if (!ignore_missing) return MISSING_VALUE; 00253 } 00254 return T(res); 00255 } 00256 00260 template<class T> 00261 T sum_of_log(const TVec<T>& vec) 00262 { 00263 double res = 0.0; 00264 T* v = vec.data(); 00265 for(int i=0; i<vec.length(); i++) 00266 res += log(v[i]); 00267 return T(res); 00268 } 00269 00270 template<class T> 00271 T product(const TVec<T>& vec) 00272 { 00273 double res = 1.0; 00274 T* v = vec.data(); 00275 for(int i=0; i<vec.length(); i++) 00276 res *= v[i]; 00277 return T(res); 00278 } 00279 00280 /* 00281 template<class T> 00282 T mean(const TVec<T>& vec) 00283 { 00284 #ifdef BOUNDCHECK 00285 if(vec.length()==0) 00286 PLERROR("IN T mean(const TVec<T>& vec) vec has zero length"); 00287 #endif 00288 double res = 0.0; 00289 T* v = vec.data(); 00290 for(int i=0; i<vec.length(); i++) 00291 res += v[i]; 00292 return T(res/vec.length()); 00293 } 00294 */ 00295 00300 template<class T> 00301 T mean(const TVec<T>& vec, bool ignore_missing=false) 00302 { 00303 #ifdef BOUNDCHECK 00304 if(vec.length()==0) 00305 PLERROR("IN T mean(const TVec<T>& vec) vec has zero length"); 00306 #endif 00307 double res = 0.0; 00308 int n = 0; 00309 T* v = vec.data(); 00310 for(int i=0; i<vec.length(); i++) 00311 { 00312 if (!is_missing(v[i])) 00313 { 00314 res += v[i]; 00315 n++; 00316 } 00317 else if (!ignore_missing) return MISSING_VALUE; 00318 } 00319 00320 if (n == 0) return MISSING_VALUE; 00321 return T(res/double(n)); 00322 } 00323 00324 template<class T> 00325 T harmonic_mean(const TVec<T>& vec, bool ignore_missing=false) 00326 { 00327 #ifdef BOUNDCHECK 00328 if(vec.length()==0) 00329 PLERROR("IN T mean(const TVec<T>& vec) vec has zero length"); 00330 #endif 00331 double res = 0.0; 00332 int n = 0; 00333 T* v = vec.data(); 00334 for(int i=0; i<vec.length(); i++) 00335 { 00336 if (!is_missing(v[i])) 00337 { 00338 res += 1.0/v[i]; 00339 n++; 00340 } 00341 else if (!ignore_missing) return MISSING_VALUE; 00342 } 00343 00344 if (n == 0) return MISSING_VALUE; 00345 return T(double(n)/res); 00346 } 00347 00348 template<class T> 00349 T avgdev(const TVec<T>& vec, T meanval) 00350 { 00351 #ifdef BOUNDCHECK 00352 if(vec.length()==0) 00353 PLERROR("IN T avgdev(const TVec<T>& vec, T meanval) vec has zero length"); 00354 #endif 00355 double res = 0.0; 00356 T* v = vec.data(); 00357 for(int i=0; i<vec.length(); i++) 00358 res += fabs(v[i]-meanval); 00359 return res/vec.length(); 00360 } 00361 00362 template<class T> 00363 T geometric_mean(const TVec<T>& vec) 00364 { 00365 #ifdef BOUNDCHECK 00366 if(vec.length()==0) 00367 PLERROR("IN T geometric_mean(const TVec<T>& vec) vec has zero length"); 00368 #endif 00369 double res = 0.0; 00370 T* v = vec.data(); 00371 for(int i=0; i<vec.length(); i++) 00372 { 00373 T vi = v[i]; 00374 if (vi<=0) 00375 PLERROR("geometric_mean(TVec<T>): argument %g <=0 at position [%d]", 00376 vi,i); 00377 res += v[i]; 00378 } 00379 return T(exp(res/vec.length())); 00380 } 00381 00382 template<class T> 00383 T weighted_mean(const TVec<T>& vec, const TVec<T>& weights, bool ignore_missing=false) 00384 { 00385 #ifdef BOUNDCHECK 00386 if(vec.length()!=weights.length() || vec.length() == 0) 00387 PLERROR("IN T weighted_mean(const TVec<T>& vec, const TVec<T>& weights) vec and weights must have equal (non-zero) lengths"); 00388 #endif 00389 double res = 0.0; 00390 T sum_weights = 0.0; 00391 T* v = vec.data(); 00392 T* w = weights.data(); 00393 for(int i=0; i<vec.length(); i++) 00394 { 00395 if (!is_missing(v[i]) && !is_missing(w[i])) 00396 { 00397 res += v[i] * w[i]; 00398 sum_weights += w[i]; 00399 } 00400 else if (!ignore_missing) return MISSING_VALUE; 00401 } 00402 if (sum_weights == 0) 00403 PLERROR("IN T weighted_mean: sum(weights) == 0"); 00404 return T(res/sum_weights); 00405 } 00406 00407 template<class T> 00408 T variance(const TVec<T>& vec, T meanval) 00409 { 00410 #ifdef BOUNDCHECK 00411 if(vec.length()<=1) 00412 PLERROR("IN T variance(const TVec<T>& vec, T meanval) vec length must be more than one"); 00413 #endif 00414 double res = 0.0; 00415 T* v = vec.data(); 00416 for(int i=0; i<vec.length(); i++) 00417 { 00418 double diff = v[i]-meanval; 00419 res += diff*diff; 00420 } 00421 return res/(vec.length()-1); 00422 } 00423 00424 template<class T> 00425 T covariance(const TVec<T>& vec1, const TVec<T>& vec2, T mean1, T mean2) 00426 { 00427 #ifdef BOUNDCHECK 00428 if(vec1.length()<=1) 00429 PLERROR("IN T covariance(const TVec<T>& vec1, const TVec<T>& vec2, T mean1, T mean2) vec1's length must be more than one"); 00430 if(vec2.length()<=1) 00431 PLERROR("IN T covariance(const TVec<T>& vec1, const TVec<T>& vec2, T mean1, T mean2) vec2's length must be more than one"); 00432 #endif 00433 if(vec1.length() != vec2.length()) 00434 PLERROR("IN T covariance(const TVec<T>& vec1, const TVec<T>& vec2, T mean1, T mean2) the lengths of vec1 and vec2 must be same"); 00435 int length = vec1.length(); 00436 double res = 0.0; 00437 T* v1 = vec1.data(); 00438 T* v2 = vec2.data(); 00439 for(int i=0; i<length; i++) 00440 { 00441 double temp = (v1[i]-mean1)*(v2[i]-mean2); 00442 res += temp; 00443 } 00444 return res/(length - 1); 00445 } 00446 00447 template<class T> 00448 T weighted_variance(const TVec<T>& vec, const TVec<T>& weights, T no_weighted_mean, T weighted_mean) 00449 { 00450 #ifdef BOUNDCHECK 00451 if(vec.length()!=weights.length() || vec.length()==0) 00452 PLERROR("IN T weighted_variance(const TVec<T>& vec, const TVec<T>& weights, T no_weighted_mean, T weighted_mean) vec and weights must have equal (non-zero) lengths"); 00453 #endif 00454 double res = 0.0; 00455 T* v = vec.data(); 00456 T* w = weights.data(); 00457 for(int i=0; i<vec.length(); i++) 00458 res += v[i] * v[i] * w[i]; 00459 T sum_weights = sum(weights); 00460 if (sum_weights == 0) 00461 PLERROR("IN T weighted_variance(const TVec<T>& vec, const TVec<T>& weights, T no_weighted_mean, T weighted_mean) sum(weights) == 0"); 00462 return (res/sum_weights - no_weighted_mean * (2*weighted_mean - no_weighted_mean))*vec.length()/(vec.length()-1); 00463 } 00464 00465 template<class T> 00466 TVec<T> histogram(const TVec<T>& vec, T minval, T maxval, int nbins) 00467 { 00468 TVec<T> histo(nbins); 00469 T deltaval = maxval-minval + 1e-6; 00470 for(int i=0; i<vec.length(); i++) 00471 { 00472 T val = vec[i]; 00473 int binpos = int((val-minval)/deltaval*nbins); 00474 if(binpos>=0 && binpos<nbins) 00475 histo[binpos]++; 00476 } 00477 return histo; 00478 } 00479 00480 00481 template<class T> 00482 T min(const TVec<T>& vec) 00483 { 00484 #ifdef BOUNDCHECK 00485 if(vec.length()==0) 00486 PLERROR("IN T min(const TVec<T>& vec) vec has zero length"); 00487 #endif 00488 T* v = vec.data(); 00489 T minval = v[0]; 00490 for(int i=1; i<vec.length(); i++) 00491 if(v[i]<minval) 00492 minval = v[i]; 00493 return minval; 00494 } 00495 00496 template<class T> 00497 T maxabs(const TVec<T>& vec) 00498 { 00499 #ifdef BOUNDCHECK 00500 if(vec.length()==0) 00501 PLERROR("IN T maxabs(const TVec<T>& vec) vec has zero length"); 00502 #endif 00503 T* v = vec.data(); 00504 T maxval = fabs(v[0]); 00505 for(int i=1; i<vec.length(); i++) 00506 { 00507 T a=fabs(v[i]); 00508 if(a>maxval) 00509 maxval = a; 00510 } 00511 return maxval; 00512 } 00513 00514 // template<class T> 00515 // T minabs(const TVec<T>& vec) 00516 // { 00517 // #ifdef BOUNDCHECK 00518 // if(vec.length()==0) 00519 // PLERROR("IN T minabs(const TVec<T>& vec) vec has zero length"); 00520 // #endif 00521 // T* v = vec.data(); 00522 // T minval = fabs(v[0]); 00523 // for(int i=1; i<vec.length(); i++) 00524 // { 00525 // T a=fabs(v[i]); 00526 // if(a<minval) 00527 // minval = a; 00528 // } 00529 // return minval; 00530 // } 00531 00532 template<class T> 00533 T minabs(const TVec<T>& vec, int index = int()) 00534 { 00535 #ifdef BOUNDCHECK 00536 if(vec.length()==0) 00537 PLERROR("IN T minabs(const TVec<T>& vec) vec has zero length"); 00538 #endif 00539 T* v = vec.data(); 00540 T minval = fabs(v[0]); 00541 for(int i=1; i<vec.length(); i++) 00542 { 00543 T a=fabs(v[i]); 00544 if(a<minval) 00545 { 00546 index = i; 00547 minval = a; 00548 } 00549 } 00550 00551 return minval; 00552 } 00553 00554 00555 template<class T> 00556 int argmax(const TVec<T>& vec) 00557 { 00558 #ifdef BOUNDCHECK 00559 if(vec.length()==0) 00560 PLERROR("IN int argmax(const TVec<T>& vec) vec has zero length"); 00561 #endif 00562 T* v = vec.data(); 00563 int indexmax = 0; 00564 T maxval = v[0]; 00565 for(int i=1; i<vec.length(); i++) 00566 if(v[i]>maxval) 00567 { 00568 maxval = v[i]; 00569 indexmax = i; 00570 } 00571 return indexmax; 00572 } 00573 00574 template<class T> 00575 int argmax(const TVec<T>& vec, bool ignore_missing) 00576 { 00577 #ifdef BOUNDCHECK 00578 if(vec.length()==0) 00579 PLERROR("IN int argmax(const TVec<T>& vec) vec has zero length"); 00580 #endif 00581 T* v = vec.data(); 00582 int indexmax = -1; 00583 T maxval = MISSING_VALUE; 00584 00585 for(int i=0; i<vec.length(); i++) 00586 { 00587 if( is_missing(v[i]) ) 00588 { 00589 if(ignore_missing) continue; 00590 else PLERROR("argmax(const TVec<T>& vec, bool ignore_missing) encountered a MISSING_VALUE\n" 00591 "at index %d and ignore_missing is false.", i); 00592 } 00593 00594 if( indexmax == -1 || 00595 v[i] > maxval ) 00596 { 00597 maxval = v[i]; 00598 indexmax = i; 00599 } 00600 } 00601 return indexmax; 00602 } 00603 00604 00605 template<class T> 00606 int argmin(const TVec<T>& vec) 00607 { 00608 #ifdef BOUNDCHECK 00609 if(vec.length()==0) 00610 PLERROR("IN int argmin(const TVec<T>& vec) vec has zero length"); 00611 #endif 00612 T* v = vec.data(); 00613 int indexmin = 0; 00614 T minval = v[0]; 00615 for(int i=1; i<vec.length(); i++) 00616 if(v[i]<minval) 00617 { 00618 minval = v[i]; 00619 indexmin = i; 00620 } 00621 return indexmin; 00622 } 00623 00624 template<class T> 00625 int argmin(const TVec<T>& vec, bool ignore_missing) 00626 { 00627 #ifdef BOUNDCHECK 00628 if(vec.length()==0) 00629 PLERROR("IN int argmin(const TVec<T>& vec) vec has zero length"); 00630 #endif 00631 T* v = vec.data(); 00632 int indexmin = -1; 00633 T minval = MISSING_VALUE; 00634 00635 for(int i=0; i<vec.length(); i++) 00636 { 00637 if( is_missing(v[i]) ) 00638 { 00639 if(ignore_missing) continue; 00640 else PLERROR("argmin(const TVec<T>& vec, bool ignore_missing) encountered a MISSING_VALUE\n" 00641 "at index %d and ignore_missing is false.", i); 00642 } 00643 00644 if( indexmin == -1 || 00645 v[i] < minval ) 00646 { 00647 minval = v[i]; 00648 indexmin = i; 00649 } 00650 } 00651 return indexmin; 00652 } 00653 00654 00655 00656 template<class T> 00657 T pownorm(const TVec<T>& vec, double n) 00658 { 00659 double result = 0.0; 00660 T* v = vec.data(); 00661 if(n==1.0) 00662 { 00663 for(int i=0; i<vec.length(); i++) 00664 { 00665 T val = v[i]; 00666 if(val>=0) 00667 result += val; 00668 else 00669 result -= val; 00670 } 00671 } 00672 else if(n==2.0) 00673 { 00674 for(int i=0; i<vec.length(); i++) 00675 { 00676 T val = v[i]; 00677 result += val*val; 00678 } 00679 } 00680 else if(n==0) 00681 { result = vec.length(); } 00682 else 00683 { 00684 for(int i=0; i<vec.length(); i++) 00685 result += mypow(fabs(v[i]),n); 00686 } 00687 return result; 00688 } 00689 00690 template<class T> 00691 inline T pownorm(const TVec<T>& vec) { return pownorm(vec,T(2.0)); } 00692 00693 template<class T> 00694 T norm(const TVec<T>& vec, double n) 00695 { 00696 if(n==T(1.0)) 00697 return pownorm(vec, T(1.0)); 00698 else if(n==T(2.0)) 00699 return sqrt(pownorm(vec,T(2.0))); 00700 else 00701 return mypow(pownorm(vec,n), T(1.0)/n); 00702 } 00703 00704 template<class T> 00705 inline T norm(const TVec<T>& vec) { return norm(vec,T(2.0)); } 00706 00707 template<class T> 00708 void normalize(const TVec<T>& vec, double n) 00709 { vec /= norm(vec,n); } 00710 00711 template<class T> 00712 T powdistance(const TVec<T>& vec1, const TVec<T>& vec2, double n) 00713 { 00714 static T result, diff; 00715 #ifdef BOUNDCHECK 00716 if(vec1.length() != vec2.length()) 00717 PLERROR("In weighted_powdistance: vec1, vec2 should have the same length"); 00718 #endif 00719 result = 0.0; 00720 T* v1 = vec1.data(); 00721 T* v2 = vec2.data(); 00722 int length = vec1.length(); 00723 if(n==1.0) // L1 distance 00724 { 00725 for(int i=0; i<length; i++) 00726 { 00727 diff = *v1++ - *v2++; 00728 if(diff>=0) 00729 result += diff; 00730 else 00731 result -= diff; 00732 } 00733 } 00734 else if(n==2.0) 00735 { 00736 for(int i=0; i<length; i++) 00737 { 00738 diff = *v1++ - *v2++; 00739 result += diff*diff; 00740 } 00741 } 00742 else 00743 { 00744 for(int i=0; i<length; i++) 00745 { 00746 diff = *v1++ - *v2++; 00747 if(diff<0) 00748 diff = -diff; 00749 result += mypow(diff,n); 00750 } 00751 } 00752 return result; 00753 } 00754 00755 template<class T> 00756 inline T powdistance(const TVec<T>& vec1, const TVec<T>& vec2) 00757 { return powdistance(vec1, vec2, 2.0); } 00758 00759 template<class T> 00760 T dist(const TVec<T>& vec1, const TVec<T>& vec2, double n) 00761 { 00762 if(n==T(1.0)) 00763 return powdistance(vec1, vec2, T(1.0)); 00764 else if(n==T(2.0)) 00765 return sqrt(powdistance(vec1, vec2, T(2.0))); 00766 else 00767 return mypow(powdistance(vec1, vec2, n), T(1.0)/n); 00768 } 00769 00770 template<class T> 00771 inline T L2distance(const TVec<T>& vec1, const TVec<T>& vec2) 00772 { return dist(vec1, vec2, 2.0); } 00773 00774 template<class T> 00775 inline T L1distance(const TVec<T>& vec1, const TVec<T>& vec2) 00776 { return dist(vec1, vec2, 1.0); } 00777 00778 00779 template<class T> 00780 T weighted_powdistance(const TVec<T>& vec1, const TVec<T>& vec2, double n, const TVec<T>& weights) 00781 { 00782 #ifdef BOUNDCHECK 00783 if(vec1.length() != weights.length() || vec2.length()!=weights.length()) 00784 PLERROR("In weighted_powdistance: vec1, vec2 and weights vector should have the same length"); 00785 #endif 00786 T result = 0.0; 00787 T* v1 = vec1.data(); 00788 T* v2 = vec2.data(); 00789 T* w = weights.data(); 00790 int length = vec1.length(); 00791 if(n==1.0) // L1 distance 00792 { 00793 for(int i=0; i<length; i++) 00794 { 00795 T diff = w[i]*(v1[i]-v2[i]); 00796 if(diff>=0) 00797 result += diff; 00798 else 00799 result -= diff; 00800 } 00801 } 00802 else if(n==2.0) 00803 { 00804 for(int i=0; i<length; i++) 00805 { 00806 T diff = w[i]*(v1[i]-v2[i]); 00807 result += diff*diff; 00808 } 00809 } 00810 else 00811 { 00812 for(int i=0; i<length; i++) 00813 { 00814 T diff = w[i]*(v1[i]-v2[i]); 00815 if(diff<0) 00816 diff = -diff; 00817 result += mypow(diff,n); 00818 } 00819 } 00820 return result; 00821 } 00822 00823 template<class T> 00824 T weighted_distance(const TVec<T>& vec1, const TVec<T>& vec2, double n, const TVec<T>& weights) 00825 { 00826 if(n==1.0) 00827 return weighted_powdistance(vec1, vec2, 1.0, weights); 00828 else if(n==2.0) 00829 return sqrt(weighted_powdistance(vec1, vec2, 2.0, weights)); 00830 else 00831 return mypow(weighted_powdistance(vec1, vec2, n, weights), 1.0/n); 00832 } 00833 00834 00835 template<class T> 00836 void operator-=(const TVec<T>& vec1, const TVec<T>& vec2) 00837 { 00838 T* v1 = vec1.data(); 00839 T* v2 = vec2.data(); 00840 for(int i=0; i<vec1.length(); i++) 00841 v1[i] -= v2[i]; 00842 } 00843 00844 template<class T> 00845 void operator*=(const TVec<T>& vec1, const TVec<T>& vec2) 00846 { 00847 T* v1 = vec1.data(); 00848 T* v2 = vec2.data(); 00849 for(int i=0; i<vec1.length(); i++) 00850 v1[i] *= v2[i]; 00851 } 00852 00853 template<class T> 00854 inline void operator/=(const TVec<T>& vec, T scalar) 00855 { vec *= T(1.0)/scalar; } 00856 00857 template<class T> 00858 inline void operator/=(const TVec<T>& vec, int scalar) 00859 { vec /= T(scalar); } 00860 00861 template<class T> 00862 void compute_log(const TVec<T>& src, const TVec<T>& dest) 00863 { 00864 #ifdef BOUNDCHECK 00865 if(src.length()!=dest.length()) 00866 PLERROR("In log, src and dest vectors must have the same length"); 00867 #endif 00868 T* ps = src.data(); 00869 T* pd = dest.data(); 00870 int n = src.length(); 00871 for(int i=0; i<n; i++) 00872 *pd++ = log(*ps++); 00873 } 00874 00875 template<class T> 00876 inline TVec<T> log(const TVec<T>& src) 00877 { TVec<T> dest(src.length()); compute_log(src,dest); return dest; } 00878 00879 template<class T> 00880 void compute_sqrt(const TVec<T>& src, const TVec<T>& dest) 00881 { 00882 #ifdef BOUNDCHECK 00883 if(src.length()!=dest.length()) 00884 PLERROR("In sqrt, src and dest vectors must have the same length"); 00885 #endif 00886 T* ps = src.data(); 00887 T* pd = dest.data(); 00888 int n = src.length(); 00889 for(int i=0; i<n; i++) 00890 *pd++ = sqrt(*ps++); 00891 } 00892 00893 template<class T> 00894 inline TVec<T> sqrt(const TVec<T>& src) 00895 { TVec<T> dest(src.length()); compute_sqrt(src,dest); return dest; } 00896 00897 template<class T> 00898 void compute_safelog(const TVec<T>& src, const TVec<T>& dest) 00899 { 00900 #ifdef BOUNDCHECK 00901 if(src.length()!=dest.length()) 00902 PLERROR("In safelog, src and dest vectors must have the same length"); 00903 #endif 00904 T* ps = src.data(); 00905 T* pd = dest.data(); 00906 int n = src.length(); 00907 for(int i=0; i<n; i++) 00908 *pd++ = safelog(*ps++); 00909 } 00910 00911 template<class T> 00912 inline TVec<T> safelog(const TVec<T>& src) 00913 { TVec<T> dest(src.length()); compute_safelog(src,dest); return dest; } 00914 00915 template<class T> 00916 void operator/=(const TVec<T>& vec1, const TVec<T>& vec2) 00917 { 00918 T* v1 = vec1.data(); 00919 T* v2 = vec2.data(); 00920 for(int i=0; i<vec1.length(); i++) 00921 v1[i] /= v2[i]; 00922 } 00923 00924 template<class T> 00925 void compute_tanh(const TVec<T>& src, const TVec<T>& dest) 00926 { 00927 #ifdef BOUNDCHECK 00928 if(src.length()!=dest.length()) 00929 PLERROR("In tanh, src and dest vectors must have the same length"); 00930 #endif 00931 T* ps = src.data(); 00932 T* pd = dest.data(); 00933 int n = src.length(); 00934 for(int i=0; i<n; i++) 00935 *pd++ = tanh(*ps++); 00936 } 00937 00938 template<class T> 00939 void bprop_tanh(const TVec<T>& tanh_x, const TVec<T>& d_tanh_x, TVec<T>& d_x) 00940 { 00941 #ifdef BOUNDCHECK 00942 if(tanh_x.length()!=d_tanh_x.length()) 00943 PLERROR("In bprop_tanh, src and dest vectors must have the same length"); 00944 #endif 00945 int n = tanh_x.length(); 00946 if (n != d_x.length()) d_x.resize(n); 00947 T* y = tanh_x.data(); 00948 T* dy = d_tanh_x.data(); 00949 T* dx = d_x.data(); 00950 for(int i=0; i<n; i++) 00951 { 00952 real yi = *y++; 00953 *dx++ = *dy++ * (1 - yi*yi); 00954 } 00955 } 00956 00957 template<class T> 00958 inline TVec<T> tanh(const TVec<T>& src) 00959 { TVec<T> dest(src.length()); compute_tanh(src,dest); return dest; } 00960 00961 00962 template<class T> 00963 void compute_fasttanh(const TVec<T>& src, const TVec<T>& dest) 00964 { 00965 #ifdef BOUNDCHECK 00966 if(src.length()!=dest.length()) 00967 PLERROR("In fasttanh, src and dest vectors must have the same length"); 00968 #endif 00969 T* ps = src.data(); 00970 T* pd = dest.data(); 00971 int n = src.length(); 00972 for(int i=0; i<n; i++) 00973 *pd++ = fasttanh(*ps++); 00974 } 00975 00976 template<class T> 00977 inline TVec<T> fasttanh(const TVec<T>& src) 00978 { TVec<T> dest(src.length()); compute_fasttanh(src,dest); return dest; } 00979 00980 template<class T> 00981 void compute_sigmoid(const TVec<T>& src, const TVec<T>& dest) 00982 { 00983 #ifdef BOUNDCHECK 00984 if(src.length()!=dest.length()) 00985 PLERROR("In sigmoid, src and dest vectors must have the same length"); 00986 #endif 00987 T* ps = src.data(); 00988 T* pd = dest.data(); 00989 int n = src.length(); 00990 for(int i=0; i<n; i++) 00991 *pd++ = sigmoid(*ps++); 00992 } 00993 00994 template<class T> 00995 inline TVec<T> sigmoid(const TVec<T>& src) 00996 { TVec<T> dest(src.length()); compute_sigmoid(src,dest); return dest; } 00997 00998 00999 template<class T> 01000 void compute_fastsigmoid(const TVec<T>& src, const TVec<T>& dest) 01001 { 01002 #ifdef BOUNDCHECK 01003 if(src.length()!=dest.length()) 01004 PLERROR("In fastsigmoid, src and dest vectors must have the same length"); 01005 #endif 01006 T* ps = src.data(); 01007 T* pd = dest.data(); 01008 int n = src.length(); 01009 for(int i=0; i<n; i++) 01010 *pd++ = fastsigmoid(*ps++); 01011 } 01012 01013 template<class T> 01014 inline TVec<T> fastsigmoid(const TVec<T>& src) 01015 { TVec<T> dest(src.length()); compute_fastsigmoid(src,dest); return dest; } 01016 01017 template<class T> 01018 void compute_inverse_sigmoid(const TVec<T>& src, const TVec<T>& dest) 01019 { 01020 #ifdef BOUNDCHECK 01021 if(src.length()!=dest.length()) 01022 PLERROR("In inverse_sigmoid, src and dest vectors must have the same length"); 01023 #endif 01024 T* ps = src.data(); 01025 T* pd = dest.data(); 01026 int n = src.length(); 01027 for(int i=0; i<n; i++) 01028 *pd++ = inverse_sigmoid(*ps++); 01029 } 01030 01031 template<class T> 01032 inline TVec<T> inverse_sigmoid(const TVec<T>& src) 01033 { TVec<T> dest(src.length()); compute_inverse_sigmoid(src,dest); return dest; } 01034 01035 01036 template<class T> 01037 void negateElements(const TVec<T>& vec) 01038 { 01039 T* v = vec.data(); 01040 for(int i=0; i<vec.length(); i++) 01041 v[i] = -v[i]; 01042 } 01043 01044 template<class T> 01045 void invertElements(const TVec<T>& vec) 01046 { 01047 T* v = vec.data(); 01048 for(int i=0; i<vec.length(); i++) 01049 v[i] = 1.0/v[i]; 01050 } 01051 01052 template<class T> 01053 TVec<T> inverted(const TVec<T>& vec) 01054 { 01055 TVec<T> ret(vec.length()); 01056 T* v = vec.data(); 01057 for(int i=0; i<vec.length(); i++) 01058 ret[i] = 1.0/v[i]; 01059 return ret; 01060 } 01061 01062 01063 template<class T> 01064 void operator+=(const TVec<T>& vec, T scalar) 01065 { 01066 T* v = vec.data(); 01067 for(int i=0; i<vec.length(); i++) 01068 v[i] += scalar; 01069 } 01070 01071 template<class T> 01072 void operator-=(const TVec<T>& vec, T scalar) 01073 { vec += -scalar; } 01074 01075 template<class T> 01076 TVec<T> operator-(TVec<T> vec) 01077 { 01078 TVec<T> opposite(vec.length()); 01079 T *v=vec.data(); 01080 T *o=opposite.data(); 01081 for (int i=0;i<vec.length();i++) 01082 o[i] = - v[i]; 01083 return opposite; 01084 } 01085 01086 template<class T> 01087 T dot(const TVec<T>& vec1, const TVec<T>& vec2) 01088 { 01089 #ifdef BOUNDCHECK 01090 if(vec1.length()!=vec2.length()) 01091 PLERROR("In T operator*(const TVec<T>& vec1, const TVec<T>& vec2) (dot product) the 2 vecs must have the same length."); 01092 #endif 01093 01094 T res = 0; 01095 T* v1 = vec1.data(); 01096 T* v2 = vec2.data(); 01097 for(int i=0; i<vec1.length(); i++) 01098 res += v1[i]*v2[i]; 01099 return res; 01100 } 01101 01107 template<class V, class T, class U> 01108 V dot(const TVec<T>& vec1, const TVec& vec2) 01109 { 01110 #ifdef BOUNDCHECK 01111 if(vec1.length()!=vec2.length()) 01112 PLERROR("In T operator*(const TVec<T>& vec1, const TVec<T>& vec2) (dot product) the 2 vecs must have the same length."); 01113 #endif 01114 01115 V res = 0; 01116 T* v1 = vec1.data(); 01117 U* v2 = vec2.data(); 01118 for(int i=0; i<vec1.length(); i++) 01119 res += v1[i]*v2[i]; 01120 return res; 01121 } 01122 01123 template<class T> 01124 T dot(const TMat<T>& m1, const TMat<T>& m2) 01125 { 01126 #ifdef BOUNDCHECK 01127 if(m1.size()!=m2.size()) 01128 PLERROR("In T operator*(const TMat<T>& m1, const TVec<T>& vec2) (dot product) the 2 matrices must have the same number of elements."); 01129 #endif 01130 01131 T res = 0; 01132 T* v1 = m1.data(); 01133 T* v2 = m2.data(); 01134 if (m1.isCompact() && m2.isCompact()) 01135 for(int i=0; i<m1.size(); i++) 01136 res += v1[i]*v2[i]; 01137 else 01138 { 01139 TMatElementIterator<T> p1 = m1.begin(); 01140 TMatElementIterator<T> p2 = m2.begin(); 01141 for (int i=0; i<m1.size(); i++,++p1,++p2) 01142 res += *p1 * *p2; 01143 } 01144 return res; 01145 } 01146 01147 template<class T> 01148 TVec<T> operator-(const TVec<T>& v1, const TVec<T>& v2) 01149 { 01150 if (v1.length() != v2.length()) 01151 PLERROR("TVec<T> - TVec<T>: different lengths %d and %d", 01152 v1.length(), v2.length()); 01153 TVec<T> v(v1.length()); 01154 v << v1; 01155 v-=v2; 01156 return v; 01157 } 01158 01159 template<class T> 01160 TVec<T> operator-(T v1, const TVec<T>& v2) 01161 { 01162 TVec<T> v(v2.length()); 01163 v = -v2; 01164 v += v1; 01165 return v; 01166 } 01167 01168 template<class T> 01169 TVec<T> operator-(const TVec<T>& v1, T v2) 01170 { 01171 TVec<T> v(v1.length()); 01172 substract(v1,v2,v); 01173 return v; 01174 } 01175 01176 template<class T> 01177 TVec<T> operator+(const TVec<T>& v1, const TVec<T>& v2) 01178 { 01179 if (v1.length() != v2.length()) 01180 PLERROR("TVec<T> + TVec<T>: different lengths %d and %d", 01181 v1.length(), v2.length()); 01182 TVec<T> v(v1.length()); 01183 v << v1; 01184 v+=v2; 01185 return v; 01186 } 01187 01188 template<class T> 01189 TVec<T> operator+(T v1, const TVec<T>& v2) 01190 { 01191 TVec<T> v(v2.length()); 01192 add(v2,v1,v); 01193 return v; 01194 } 01195 01196 template<class T> 01197 TVec<T> operator+(const TVec<T>& v1, T v2) 01198 { 01199 TVec<T> v(v1.length()); 01200 add(v1,v2,v); 01201 return v; 01202 } 01203 01204 template<class T> 01205 TVec<T> operator%(const TVec<T>& v1, const TVec<T>& v2) 01206 { 01207 if (v1.length() != v2.length()) 01208 PLERROR("TVec<T> + TVec<T>: different lengths %d and %d", 01209 v1.length(), v2.length()); 01210 TVec<T> v(v1.length()); 01211 v << v1; 01212 v*=v2; 01213 return v; 01214 } 01215 01216 template<class T> 01217 TVec<T> operator*(T scalar, const TVec<T>& v) 01218 { 01219 TVec<T> result(v.length()); 01220 multiply(v,scalar,result); 01221 return result; 01222 } 01223 01224 template<class T> 01225 TVec<T> operator*(const TVec<T>& v1, T v2) 01226 { 01227 TVec<T> v(v1.length()); 01228 multiply(v1,v2,v); 01229 return v; 01230 } 01231 01232 template<class T> 01233 TVec<T> operator/(const TVec<T>& v1, const TVec<T>& v2) 01234 { 01235 if (v1.length() != v2.length()) 01236 PLERROR("TVec<T> + TVec<T>: different lengths %d and %d", 01237 v1.length(), v2.length()); 01238 TVec<T> v(v1.length()); 01239 v << v1; 01240 v/=v2; 01241 return v; 01242 } 01243 01244 template<class T> 01245 TVec<T> operator/(T v1, const TVec<T>& v2) 01246 { 01247 int n=v2.length(); 01248 TVec<T> v(n); 01249 T* s2=v2.data(); 01250 T* d=v.data(); 01251 for (int i=0;i<n;i++) 01252 d[i] = v1/s2[i]; 01253 return v; 01254 } 01255 01256 // norman: changed to unharmful declaration (see below old style) 01257 template<class T1, class T2> 01258 TVec<T1> operator/(const TVec<T1>& v1, T2 scalar) 01259 { 01260 TVec<T1> v(v1.length()); 01261 multiply(v1,T1(1.0)/(T1)scalar,v); 01262 return v; 01263 } 01264 01265 // norman: harmful declarations 01266 // Replaced with a better declaration above 01267 //template<class T> 01268 //TVec<T> operator/(const TVec<T>& v1, T scalar) 01269 //{ 01270 // TVec<T> v(v1.length()); 01271 // multiply(v1,T(1.0)/scalar,v); 01272 // return v; 01273 //} 01274 01275 // norman: This will cause problems if T = int (recursive declaration) 01276 // Replaced with a better declaration above 01277 //template<class T> 01278 //TVec<T> operator/(const TVec<T>& v1, int scalar) 01279 //{ return v1/T(scalar); } 01280 01281 template<class T> 01282 T logadd(const TVec<T>& vec) 01283 { 01284 int l = vec.length(); 01285 if(l==0) 01286 return LOG_INIT; 01287 01288 T *p_x = vec.data(); 01289 T sum = *p_x++; 01290 for (int i=1; i<l; i++, p_x++) 01291 sum = logadd(sum, *p_x); 01292 return sum; 01293 } 01294 01295 template<class T> 01296 T output_margin(const TVec<T>& class_scores, int correct_class) 01297 { 01298 T maxother = -FLT_MAX; 01299 for(int i=0; i<class_scores.length(); i++) 01300 { 01301 if(i!=correct_class && class_scores[i]>maxother) 01302 maxother = class_scores[i]; 01303 } 01304 return class_scores[correct_class]-maxother; 01305 } 01306 01307 template<class T> 01308 void fill_one_hot(const TVec<T>& vec, int hotpos, T coldvalue, T hotvalue) 01309 { 01310 #ifdef BOUNDCHECK 01311 if(!vec) 01312 PLERROR("In fill_one_hot given vec must have the correct size"); 01313 if(hotpos<0 || (vec.length()==1 && hotpos>1) || (vec.length()>1 && hotpos>=vec.length())) 01314 PLERROR("In fill_one_hot given hotpos out of vec range"); 01315 #endif 01316 if(vec.length()==1) 01317 vec[0] = (hotpos==0 ?coldvalue :hotvalue); 01318 else 01319 { 01320 vec.fill(coldvalue); 01321 vec[hotpos] = hotvalue; 01322 } 01323 } 01324 01325 template<class T> 01326 TVec<T> one_hot(int length, int hotpos, T coldvalue, T hotvalue) 01327 { 01328 TVec<T> result(length); 01329 fill_one_hot(result, hotpos, coldvalue, hotvalue); 01330 return result; 01331 } 01332 01333 template<class T> 01334 TVec<T> square(const TVec<T>& vec) 01335 { 01336 int n = vec.length(); 01337 TVec<T> result(n); 01338 T* v = vec.data(); 01339 T* r = result.data(); 01340 for(int i=0; i<n; i++) 01341 r[i] = v[i]*v[i]; 01342 return result; 01343 } 01344 01345 template<class T> 01346 TVec<T> squareroot(const TVec<T>& vec) 01347 { 01348 int n = vec.length(); 01349 TVec<T> result(n); 01350 T* v = vec.data(); 01351 T* r = result.data(); 01352 for(int i=0; i<n; i++) 01353 r[i] = sqrt(v[i]); 01354 return result; 01355 } 01356 01357 template<class T> 01358 TVec<T> remove_missing(const TVec<T>& vec) 01359 { 01360 int n = vec.length(); 01361 int n_non_missing = 0; 01362 TVec<T> result(n); 01363 T* v = vec.data(); 01364 T* r = result.data(); 01365 for(int i=0; i<n; i++) { 01366 if (!is_missing(v[i])) 01367 r[n_non_missing++] = v[i]; 01368 } 01369 result.resize(n_non_missing); 01370 return result; 01371 } 01372 01375 template<class T, class U, class V> 01376 TVec apply(const TVec<T>& vec, U (*func)(V)) 01377 { 01378 TVec destination(vec.length()); 01379 apply(vec,destination,func); 01380 return destination; 01381 } 01382 01384 template<class T, class U> 01385 void apply(const TVec<T>& source, TVec& destination, U (*func)(T)) 01386 { 01387 int n=source.length(); 01388 if (n!=destination.length()) 01389 PLERROR("apply: source(%d) and destination(%d) TVec<T>'s must have same length", 01390 n,destination.length()); 01391 T* s = source.data(); 01392 U* d = destination.data(); 01393 for(int i=0; i<n; i++) 01394 d[i]=func(s[i]); 01395 } 01396 01399 template<class T, class U, class V> 01400 void apply(const TVec<T>& src1,const TVec& src2, TVec<V>& dest, 01401 V (*func)(T,U)) 01402 { 01403 int n=src1.length(); 01404 if (n!=dest.length() || n!=src2.length()) 01405 PLERROR("apply: src1, src2 and destination TVec<T>'s must have same length"); 01406 T* s1 = src1.data(); 01407 U* s2 = src2.data(); 01408 V* d = dest.data(); 01409 for(int i=0; i<n; i++) 01410 d[i]=func(s1[i],s2[i]); 01411 } 01412 01413 01414 // Efficient mathematical operations (without memory allocation) 01415 01416 // destination[i] = source1[i]*source2[i] 01417 template<class T> 01418 void multiply(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 01419 { 01420 int n=source1.length(); 01421 if (n!=source2.length()) 01422 PLERROR("multiply: two sources (l=%d and %d) must have same length", 01423 n,source2.length()); 01424 if (n!=destination.length()) 01425 destination.resize(n); 01426 T* s1=source1.data(); 01427 T* s2=source2.data(); 01428 T* d=destination.data(); 01429 for (int i=0;i<n;i++) 01430 d[i] = s1[i]*s2[i]; 01431 } 01432 01433 // destination[i] = source1[i] + source2[i]*source3 01434 template<class T> 01435 void multiplyAdd(const TVec<T>& source1, const TVec<T>& source2, 01436 T source3, TVec<T>& destination) 01437 { 01438 int n=source1.length(); 01439 if (n!=source2.length()) 01440 PLERROR("multiply: two sources (l=%d and %d) must have same length", 01441 n,source2.length()); 01442 if (n!=destination.length()) 01443 destination.resize(n); 01444 T* s1=source1.data(); 01445 T* s2=source2.data(); 01446 T* d=destination.data(); 01447 for (int i=0;i<n;i++) 01448 d[i] = s1[i]+s2[i]*source3; 01449 } 01450 01451 // destination[i] = a*destination[i] + b*source[i] 01452 template<class T> 01453 void multiplyScaledAdd(const TVec<T>& source, T a, T b, TVec<T>& destination) 01454 { 01455 int n=source.length(); 01456 if (n!=destination.length()) 01457 PLERROR("multiply: source and destination (l=%d and %d) must have same length", 01458 n,destination.length()); 01459 T* s=source.data(); 01460 T* d=destination.data(); 01461 for (int i=0;i<n;i++) 01462 d[i] = a*d[i] + b*s[i]; 01463 } 01464 01465 // destination[i] = source1[i]+source2[i] 01466 template<class T> 01467 void add(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 01468 { 01469 int n=source1.length(); 01470 if (n!=source2.length()) 01471 PLERROR("add: two sources (l=%d and %d) must have same length", 01472 n,source2.length()); 01473 if (n!=destination.length()) 01474 destination.resize(n); 01475 T* s1=source1.data(); 01476 T* s2=source2.data(); 01477 T* d=destination.data(); 01478 for (int i=0;i<n;i++) 01479 d[i] = s1[i]+s2[i]; 01480 } 01481 01482 // destination[i] = source1[i]+source2 01483 template<class T> 01484 void add(const TVec<T>& source1, T source2, TVec<T>& destination) 01485 { 01486 int n=source1.length(); 01487 if (n!=destination.length()) 01488 destination.resize(n); 01489 T* s1=source1.data(); 01490 T* d=destination.data(); 01491 for (int i=0;i<n;i++) 01492 d[i] = s1[i]+source2; 01493 } 01494 01495 template<class T> 01496 inline void substract(const TVec<T>& source1, T source2, TVec<T>& destination) 01497 { add(source1,-source2,destination); } 01498 01499 // destination[i] = source1[i]-source2[i] 01500 template<class T> 01501 void substract(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 01502 { 01503 int n=source1.length(); 01504 if (n!=source2.length()) 01505 PLERROR("substract: two sources (l=%d and %d) must have same length", 01506 n,source2.length()); 01507 if (n!=destination.length()) 01508 destination.resize(n); 01509 T* s1=source1.data(); 01510 T* s2=source2.data(); 01511 T* d=destination.data(); 01512 for (int i=0;i<n;i++) 01513 d[i] = s1[i]-s2[i]; 01514 } 01515 01516 template<class T> 01517 inline void divide(const TVec<T>& source1, T source2, TVec<T>& destination) 01518 { multiply(source1,1.0/source2,destination); } 01519 01520 // destination[i] = source1[i]/source2[i] 01521 template<class T> 01522 void divide(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 01523 { 01524 int n=source1.length(); 01525 if (n!=source2.length()) 01526 PLERROR("divide: two sources (l=%d and %d) must have same length", 01527 n,source2.length()); 01528 if (n!=destination.length()) 01529 destination.resize(n); 01530 T* s1=source1.data(); 01531 T* s2=source2.data(); 01532 T* d=destination.data(); 01533 for (int i=0;i<n;i++) 01534 d[i] = s1[i]/s2[i]; 01535 } 01536 01537 // destination[i] = source1/source2[i] 01538 template<class T> 01539 void divide(T source1, const TVec<T>& source2, TVec<T>& destination) 01540 { 01541 int n=source2.length(); 01542 if (n!=destination.length()) 01543 destination.resize(n); 01544 T* s2=source2.data(); 01545 T* d=destination.data(); 01546 for (int i=0;i<n;i++) 01547 d[i] = source1/s2[i]; 01548 } 01549 01550 // destination[i] = max(source1[i],source2[i]) 01551 template<class T> 01552 void max(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 01553 { 01554 int n=source1.length(); 01555 if (n!=source2.length()) 01556 PLERROR("max: two sources (l=%d and %d) must have same length", 01557 n,source2.length()); 01558 if (n!=destination.length()) 01559 destination.resize(n); 01560 T* s1=source1.data(); 01561 T* s2=source2.data(); 01562 T* d=destination.data(); 01563 for (int i=0;i<n;i++) 01564 d[i] = MAX(s1[i],s2[i]); 01565 } 01566 01567 // destination[i] = max(source1[i],source2) 01568 template<class T> 01569 void max(const TVec<T>& source1, T source2, TVec<T>& destination) 01570 { 01571 int n=source1.length(); 01572 if (n!=destination.length()) 01573 destination.resize(n); 01574 T* s1=source1.data(); 01575 T* d=destination.data(); 01576 for (int i=0;i<n;i++) 01577 d[i] = MAX(s1[i],source2); 01578 } 01579 01580 01581 // destination[i] = min(source1[i],source2[i]) 01582 template<class T> 01583 void min(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 01584 { 01585 int n=source1.length(); 01586 if (n!=source2.length()) 01587 PLERROR("min: two sources (l=%d and %d) must have same length", 01588 n,source2.length()); 01589 if (n!=destination.length()) 01590 destination.resize(n); 01591 T* s1=source1.data(); 01592 T* s2=source2.data(); 01593 T* d=destination.data(); 01594 for (int i=0;i<n;i++) 01595 d[i] = MIN(s1[i],s2[i]); 01596 } 01597 01598 // destination[i] = min(source1[i],source2) 01599 template<class T> 01600 void min(const TVec<T>& source1, T source2, TVec<T>& destination) 01601 { 01602 int n=source1.length(); 01603 if (n!=destination.length()) 01604 destination.resize(n); 01605 T* s1=source1.data(); 01606 T* d=destination.data(); 01607 for (int i=0;i<n;i++) 01608 d[i] = MIN(s1[i],source2); 01609 } 01610 01611 01612 template<class T> 01613 TVec<T> softmax(const TVec<T>& x) 01614 { 01615 TVec<T> y(x.length()); 01616 softmax(x,y); 01617 return y; 01618 } 01619 01620 template<class T> 01621 void tanh(const TVec<T>& x, TVec<T>& y) 01622 { 01623 int n = x.length(); 01624 #ifdef BOUNDCHECK 01625 if (y.length()!=n) 01626 PLERROR("tanh(TVec<T>,TVec<T>), second argument of length %d, first of length %d, should be =", 01627 n,y.length()); 01628 #endif 01629 if (n>0) 01630 { 01631 T* yp = y.data(); 01632 T* xp = x.data(); 01633 for (int i=0;i<n;i++) 01634 yp[i] = tanh(xp[i]); 01635 } 01636 } 01637 01638 template<class T> 01639 TVec<T> exp(TVec<T> vec) { return apply(vec,safeexp); } 01640 01641 // return indices of non-zero elements 01642 template<class T> 01643 TVec<T> nonZeroIndices(TVec<T> v) 01644 { 01645 int n=v.length(); 01646 TVec<T> indices(n); 01647 int ni=0; 01648 T* val = v.data(); 01649 T* indx= indices.data(); 01650 for (int i=0;i<n;i++) 01651 if (val[i]!=0) 01652 indx[ni++]=i; 01653 indices.resize(ni); 01654 return indices; 01655 } 01656 01657 // return indices of non-zero elements 01658 template<class T> 01659 TVec<T> nonZeroIndices(TVec<bool> v) 01660 { 01661 int n=v.length(); 01662 TVec<T> indices(n); 01663 int ni=0; 01664 bool* val = v.data(); 01665 T* indx= indices.data(); 01666 for (int i=0;i<n;i++) 01667 if (val[i]) 01668 indx[ni++]=i; 01669 indices.resize(ni); 01670 return indices; 01671 } 01672 01673 // Set the complement indices, i.e. if 0<=i<n is not an element 01674 // of the indices vector it is put in the complement_indices vector. 01675 template<class T> 01676 void complement_indices(TVec<T>& indices, int n, 01677 TVec<T>& complement_indices, 01678 TVec<T>& buffer) 01679 { 01680 int ni=indices.length(); 01681 T* ind = indices.data(); 01682 T* cind = complement_indices.data(); 01683 buffer.resize(n); 01684 buffer.fill(0); 01685 T* buf=buffer.data(); 01686 for (int i=0;i<ni;i++) 01687 buf[(int)ind[i]]=1.0; 01688 for (int i=0,j=0;i<n;i++) 01689 if (buf[i]==0.0) 01690 cind[j++]=i; 01691 } 01692 01693 // dest[i] = 1 if src[i]==v, 0 otherwise 01694 template<class T> 01695 void equals(const TVec<T>& src, T v, TVec<T>& dest) 01696 { 01697 T* s=src.data(); 01698 T* d=dest.data(); 01699 int n=src.length(); 01700 #ifdef BOUNDCHECK 01701 if (n!=dest.length()) 01702 PLERROR("equals(TVec<T>(%d),T,TVec<T>(%d)) args of unequal lengths", 01703 n,dest.length()); 01704 #endif 01705 for (int i=0;i<n;i++) 01706 if (s[i]==v) d[i]=1.0; else d[i]=0.0; 01707 } 01708 01709 // dest[i] = 1 if first[i] > second[i], 0 otherwise 01710 template<class T> 01711 void isLargerThan(const TVec<T>& first, const TVec<T>& second, TVec<T>& dest) 01712 { 01713 T* f=first.data(); 01714 T* s=second.data(); 01715 T* d=dest.data(); 01716 int n=first.length(); 01717 if(n!=second.length() || n!=dest.length()) 01718 PLERROR("isLargerThan(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal length", 01719 n, second.length(), dest.length()); 01720 for (int i=0; i<n; i++) 01721 d[i] = f[i] > s[i]; 01722 } 01723 01724 // dest[i] = 1 if first[i] >= second[i], 0 otherwise 01725 template<class T> 01726 void isLargerThanOrEqualTo(const TVec<T>& first, const TVec<T>& second, TVec<T>& dest) 01727 { 01728 T* f=first.data(); 01729 T* s=second.data(); 01730 T* d=dest.data(); 01731 int n=first.length(); 01732 if(n!=second.length() || n!=dest.length()) 01733 PLERROR("isLargerThan(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal length", 01734 n, second.length(), dest.length()); 01735 for (int i=0; i<n; i++) 01736 d[i] = f[i] >= s[i]; 01737 } 01738 01739 // dest[i] = 1 if first[i] < second[i], 0 otherwise 01740 template<class T> 01741 void isSmallerThan(const TVec<T>& first, const TVec<T>& second, TVec<T>& dest) 01742 { 01743 T* f=first.data(); 01744 T* s=second.data(); 01745 T* d=dest.data(); 01746 int n=first.length(); 01747 if(n!=second.length() || n!=dest.length()) 01748 PLERROR("isLargerThan(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal length", 01749 n, second.length(), dest.length()); 01750 for (int i=0; i<n; i++) 01751 d[i] = f[i] < s[i]; 01752 } 01753 01754 // dest[i] = 1 if first[i] <= second[i], 0 otherwise 01755 template<class T> 01756 void isSmallerThanOrEqualTo(const TVec<T>& first, const TVec<T>& second, TVec<T>& dest) 01757 { 01758 T* f=first.data(); 01759 T* s=second.data(); 01760 T* d=dest.data(); 01761 int n=first.length(); 01762 if(n!=second.length() || n!=dest.length()) 01763 PLERROR("isLargerThan(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal length", 01764 n, second.length(), dest.length()); 01765 for (int i=0; i<n; i++) 01766 d[i] = f[i] <= s[i]; 01767 } 01768 01769 // dest[i] = if_vec[i] ? then_vec[i] : else_vec[i]; 01770 template<class T> 01771 void ifThenElse(const TVec<T>& if_vec, const TVec<T>& then_vec, const TVec<T>& else_vec, TVec<T>& dest) 01772 { 01773 T* i_=if_vec.data(); 01774 T* t_=then_vec.data(); 01775 T* e_=else_vec.data(); 01776 T* d_=dest.data(); 01777 int n=if_vec.length(); 01778 if (n!=then_vec.length() || n!=else_vec.length()) 01779 PLERROR("ifThenElse(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal lengths", 01780 n, then_vec.length(), else_vec.length(), dest.length()); 01781 for (int i=0;i<n;i++) 01782 d_[i] = i_[i] ? t_[i] : e_[i]; 01783 } 01784 01785 // returns the number of times that src[i] == value 01786 template<class T> 01787 int vec_counts(const TVec<T>& src, T value) 01788 { 01789 int len = src.length(); 01790 int n = 0; 01791 T *p = src.data(); 01792 for (int i=0; i<len; i++, p++) 01793 if (*p == value) 01794 n++; 01795 return n; 01796 } 01797 01798 // returns the position of f in src (-1 if f is not found) 01799 template<class T> 01800 int vec_find(const TVec<T>& src, T f) 01801 { 01802 int len = src.length(); 01803 T *p = src.data(); 01804 for (int i=0; i<len; i++, p++) 01805 if (*p == f) 01806 return(i); 01807 return -1; 01808 } 01809 01810 01811 template<class T> 01812 T estimatedCumProb(T x, TVec<T> bins) 01813 { 01814 const int nbins = bins.length()-1; 01815 if (nbins<1) PLERROR("estimatedCumProb:: there should be at least two elements in the bins vector"); 01816 // +0.5 because we allocate mass 0.25 at the left and 0.25 at the right of the interval (bins(0),bins(nbins)) 01817 const T one_over_nbins = 1.0/(T)(nbins+0.5); 01818 01819 int k = binary_search(bins, x); 01820 01821 if (k == -1) 01822 return 0.25*one_over_nbins; 01823 else if (k == nbins-1) 01824 return 1.0 - 0.25*one_over_nbins; 01825 else if (bins[k] != bins[k+1]) 01826 return one_over_nbins*(0.25 + k + (x-bins[k])/(bins[k+1]-bins[k])); 01827 else 01828 return one_over_nbins*(0.75 + k); 01829 } 01830 01831 // returns the index of the kth ordered element of v 01832 // (dumb algorithm, takes time in k*n ) 01833 template<class T> 01834 int positionOfkthOrderedElement(const TVec<T>& vec, int k) 01835 { 01836 #ifdef BOUNDCHECK 01837 if(k<0 || k>=vec.length()) 01838 PLERROR("In positionOfkthOrderedElement, k out of bounds"); 01839 #endif 01840 01841 T* v = vec.data(); 01842 01843 T minval = -FLT_MAX; 01844 int pos = -1; 01845 int l=0; 01846 01847 while(l<=k) 01848 { 01849 int nelements_equal_to_newminval = 0; 01850 T newminval = FLT_MAX; 01851 for(int i=0; i<vec.length(); i++) 01852 { 01853 if(v[i]>minval) 01854 { 01855 if(v[i]<newminval) 01856 { 01857 newminval = v[i]; 01858 nelements_equal_to_newminval = 1; 01859 pos = i; 01860 } 01861 else if(v[i]==newminval) 01862 nelements_equal_to_newminval++; 01863 } 01864 } 01865 l += nelements_equal_to_newminval; 01866 minval = newminval; 01867 } 01868 01869 return pos; 01870 } 01871 01874 template<class T> 01875 inline T kthOrderedElement(const TVec<T>& vec, int k) 01876 { return vec[positionOfkthOrderedElement(vec,k)]; } 01877 01879 template<class T> 01880 inline T median(const TVec<T>& vec) 01881 { return kthOrderedElement(vec, (vec.length()-1)/2); } 01882 01883 01884 //-------------- These were previouslty methods of TVec ---------------------------------- 01885 01886 01889 template<class T> 01890 T selectAndOrder(const TVec<T>& vec, int pos) 01891 { 01892 if (pos<0 || pos>=vec.length()) PLERROR("Bad position (%d)", pos); 01893 01894 int l=0; 01895 int h=vec.length()-1; 01896 T* v = vec.data(); 01897 01898 while (l<h) 01899 { 01900 T p = v[(l+h)/2]; 01901 int x = l; 01902 int y = h; 01903 01904 do 01905 { 01906 while (v[x]<p) x++; 01907 while (p<v[y]) y--; 01908 if (x<=y) 01909 { 01910 PLearn::swap(v[x],v[y]); 01911 x++; 01912 y--; 01913 } 01914 } while (x<=y); 01915 01916 if (pos>=x) l=x; 01917 else h=x-1; 01918 } 01919 01920 return v[l]; 01921 } 01922 01927 template<class T> 01928 TVec<T> getQuantiles(const TVec<T>& vec, int q) 01929 { 01930 int l = vec.length(); 01931 T* v = vec.data(); 01932 TVec<T> w(q+1); 01933 T linvq = T(l)/q; 01934 for(int i=0;i<q;i++) w[i] = v[int(linvq*i)]; 01935 w[q]=v[l-1]; 01936 return w; 01937 } 01938 01941 template<class T> 01942 TVec<T> nonZero(const TVec<T>& vec) 01943 { 01944 T *v =vec.data(); 01945 int n=0; 01946 for(int i=0;i<vec.length(); i++) if (v[i]!=0) n++; 01947 TVec<T> w(n); 01948 int j=0; 01949 for(int i=0;i<vec.length(); i++) if (v[i]!=0) w[j++]=v[i]; 01950 return(w); 01951 } 01952 01955 template<class T> 01956 TVec<T> positiveValues(const TVec<T>& vec) 01957 { 01958 T *v =vec.data(); 01959 int n=0; 01960 for(int i=0;i<vec.length(); i++) if (v[i]>0) n++; 01961 TVec<T> w(n); 01962 int j=0; 01963 for(int i=0;i<vec.length(); i++) if (v[i]>0) w[j++]=v[i]; 01964 return(w); 01965 } 01966 01971 template<class T> 01972 int positionOfClosestElement(const TVec<T>& vec, const T& value, bool is_sorted_vec=false) 01973 { 01974 T* v = vec.data(); 01975 if (is_sorted_vec) // dichotomy search 01976 { 01977 int pos = binary_search(vec, value); 01978 if (pos == -1) return 0; 01979 else if (pos == vec.length()-1) return pos; 01980 T dist1 = fabs(v[pos]-value); 01981 T dist2 = fabs(v[pos+1]-value); 01982 if (dist1 <= dist2) return pos; 01983 else return pos+1; 01984 } 01985 else // linear search 01986 { 01987 int pos_of_closest = 0; 01988 T dist_to_closest = fabs(v[0]-value); 01989 for(int i=1; i<vec.length(); i++) 01990 { 01991 T dist = fabs(v[i]-value); 01992 if(dist<dist_to_closest) 01993 { 01994 pos_of_closest = i; 01995 dist_to_closest = dist; 01996 } 01997 } 01998 return pos_of_closest; 01999 } 02000 } 02001 02002 02009 template <class T> 02010 void projectOnOrthogonalSubspace(const TVec<T>& vec, const TMat<T>& orthonormal_subspace) 02011 { 02012 for (int i=0;i<orthonormal_subspace.length();i++) 02013 { 02014 TVec<T> vi = orthonormal_subspace(i); 02015 T dp = dot(vec,vi); 02016 multiplyAcc(vec, vi,-dp); 02017 } 02018 } 02019 02020 template<class T> 02021 void operator*=(const TVec<T>& vec, T factor) 02022 { 02023 T* p = vec.data(); 02024 int l = vec.length(); 02025 for (int i=0;i<l;i++) 02026 *p++ *= factor; 02027 } 02028 02030 template<class T> 02031 inline void operator+=(const TVec<T>& vec1, const TVec<T>& vec2) 02032 { 02033 T* v1 = vec1.data(); 02034 T* v2 = vec2.data(); 02035 int l = vec1.length(); 02036 for(int i=0; i<l; i++) 02037 *v1++ += *v2++; 02038 } 02039 02040 02042 template<class T> 02043 void multiplyAcc(const TVec<T>& vec, const TVec<T>& x, T scale) 02044 { 02045 int n=x.length(); 02046 if (vec.length()!=n) 02047 PLERROR("TVec::multiplyAcc this has length_=%d and x has length_=%d", vec.length(),n); 02048 T* p=vec.data(); 02049 T* xp=x.data(); 02050 for (int i=0;i<n;i++) 02051 *p++ += scale * *xp++; 02052 } 02053 02055 template<class T> 02056 void exponentialMovingAverageUpdate(const TVec<T>& vec, const TVec<T>& x, T alpha) 02057 { 02058 int n=x.length(); 02059 if (vec.length()!=n) 02060 PLERROR("TVec::exponentialMovingAverageUpdate length_=%d and x has length_=%d", 02061 vec.length(),n); 02062 T* p=vec.data(); 02063 T* xp=x.data(); 02064 T one_minus_alpha = 1-alpha; 02065 for (int i=0;i<n;i++) 02066 p[i] = one_minus_alpha*p[i] + alpha*xp[i]; 02067 } 02068 02070 template<class T> 02071 void exponentialMovingVarianceUpdate(const TVec<T>& vec, const TVec<T>& x, const TVec<T>& mu, T alpha) 02072 { 02073 int n=x.length(); 02074 if (vec.length()!=n || vec.length()!=mu.length()) 02075 PLERROR("TVec::exponentialVarianceAverageUpdate length_=%d and" 02076 "x has length_=%d, mu has length() %d", 02077 vec.length(),n,mu.length()); 02078 T* p=vec.data(); 02079 T* xp=x.data(); 02080 T* mp=mu.data(); 02081 T one_minus_alpha = 1-alpha; 02082 for (int i=0;i<n;i++) 02083 { 02084 T dif = (xp[i]-mp[i]); 02085 p[i] = one_minus_alpha*p[i] + alpha*dif*dif; 02086 } 02087 } 02088 02089 02091 template<class T> 02092 void exponentialMovingSquareUpdate(const TVec<T>& vec, const TVec<T>& x, T alpha) 02093 { 02094 int n=x.length(); 02095 if (vec.length()!=n) 02096 PLERROR("TVec::exponentialMovingAverageUpdate length_=%d and x has length_=%d", 02097 vec.length(),n); 02098 T* p=vec.data(); 02099 T* xp=x.data(); 02100 T one_minus_alpha = 1-alpha; 02101 for (int i=0;i<n;i++) 02102 { 02103 T xpi = xp[i]; 02104 p[i] = one_minus_alpha*p[i] + alpha*xpi*xpi; 02105 } 02106 } 02107 02109 template<class T> 02110 void multiplyAcc(const TVec<T>& vec, const TVec<T>& x, const TVec<T>& y) 02111 { 02112 int n=x.length(); 02113 if (vec.length()!=n || y.length()!=n) 02114 PLERROR("TVec::multiplyAcc, this+=x*y: length_=%d, x.length_=%d, y.length_=%d", 02115 vec.length(),n,y.length()); 02116 T* p=vec.data(); 02117 T* xp=x.data(); 02118 T* yp=y.data(); 02119 for (int i=0;i<n;i++) 02120 p[i] += xp[i] * yp[i]; 02121 } 02122 02124 template<class T> 02125 void squareMultiplyAcc(const TVec<T>& vec, const TVec<T>& x, T scale) 02126 { 02127 int n=x.length(); 02128 if (vec.length()!=n) 02129 PLERROR("TVec::squareMultiplyAcc this has length_=%d and x has length_=%d", vec.length(),n); 02130 T* p=vec.data(); 02131 T* xp=x.data(); 02132 for (int i=0;i<n;i++) 02133 { 02134 T xpi = xp[i]; 02135 p[i] += scale * xpi * xpi; 02136 } 02137 } 02138 02140 template<class T> 02141 void squareAcc(const TVec<T>& vec, const TVec<T>& x) 02142 { 02143 int n=x.length(); 02144 if (vec.length()!=n) 02145 PLERROR("TVec::squareAcc this has length_=%d and x has length_=%d", vec.length(),n); 02146 T* p=vec.data(); 02147 T* xp=x.data(); 02148 for (int i=0;i<n;i++) 02149 { 02150 T xpi = xp[i]; 02151 p[i] += xpi * xpi; 02152 } 02153 } 02154 02156 template<class T> 02157 void squareSubtract(const TVec<T>& vec, const TVec<T>& x) 02158 { 02159 int n=x.length(); 02160 if (vec.length()!=n) 02161 PLERROR("TVec::squareDiff this has length_=%d and x has length_=%d", vec.length(),n); 02162 T* p=vec.data(); 02163 T* xp=x.data(); 02164 for (int i=0;i<n;i++) 02165 { 02166 T xpi = xp[i]; 02167 p[i] -= xpi * xpi; 02168 } 02169 } 02170 02172 template<class T> 02173 void diffSquareMultiplyAcc(const TVec<T>& vec, const TVec<T>& x, const TVec<T>& y, T scale) 02174 { 02175 int n=x.length(); 02176 if (vec.length()!=n || y.length()!=n) 02177 PLERROR("TVec::diffSquareMultiplyAcc this.length_=%d, x.length_=%d, y.length_=%d", 02178 vec.length(),n,y.length()); 02179 T* p=vec.data(); 02180 T* xp=x.data(); 02181 T* yp=y.data(); 02182 for (int i=0;i<n;i++) 02183 { 02184 T diff = xp[i]-yp[i]; 02185 p[i] += scale * diff * diff; 02186 } 02187 } 02188 02190 template<class T> 02191 void diffSquareMultiplyScaledAcc(const TVec<T>& vec, const TVec<T>& x, const TVec<T>& y, T fact1, T fact2) 02192 { 02193 int n=x.length(); 02194 if (vec.length()!=n || y.length()!=n) 02195 PLERROR("TVec::diffSquareMultiplyAcc this.length_=%d, x.length_=%d, y.length_=%d", 02196 vec.length(),n,y.length()); 02197 T* p=vec.data(); 02198 T* xp=x.data(); 02199 T* yp=y.data(); 02200 for (int i=0;i<n;i++) 02201 { 02202 T diff = xp[i]-yp[i]; 02203 p[i] = fact1 * p[i] + fact2 * diff * diff; 02204 } 02205 } 02206 02208 template <class T> 02209 void product(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v) 02210 { 02211 int l=m.width(); 02212 if (l!=v.length() || m.length()!=result.length()) 02213 PLERROR("product(TVec, TMat,TVec), incompatible arguments %dx%d times %d -> %d", 02214 m.length(),m.width(), v.length(),result.length()); 02215 T *rp = result.data(); 02216 T *vp = v.data(); 02217 for (int i=0;i<result.length();i++) 02218 { 02219 const T* mi = m[i]; 02220 T s = 0; 02221 for (int j=0;j<l;j++) 02222 s += mi[j] * vp[j]; 02223 rp[i] = s; 02224 } 02225 } 02226 02228 template <class T> 02229 void productAcc(const TVec<T>& vec, const TMat<T>& m, const TVec<T>& v) 02230 { 02231 int l = m.length(); 02232 int w = m.width(); 02233 #ifdef BOUNDCHECK 02234 if (w!=v.length() || l!=vec.length()) 02235 PLERROR("TVec::productAcc(TMat,TVec), incompatible arguments %dx%d times %d -> %d", 02236 m.length(),m.width(), v.length(),vec.length()); 02237 #endif 02238 T* rp = vec.data(); 02239 T* mp = m.data(); 02240 T* vdata = v.data(); 02241 int deltam = m.mod()-m.width(); 02242 for (int i=0;i<l;i++) 02243 { 02244 T *vp = vdata; 02245 T s = *rp; 02246 for (int j=0;j<w;j++) 02247 s += *mp++ * *vp++; 02248 *rp++ = s; 02249 mp += deltam; 02250 } 02251 } 02252 02256 template <class T> 02257 void transposeProduct(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v) 02258 { 02259 int l=m.length(); 02260 if (l!=v.length() || m.width()!=result.length()) 02261 PLERROR("TVec::transposeProduct(TMat,TVec), incompatible arguments %dx%d' times %d -> %d", 02262 m.length(),m.width(), v.length(),result.length()); 02263 T *rp = result.data(); 02264 T *vp = v.data(); 02265 result.clear(); 02266 for (int j=0;j<l;j++) 02267 { 02268 const T* mj = m[j]; 02269 T vj = vp[j]; 02270 for (int i=0;i<result.length();i++) 02271 rp[i] += mj[i] * vj; 02272 } 02273 } 02274 02276 template <class T> 02277 void transposeProductAcc(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v) 02278 { 02279 int l=m.length(); 02280 int w=m.width(); 02281 if (l!=v.length() || w!=result.length()) 02282 PLERROR("TVec::transposeProductAcc(TMat,TVec), incompatible arguments %dx%d' times %d -> %d", 02283 m.length(),m.width(), v.length(),result.length()); 02284 T* vecdata = result.data(); 02285 T* vp = v.data(); 02286 T* mp = m.data(); 02287 int deltam = m.mod()-m.width(); 02288 for (int j=0;j<l;j++) 02289 { 02290 T vj = *vp++; 02291 02292 /* 02293 T* rp = vecdata; 02294 for (int i=0;i<w;i++) 02295 *rp++ += vj * *mp++; 02296 mp += deltam; 02297 */ 02298 02299 if(vj!=0) 02300 { 02301 if(vj==1) 02302 { 02303 T* rp = vecdata; 02304 for (int i=0;i<w;i++) 02305 *rp++ += *mp++; 02306 mp += deltam; 02307 } 02308 else 02309 { 02310 T* rp = vecdata; 02311 for (int i=0;i<w;i++) 02312 *rp++ += vj * *mp++; 02313 mp += deltam; 02314 } 02315 } 02316 else mp += w + deltam; 02317 } 02318 } 02319 02321 template <class T> 02322 void transposeProductAcc(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v, T alpha) 02323 { 02324 int l=m.length(); 02325 int w=m.width(); 02326 if (l!=v.length() || w!=result.length()) 02327 PLERROR("TVec::transposeProductAcc(TMat,TVec), incompatible arguments %dx%d' times %d -> %d", 02328 m.length(),m.width(), v.length(),result.length()); 02329 T* vecdata = result.data(); 02330 T* vp = v.data(); 02331 T* mp = m.data(); 02332 int deltam = m.mod()-m.width(); 02333 for (int j=0;j<l;j++) 02334 { 02335 T vj = *vp++; 02336 02337 /* 02338 T* rp = vecdata; 02339 for (int i=0;i<w;i++) 02340 *rp++ += vj * *mp++; 02341 mp += deltam; 02342 */ 02343 02344 if(vj!=0) 02345 { 02346 if(vj==1) 02347 { 02348 T* rp = vecdata; 02349 for (int i=0;i<w;i++) 02350 *rp++ += alpha * *mp++; 02351 mp += deltam; 02352 } 02353 else 02354 { 02355 T* rp = vecdata; 02356 for (int i=0;i<w;i++) 02357 *rp++ += alpha * vj * *mp++; 02358 mp += deltam; 02359 } 02360 } 02361 else mp += w + deltam; 02362 } 02363 } 02364 02365 template <class T> 02366 void compressedTransposeProductAcc(const TVec<T>& result, const TMat<T>& m, char* comprbufvec) 02367 { 02368 cout<<"using kasjdlkadja"<<endl; 02369 union { double d; char c[8]; } uni; 02370 int l=m.length(),n, idx=0; 02371 unsigned char mode; 02372 cout<<"l="<<l<<endl; 02373 for(int i=0;i<l;i++) 02374 cout<<i<<":"<<char(comprbufvec[i])<<endl; 02375 while(l>0) 02376 { 02377 read_compr_mode_and_size_ptr(comprbufvec, mode, n); 02378 if(mode==0 || mode==1) 02379 { 02380 idx+=n; 02381 cout<<"0x"<<n<<" "; 02382 l-=n; 02383 if(mode==1) 02384 { 02385 --l; 02386 result+=m(idx++); // !!!!!! 02387 cout<<"1 "; 02388 } 02389 } 02390 else if(mode==2) 02391 { 02392 while(n--) 02393 { 02394 cout<<double(*comprbufvec)<<" "<<endl; 02395 result+= m(idx++) * double(*comprbufvec++); // !!!!!! 02396 02397 --l; 02398 } 02399 } 02400 else if(mode==3) 02401 { 02402 while(n--) 02403 { 02404 memcpy(uni.c,comprbufvec,sizeof(double)); 02405 cout<<double(uni.d)<<" "<<endl; 02406 comprbufvec+=8; 02407 result+= m(idx++) * uni.d; // !!!!!!! 02408 --l; 02409 } 02410 } 02411 else 02412 PLERROR("BUG IN binread_compressed: mode is only 2 bits, so how can it be other than 0,1,2,3 ?"); 02413 } 02414 02415 if(l!=0) 02416 PLERROR("In compressed_dot_product : l is not 0 at exit of function, wrong data?"); 02417 } 02418 02420 template<class T> 02421 void diagonalizedFactorsProduct(TMat<T>& result, const TMat<T>& U, const TVec<T> d, const TMat<T> V, bool accumulate=false) 02422 { 02423 #ifdef BOUNDCHECK 02424 if (result.length()!=U.length() || result.width()!=V.width() || U.width()!=d.length() || V.length()!=d.length()) 02425 PLERROR("diagonalizedFactorsProduct: incompatible arguments: (%dx%d)*(%d)*(%dx%d) --> (%dx%d)", 02426 U.length(),U.width(),d.length(),V.length(),V.width(),result.length(),result.width()); 02427 #endif 02428 int n1=U.length(); 02429 int n2=U.width(); 02430 int n3=V.width(); 02431 T *r_ij = result.data(); 02432 for (int i=0;i<n1;i++) 02433 { 02434 T *u_i = U[i]; 02435 for (int j=0;j<n3;j++,r_ij++) 02436 { 02437 T* d_k = d.data(); 02438 T res=0; 02439 for (int k=0;k<n2;k++,d_k++) 02440 res += *d_k * u_i[k] * V(k,j); 02441 if (accumulate) 02442 *r_ij += res; 02443 else 02444 *r_ij = res; 02445 } 02446 } 02447 } 02448 02454 template<class T> 02455 void diagonalizedFactorsProductBprop(const TMat<T>& dCdresult, const TMat<T>& U, const TVec<T> d, 02456 const TMat<T> V, TMat<T>& dCdU, TVec<T>& dCdd, TMat<T>& dCdV) 02457 { 02458 #ifdef BOUNDCHECK 02459 if (dCdU.length()!=U.length() || dCdU.width()!=U.width() || dCdd.length()!=d.length() 02460 || dCdV.length()!=V.length() || dCdV.width()!=V.width() || 02461 U.width()!=d.length() || V.length()!=d.length()) 02462 PLERROR("diagonalizedFactorsProductBprop: incompatible arguments"); 02463 #endif 02464 int n1=U.length(); 02465 int n2=U.width(); 02466 int n3=V.width(); 02467 T *dCdr_ij = dCdresult.data(); 02468 for (int i=0;i<n1;i++) 02469 { 02470 T *u_i = U[i]; 02471 T *dCdu_i = dCdU[i]; 02472 for (int j=0;j<n3;j++,dCdr_ij++) 02473 { 02474 T dcdr = *dCdr_ij; 02475 T* d_k = d.data(); 02476 T* dCdd_k = dCdd.data(); 02477 for (int k=0;k<n2;k++,d_k++,dCdd_k++) 02478 { 02479 T dk = *d_k; 02480 T u_ik = u_i[k]; 02481 T v_kj = V(k,j); 02482 dCdu_i[k] += dcdr * dk * v_kj; 02483 *dCdd_k += dcdr * u_ik * v_kj; 02484 dCdV(k,j) += dk * u_ik * dcdr; 02485 } 02486 } 02487 } 02488 } 02489 02491 template<class T> 02492 void diagonalizedFactorsProductTranspose(TMat<T>& result, const TMat<T>& U, const TVec<T> d, const TMat<T> V, bool accumulate=false) 02493 { 02494 #ifdef BOUNDCHECK 02495 if (result.length()!=U.length() || result.width()!=V.length() || U.width()!=d.length() || V.width()!=d.length()) 02496 PLERROR("diagonalizedFactorsProductTranspose: incompatible arguments: (%dx%d)*(%d)*(%dx%d)' --> (%dx%d)", 02497 U.length(),U.width(),d.length(),V.length(),V.width(),result.length(),result.width()); 02498 #endif 02499 int n1=U.length(); 02500 int n2=U.width(); 02501 int n3=V.length(); 02502 T *r_ij = result.data(); 02503 for (int i=0;i<n1;i++) 02504 { 02505 T *u_i = U[i]; 02506 for (int j=0;j<n3;j++,r_ij++) 02507 { 02508 T* d_k = d.data(); 02509 T* v_j = V[j]; 02510 T res=0; 02511 for (int k=0;k<n2;k++,d_k++) 02512 res += *d_k * u_i[k] * v_j[k]; 02513 if (accumulate) 02514 *r_ij += res; 02515 else 02516 *r_ij = res; 02517 } 02518 } 02519 } 02520 02521 // SINCE res[i,j] = sum_k U[i,k] d[k] V[j,k] ==> 02522 // gradients on dC/dU, dC/dd and dC/dV: 02523 // dC/dU[i,k] = sum_j dC/dres[i,j] d_k V[j,k] 02524 // dC/dd[k] = sum_{ij} dC/dres[i,j] U[i,k] V[j,k] 02525 // dC/dV[j,k] = sum_i dC/dres[i,j] d_k U[i,k] 02526 template<class T> 02527 void diagonalizedFactorsProductTransposeBprop(const TMat<T>& dCdresult, const TMat<T>& U, 02528 const TVec<T> d, const TMat<T> V, TMat<T>& dCdU, 02529 TVec<T>& dCdd, TMat<T>& dCdV) 02530 { 02531 #ifdef BOUNDCHECK 02532 if (dCdU.length()!=U.length() || dCdU.width()!=U.width() || dCdd.length()!=d.length() 02533 || dCdV.length()!=V.length() || dCdV.width()!=V.width() || 02534 U.width()!=d.length() || V.width()!=d.length()) 02535 PLERROR("diagonalizedFactorsProductTransposeBprop: incompatible arguments"); 02536 #endif 02537 int n1=U.length(); 02538 int n2=U.width(); 02539 int n3=V.length(); 02540 T *dCdr_ij = dCdresult.data(); 02541 for (int i=0;i<n1;i++) 02542 { 02543 T *u_i = U[i]; 02544 T *dCdu_i = dCdU[i]; 02545 for (int j=0;j<n3;j++,dCdr_ij++) 02546 { 02547 T* d_k = d.data(); 02548 T* dCdd_k = dCdd.data(); 02549 T* v_j = V[j]; 02550 T* dCdv_j = dCdV[j]; 02551 for (int k=0;k<n2;k++,d_k++,dCdd_k++) 02552 { 02553 T dcdr = *dCdr_ij; 02554 T dk = *d_k; 02555 T v_jk = v_j[k]; 02556 T u_ik = u_i[k]; 02557 dCdu_i[k] += dcdr * dk * v_jk; 02558 *dCdd_k += dcdr * u_ik * v_jk; 02559 dCdv_j[k] += dcdr * dk * u_ik; 02560 } 02561 } 02562 } 02563 } 02564 02566 template<class T> 02567 void diagonalizedFactorsTransposeProduct(TMat<T>& result, const TMat<T>& U, const TVec<T> d, const TMat<T> V, bool accumulate=false) 02568 { 02569 #ifdef BOUNDCHECK 02570 if (result.length()!=U.width() || result.width()!=V.width() || U.length()!=d.length() || V.length()!=d.length()) 02571 PLERROR("diagonalizedFactorsTransposeProduct: incompatible arguments: (%dx%d)'*(%d)*(%dx%d) --> (%dx%d)", 02572 U.length(),U.width(),d.length(),V.length(),V.width(),result.length(),result.width()); 02573 #endif 02574 int n1=U.width(); 02575 int n2=U.length(); 02576 int n3=V.width(); 02577 if (!accumulate) 02578 result.clear(); 02579 T* d_k = d.data(); 02580 for (int k=0;k<n2;k++,d_k++) 02581 { 02582 T *u_k = U[k]; 02583 T *v_k = V[k]; 02584 T *r_ij = result.data(); 02585 for (int i=0;i<n1;i++) 02586 { 02587 T u_ki = u_k[i]; 02588 for (int j=0;j<n3;j++,r_ij++) 02589 *r_ij += *d_k * u_ki * v_k[j]; 02590 } 02591 } 02592 } 02593 02594 // SINCE res[i,j] = sum_k U[k,i] d[k] V[k,j] ==> 02595 // gradients on dC/dU, dC/dd and dC/dV: 02596 // dC/dU[k,i] = d_k * sum_j dC/dres[i,j] V[k,j] 02597 // dC/dd[k] = sum_{ij} dC/dres[i,j] U[k,i] V[k,j] 02598 // dC/dV[k,j] = d_k sum_i dC/dres[i,j] U[k,i] 02599 template<class T> 02600 void diagonalizedFactorsTransposeProductBprop(const TMat<T>& dCdresult, const TMat<T>& U, 02601 const TVec<T> d, const TMat<T> V, TMat<T>& dCdU, 02602 TVec<T>& dCdd, TMat<T>& dCdV) 02603 { 02604 #ifdef BOUNDCHECK 02605 if (dCdU.length()!=U.length() || dCdU.width()!=U.width() || dCdd.length()!=d.length() 02606 || dCdV.length()!=V.length() || dCdV.width()!=V.width() || 02607 U.length()!=d.length() || V.length()!=d.length()) 02608 PLERROR("diagonalizedFactorsTransposeProductBprop: incompatible arguments"); 02609 #endif 02610 int n1=U.width(); 02611 int n2=U.length(); 02612 int n3=V.width(); 02613 T* d_k = d.data(); 02614 T* dCdd_k = dCdd.data(); 02615 for (int k=0;k<n2;k++,d_k++,dCdd_k++) 02616 { 02617 T dk = *d_k; 02618 T *u_k = U[k]; 02619 T *dCdu_k = dCdU[k]; 02620 T *v_k = V[k]; 02621 T *dCdv_k = dCdV[k]; 02622 T *dCdr_ij = dCdresult.data(); 02623 for (int i=0;i<n1;i++) 02624 { 02625 T u_ki = u_k[i]; 02626 T& dCdu_ki = dCdu_k[i]; 02627 for (int j=0;j<n3;j++,dCdr_ij++) 02628 { 02629 T dcdr = *dCdr_ij; 02630 T v_kj = v_k[j]; 02631 dCdu_ki += dcdr * dk * v_kj; 02632 *dCdd_k += dcdr * u_ki * v_kj; 02633 dCdv_k[j] += dcdr * dk * u_ki; 02634 } 02635 } 02636 } 02637 } 02638 02640 template<class T> 02641 void diagonalizedFactorsTransposeProductTranspose(TMat<T>& result, const TMat<T>& U, const TVec<T> d, const TMat<T> V, bool accumulate=false) 02642 { 02643 #ifdef BOUNDCHECK 02644 if (result.length()!=U.width() || result.width()!=V.length() || U.length()!=d.length() || V.width()!=d.length()) 02645 PLERROR("diagonalizedFactorsTransposeProductTranspose: incompatible arguments: (%dx%d)'*(%d)*(%dx%d)' --> (%dx%d)", 02646 U.length(),U.width(),d.length(),V.length(),V.width(),result.length(),result.width()); 02647 #endif 02648 int n1=U.width(); 02649 int n2=U.length(); 02650 int n3=V.length(); 02651 if (!accumulate) 02652 result.clear(); 02653 T* d_k = d.data(); 02654 for (int k=0;k<n2;k++,d_k++) 02655 { 02656 T *u_k = U[k]; 02657 T *r_ij = result.data(); 02658 for (int i=0;i<n1;i++) 02659 { 02660 T u_ki = u_k[i]; 02661 for (int j=0;j<n3;j++,r_ij++) 02662 *r_ij += *d_k * u_ki * V(j,k); 02663 } 02664 } 02665 } 02666 02667 // SINCE res[i,j] = sum_k U[k,i] d[k] V[j,k] ==> 02668 // gradients on dC/dU, dC/dd and dC/dV: 02669 // dC/dU[k,i] = d_k * sum_j dC/dres[i,j] V[j,k] 02670 // dC/dd[k] = sum_{ij} dC/dres[i,j] U[k,i] V[j,k] 02671 // dC/dV[j,k] = d_k * sum_i dC/dres[i,j] U[k,i] 02672 template<class T> 02673 void diagonalizedFactorsTransposeProductTransposeBprop(const TMat<T>& dCdresult, const TMat<T>& U, 02674 const TVec<T> d, const TMat<T> V, TMat<T>& dCdU, 02675 TVec<T>& dCdd, TMat<T>& dCdV) 02676 { 02677 #ifdef BOUNDCHECK 02678 if (dCdU.length()!=U.length() || dCdU.width()!=U.width() || dCdd.length()!=d.length() 02679 || dCdV.length()!=V.length() || dCdV.width()!=V.width() || 02680 U.length()!=d.length() || V.width()!=d.length()) 02681 PLERROR("diagonalizedFactorsTransposeProductTransposeBprop: incompatible arguments"); 02682 #endif 02683 int n1=U.width(); 02684 int n2=U.length(); 02685 int n3=V.length(); 02686 T* d_k = d.data(); 02687 T* dCdd_k = dCdd.data(); 02688 for (int k=0;k<n2;k++,d_k++,dCdd_k++) 02689 { 02690 T dk = *d_k; 02691 T *u_k = U[k]; 02692 T *dCdu_k = dCdU[k]; 02693 T *dCdr_ij = dCdresult.data(); 02694 for (int i=0;i<n1;i++) 02695 { 02696 T u_ki = u_k[i]; 02697 T& dCdu_ki = dCdu_k[i]; 02698 for (int j=0;j<n3;j++,dCdr_ij++) 02699 { 02700 T dcdr = *dCdr_ij; 02701 T v_jk = V(j,k); 02702 dCdu_ki += dcdr * dk * v_jk; 02703 *dCdd_k += dcdr * u_ki * v_jk; 02704 dCdV(j,k) += dcdr * dk * u_ki; 02705 } 02706 } 02707 } 02708 } 02709 02710 // ---------- these were previously methods of TMat --------------- 02711 02713 template<class T> 02714 T matRowDotVec(const TMat<T>& mat, int i, const TVec<T> v) 02715 { 02716 #ifdef BOUNDCHECK 02717 if (v.length()!=mat.width()) 02718 PLERROR("dotRow(%d,v), v.length_=%d != matrix width_=%d", 02719 i,v.length(),mat.width()); 02720 #endif 02721 T s = 0; 02722 T* rowi = mat.rowdata(i); 02723 T* v_=v.data(); 02724 for (int j=0;j<mat.width();j++) 02725 s += rowi[j] * v_[j]; 02726 return s; 02727 } 02728 02730 template<class T> 02731 T matColumnDotVec(const TMat<T>& mat, int j, const TVec<T> v) 02732 { 02733 #ifdef BOUNDCHECK 02734 if (v.length()!=mat.length()) 02735 PLERROR("dotColumn(%d,v), v.length_=%d != matrix length_=%d", 02736 j,v.length(),mat.length()); 02737 #endif 02738 T s = 0; 02739 T* colj = mat.data()+j; 02740 T* v_=v.data(); 02741 for (int i=0;i<mat.length();i++, colj+=mat.mod()) 02742 s += *colj * v_[i]; 02743 return s; 02744 } 02745 02747 template<class T> 02748 void matRowsDots(TVec<T> v, const TMat<T>& A, const TMat<T>& B) 02749 { 02750 #ifdef BOUNDCHECK 02751 if (A.length()!=v.length()) 02752 PLERROR("matRowDotsVec(v,A,B): v.length_=%d != A.length_=%d", 02753 v.length(),A.length()); 02754 if (A.length()!=B.length()) 02755 PLERROR("matRowDotsVec(v,A,B): A.length_=%d != B.length_=%d", 02756 A.length(),B.length()); 02757 if (A.width()!=B.width()) 02758 PLERROR("matRowDotsVec(v,A,B): A.width_=%d != B.width_=%d", 02759 A.width(),B.width()); 02760 #endif 02761 int l=A.length(), w=A.width(); 02762 T* vi = v.data(); 02763 for (int i=0;i<l;i++) 02764 { 02765 T s = 0; 02766 T* Aij = A[i]; 02767 T* Bij = B[i]; 02768 for (int j=0;j<w;j++) 02769 s += *Aij++ * *Bij++; 02770 *vi++ = s; 02771 } 02772 } 02773 02775 template<class T> 02776 void matRowsDotsAcc(TVec<T> v, const TMat<T>& A, const TMat<T>& B) 02777 { 02778 #ifdef BOUNDCHECK 02779 if (A.length()!=v.length()) 02780 PLERROR("matRowDotsVec(v,A,B): v.length_=%d != A.length_=%d", 02781 v.length(),A.length()); 02782 if (A.length()!=B.length()) 02783 PLERROR("matRowDotsVec(v,A,B): A.length_=%d != B.length_=%d", 02784 A.length(),B.length()); 02785 if (A.width()!=B.width()) 02786 PLERROR("matRowDotsVec(v,A,B): A.width_=%d != B.width_=%d", 02787 A.width(),B.width()); 02788 #endif 02789 int l=A.length(), w=A.width(); 02790 T* vi = v.data(); 02791 for (int i=0;i<l;i++) 02792 { 02793 T s = 0; 02794 T* Aij = A[i]; 02795 T* Bij = B[i]; 02796 for (int j=0;j<w;j++) 02797 s += *Aij++ * *Bij++; 02798 *vi++ += s; 02799 } 02800 } 02801 02804 template<class T> 02805 void fillItSymmetric(const TMat<T>& mat) { 02806 int m = mat.mod(); 02807 T* mat_data_to_fill; 02808 T* mat_data_to_copy; 02809 for (int i = 0; i < mat.length(); i++) { 02810 mat_data_to_fill = mat[i]; 02811 mat_data_to_copy = &mat[0][i]; 02812 for (int j = 0; j < i; j++) { 02813 *(mat_data_to_fill++) = *mat_data_to_copy; 02814 mat_data_to_copy += m; 02815 } 02816 } 02817 } 02818 02819 template<class T> 02820 void makeItSymmetric(const TMat<T>& mat, T max_dif) 02821 { 02822 if (!mat.isSquare()) 02823 PLERROR("at void makeItSymmetric, the matrix is not even square\n"); 02824 T dif; 02825 T value; 02826 bool warning_flag = false; 02827 for (int i=0; i<mat.length()-1 ; i++) 02828 for (int j=i+1; j<mat.width(); j++) 02829 { 02830 dif = std::abs(mat[i][j] - mat[j][i]); 02831 if (dif > max_dif) 02832 { 02833 max_dif = dif; 02834 warning_flag = true; 02835 } 02836 value = (mat[i][j] + mat[j][i])/2; 02837 mat[i][j] = value; mat[j][i] = value; 02838 } 02839 if (warning_flag) 02840 PLWARNING("At void makeItSymmetric, the maximum difference %f is not affordable\n", max_dif); 02841 } 02842 02843 // y[i] = sum_j A[i,j] x[j] 02844 02845 template<class T> 02846 void product(const TMat<T>& mat, const TVec<T>& x, TVec<T>& y) 02847 { 02848 #ifdef BOUNDCHECK 02849 if (mat.length()!=y.length() || mat.width()!=x.length()) 02850 PLERROR("TMat(%d,%d)::product(TVec& x(%d),TVec& y(%d)), incompatible arguments", 02851 mat.length(),mat.width(),x.length(),y.length()); 02852 #endif 02853 T* x_=x.data(); 02854 T* y_=y.data(); 02855 for (int i=0;i<mat.length();i++) 02856 { 02857 T* Ai = mat[i]; 02858 T yi = 0; 02859 for (int j=0;j<mat.width();j++) 02860 yi += Ai[j] * x_[j]; 02861 y_[i]=yi; 02862 } 02863 } 02864 02865 // result[i,j] = sum_k m1[i,k] * m2[k,j] 02866 template<class T> 02867 void product(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 02868 { 02869 #ifdef BOUNDCHECK 02870 if (m1.width()!=m2.length() || mat.length()!=m1.length() || mat.width()!=m2.width()) 02871 PLERROR("product(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d", 02872 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 02873 #endif 02874 int n=m1.length(); 02875 int m=m1.width(); 02876 int l=m2.width(); 02877 for (int i=0;i<n;i++) 02878 { 02879 const T* m1i = m1[i]; 02880 T* mi = mat[i]; 02881 for (int j=0;j<l;j++) 02882 { 02883 T s=0; 02884 const T* m2kj = m2.data()+j; 02885 for (int k=0;k<m;k++,m2kj+=m2.mod()) 02886 s += m1i[k] * (*m2kj); 02887 mi[j] = s; 02888 } 02889 } 02890 } 02891 02892 // result[i,j] += sum_k m1[i,k] * m2[k,j] 02893 template<class T> 02894 void productAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 02895 { 02896 #ifdef BOUNDCHECK 02897 if (m1.width()!=m2.length() || mat.length()!=m1.length() || mat.width()!=m2.width()) 02898 PLERROR("productAcc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d", 02899 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 02900 #endif 02901 int n=m1.length(); 02902 int m=m1.width(); 02903 int l=m2.width(); 02904 for (int i=0;i<n;i++) 02905 { 02906 const T* m1i = m1[i]; 02907 T* mi = mat[i]; 02908 for (int j=0;j<l;j++) 02909 { 02910 T s=0; 02911 T* m2kj = m2.data()+j; 02912 for (int k=0;k<m;k++,m2kj+=m2.mod()) 02913 s += m1i[k] * (*m2kj); 02914 mi[j] += s; 02915 } 02916 } 02917 } 02918 02919 // result[i,j] += sum_k m1[i,k] * m2[k,j]^2 02920 template<class T> 02921 void product2Acc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 02922 { 02923 #ifdef BOUNDCHECK 02924 if (m1.width()!=m2.length() || mat.length()!=m1.length() || mat.width()!=m2.width()) 02925 PLERROR("product2Acc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d", 02926 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 02927 #endif 02928 int n=m1.length(); 02929 int m=m1.width(); 02930 int l=m2.width(); 02931 for (int i=0;i<n;i++) 02932 { 02933 const T* m1i = m1[i]; 02934 T* mi = mat[i]; 02935 for (int j=0;j<l;j++) 02936 { 02937 T s=0; 02938 T* m2kj = m2.data()+j; 02939 for (int k=0;k<m;k++,m2kj+=m2.mod()) 02940 s += m1i[k] * (*m2kj) * (*m2kj); 02941 mi[j] += s; 02942 } 02943 } 02944 } 02945 02946 // result[i,j] += sum_k m1[i,k]^2 * m2[k,j] 02947 template<class T> 02948 void squareProductAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 02949 { 02950 #ifdef BOUNDCHECK 02951 if (m1.width()!=m2.length() || mat.length()!=m1.length() || mat.width()!=m2.width()) 02952 PLERROR("squareProductAcc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d", 02953 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 02954 #endif 02955 int n=m1.length(); 02956 int m=m1.width(); 02957 int l=m2.width(); 02958 for (int i=0;i<n;i++) 02959 { 02960 const T* m1i = m1[i]; 02961 T* mi = mat[i]; 02962 for (int j=0;j<l;j++) 02963 { 02964 T s=0; 02965 T* m2kj = m2.data()+j; 02966 for (int k=0;k<m;k++,m2kj+=m2.mod()) 02967 { 02968 T m1ik=m1i[k]; 02969 s += m1ik*m1ik * (*m2kj); 02970 } 02971 mi[j] += s; 02972 } 02973 } 02974 } 02975 02976 // result[i][j] = v1[i] * v2[j] 02977 02978 template<class T> 02979 void externalProduct(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2) 02980 { 02981 #ifdef BOUNDCHECK 02982 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 02983 PLERROR("externalProduct(Vec,Vec), incompatible arguments %dx%d= %d times %d", 02984 mat.length(),mat.width(),v1.length(), v2.length()); 02985 #endif 02986 const T* v_1=v1.data(); 02987 const T* v_2=v2.data(); 02988 for (int i=0;i<mat.length();i++) 02989 { 02990 T* mi = mat[i]; 02991 T v1i = v_1[i]; 02992 for (int j=0;j<mat.width();j++) 02993 mi[j] = v1i * v_2[j]; 02994 } 02995 } 02996 02997 // mat[i][j] += v1[i] * v2[j] 02998 template<class T> 02999 void externalProductAcc(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2) 03000 { 03001 #ifdef BOUNDCHECK 03002 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03003 PLERROR("externalProductAcc(Vec,Vec), incompatible arguments %dx%d= %d times %d", 03004 mat.length(),mat.width(),v1.length(), v2.length()); 03005 #endif 03006 03007 T* v_1=v1.data(); 03008 T* v_2=v2.data(); 03009 T* mp = mat.data(); 03010 int l = mat.length(); 03011 int w = mat.width(); 03012 03013 if(mat.isCompact()) 03014 { 03015 T* pv1 = v_1; 03016 for(int i=0; i<l; i++) 03017 { 03018 T* pv2 = v_2; 03019 T val = *pv1++; 03020 for(int j=0; j<w; j++) 03021 *mp++ += val * *pv2++; 03022 } 03023 } 03024 else 03025 { 03026 cerr << "!"; 03027 for (int i=0;i<l;i++) 03028 { 03029 T* mi = mat[i]; 03030 T v1i = v_1[i]; 03031 for (int j=0;j<w;j++) 03032 mi[j] += v1i * v_2[j]; 03033 } 03034 } 03035 } 03036 03037 // mat[i][j] += gamma * v1[i] * v2[j] 03038 template<class T> 03039 void externalProductScaleAcc(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2, T gamma) 03040 { 03041 #ifdef BOUNDCHECK 03042 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03043 PLERROR("externalProductScaleAcc(Vec,Vec), incompatible arguments %dx%d= %d times %d", 03044 mat.length(),mat.width(),v1.length(), v2.length()); 03045 #endif 03046 const T* v_1=v1.data(); 03047 const T* v_2=v2.data(); 03048 for (int i=0;i<mat.length();i++) 03049 { 03050 T* mi = mat[i]; 03051 T v1i = v_1[i]; 03052 for (int j=0;j<mat.width();j++) 03053 mi[j] += gamma * v1i * v_2[j]; 03054 } 03055 } 03056 03057 // mat[i][j] = alpha * mat[i][j] + gamma * v1[i] * v2[j] 03058 template<class T> 03059 void externalProductScaleAcc(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2, T gamma, T alpha) 03060 { 03061 #ifdef BOUNDCHECK 03062 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03063 PLERROR("externalProductScaleAcc(Vec,Vec), incompatible arguments %dx%d= %d times %d", 03064 mat.length(),mat.width(),v1.length(), v2.length()); 03065 #endif 03066 const T* v_1=v1.data(); 03067 const T* v_2=v2.data(); 03068 for (int i=0;i<mat.length();i++) 03069 { 03070 T* mi = mat[i]; 03071 T v1i = v_1[i]; 03072 for (int j=0;j<mat.width();j++) 03073 mi[j] = alpha*mi[j] + gamma * v1i * v_2[j]; 03074 } 03075 } 03076 03077 // result[i,j] = sum_k m1[i,k] * m2[j,k] 03078 template<class T> 03079 void productTranspose(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03080 { 03081 #ifdef BOUNDCHECK 03082 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 03083 PLERROR("productTranspose(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 03084 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03085 #endif 03086 int n=m1.length(); 03087 int m=m1.width(); 03088 int l=m2.length(); 03089 for (int i=0;i<n;i++) 03090 { 03091 const T* m1i = m1[i]; 03092 T* mi = mat[i]; 03093 for (int j=0;j<l;j++) 03094 { 03095 T s=0; 03096 const T* m2j = m2[j]; 03097 for (int k=0;k<m;k++) 03098 s += m1i[k] * m2j[k]; 03099 mi[j] = s; 03100 } 03101 } 03102 } 03103 03104 // result[i,j] = sum_k m1[i,k]^2 * m2[j,k] 03105 template<class T> 03106 void squareProductTranspose(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03107 { 03108 #ifdef BOUNDCHECK 03109 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 03110 PLERROR("squareProductTranspose(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 03111 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03112 #endif 03113 int n=m1.length(); 03114 int m=m1.width(); 03115 int l=m2.length(); 03116 for (int i=0;i<n;i++) 03117 { 03118 const T* m1i = m1[i]; 03119 T* mi = mat[i]; 03120 for (int j=0;j<l;j++) 03121 { 03122 T s=0; 03123 const T* m2j = m2[j]; 03124 for (int k=0;k<m;k++) 03125 { 03126 T m1ik=m1i[k]; 03127 s += m1ik*m1ik * m2j[k]; 03128 } 03129 mi[j] = s; 03130 } 03131 } 03132 } 03133 03134 // result[i,j] = sum_k m1[i,k] * m2[j,k]^2 03135 template<class T> 03136 void product2Transpose(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03137 { 03138 #ifdef BOUNDCHECK 03139 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 03140 PLERROR("product2Transpose(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 03141 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03142 #endif 03143 int n=m1.length(); 03144 int m=m1.width(); 03145 int l=m2.length(); 03146 for (int i=0;i<n;i++) 03147 { 03148 const T* m1i = m1[i]; 03149 T* mi = mat[i]; 03150 for (int j=0;j<l;j++) 03151 { 03152 T s=0; 03153 const T* m2j = m2[j]; 03154 for (int k=0;k<m;k++) 03155 { 03156 T m2jk=m2j[k]; 03157 s += m1i[k] * m2jk*m2jk; 03158 } 03159 mi[j] = s; 03160 } 03161 } 03162 } 03163 03164 // result[i,j] += sum_k m1[i,k] * m2[j,k] 03165 template<class T> 03166 void productTransposeAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03167 { 03168 #ifdef BOUNDCHECK 03169 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 03170 PLERROR("productTransposeAcc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 03171 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03172 #endif 03173 int n=m1.length(); 03174 int m=m1.width(); 03175 int l=m2.length(); 03176 for (int i=0;i<n;i++) 03177 { 03178 const T* m1i = m1[i]; 03179 T* mi = mat[i]; 03180 for (int j=0;j<l;j++) 03181 { 03182 T s=0; 03183 const T* m2j = m2[j]; 03184 for (int k=0;k<m;k++) 03185 s += m1i[k] * m2j[k]; 03186 mi[j] += s; 03187 } 03188 } 03189 } 03190 03191 // result[i,j] += sum_k m1[i,k] * m2[j,k]^2 03192 template<class T> 03193 void product2TransposeAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03194 { 03195 #ifdef BOUNDCHECK 03196 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 03197 PLERROR("product2TransposeAcc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 03198 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03199 #endif 03200 int n=m1.length(); 03201 int m=m1.width(); 03202 int l=m2.length(); 03203 for (int i=0;i<n;i++) 03204 { 03205 const T* m1i = m1[i]; 03206 T* mi = mat[i]; 03207 for (int j=0;j<l;j++) 03208 { 03209 T s=0; 03210 const T* m2j = m2[j]; 03211 for (int k=0;k<m;k++) 03212 { 03213 T m2jk=m2j[k]; 03214 s += m1i[k] * m2jk*m2jk; 03215 } 03216 mi[j] += s; 03217 } 03218 } 03219 } 03220 03221 // result[i,j] += sum_k m1[i,k]^2 * m2[j,k] 03222 template<class T> 03223 void squareProductTransposeAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03224 { 03225 #ifdef BOUNDCHECK 03226 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 03227 PLERROR("squareProductTransposeAcc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 03228 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03229 #endif 03230 int n=m1.length(); 03231 int m=m1.width(); 03232 int l=m2.length(); 03233 for (int i=0;i<n;i++) 03234 { 03235 const T* m1i = m1[i]; 03236 T* mi = mat[i]; 03237 for (int j=0;j<l;j++) 03238 { 03239 T s=0; 03240 const T* m2j = m2[j]; 03241 for (int k=0;k<m;k++) 03242 { 03243 T m1ik=m1i[k]; 03244 s += m1ik*m1ik * m2j[k]; 03245 } 03246 mi[j] += s; 03247 } 03248 } 03249 } 03250 03251 // result[i,j] = sum_k m1[k,i] * m2[k,j] 03252 template<class T> 03253 void transposeProduct(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03254 { 03255 int n=m1.width(); 03256 int m=m1.length(); 03257 int l=m2.width(); 03258 #ifdef BOUNDCHECK 03259 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 03260 PLERROR("transposeProduct(Mat,Mat), incompatible arguments " 03261 "%dx%d' times %dx%d into %dx%d", 03262 m1.length(),m1.width(), m2.length(),m2.width(), mat.length(), mat.width()); 03263 #endif 03264 mat.clear(); 03265 for (int i=0;i<n;i++) 03266 { 03267 T* m1ki = m1.data()+i; 03268 T* mi = mat[i]; 03269 for (int k=0;k<m;k++,m1ki+=m1.mod()) 03270 { 03271 const T* m2k = m2[k]; 03272 T m1_ki = *m1ki; 03273 for (int j=0;j<l;j++) 03274 mi[j] += m1_ki * m2k[j]; 03275 } 03276 } 03277 } 03278 03279 // result[i,j] = sum_k m1[k,i] * m2[k,j]^2 03280 template<class T> 03281 void transposeProduct2(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03282 { 03283 int n=m1.width(); 03284 int m=m1.length(); 03285 int l=m2.width(); 03286 #ifdef BOUNDCHECK 03287 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 03288 PLERROR("transposeProduct2(Mat,Mat), incompatible arguments " 03289 "%dx%d' times %dx%d into %dx%d", 03290 m1.length(),m1.width(), m2.length(),m2.width(), mat.length(), mat.width()); 03291 #endif 03292 mat.clear(); 03293 for (int i=0;i<n;i++) 03294 { 03295 T* m1ki = m1.data()+i; 03296 T* mi = mat[i]; 03297 for (int k=0;k<m;k++,m1ki+=m1.mod()) 03298 { 03299 const T* m2k = m2[k]; 03300 T m1_ki = *m1ki; 03301 for (int j=0;j<l;j++) 03302 { 03303 T m2kj=m2k[j]; 03304 mi[j] += m1_ki * m2kj*m2kj; 03305 } 03306 } 03307 } 03308 } 03309 03310 // result[i,j] += sum_k m1[k,i] * m2[k,j] 03311 template<class T> 03312 void transposeProductAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03313 { 03314 int n=m1.width(); 03315 int m=m1.length(); 03316 int l=m2.width(); 03317 #ifdef BOUNDCHECK 03318 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 03319 PLERROR("transposeProductAcc(Mat,Mat), incompatible arguments " 03320 "%dx%d' times %dx%d into %dx%d", 03321 m1.length(),m1.width(), m2.length(),m2.width(), mat.length(), mat.width()); 03322 #endif 03323 for (int i=0;i<n;i++) 03324 { 03325 T* m1ki = m1.data()+i; 03326 T* mi = mat[i]; 03327 for (int k=0;k<m;k++,m1ki+=m1.mod()) 03328 { 03329 const T* m2k = m2[k]; 03330 T m1_ki = *m1ki; 03331 for (int j=0;j<l;j++) 03332 mi[j] += m1_ki * m2k[j]; 03333 } 03334 } 03335 } 03336 03337 // result[i,j] += sum_k m1[k,i] * m2[k,j]^2 03338 template<class T> 03339 void transposeProduct2Acc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03340 { 03341 int n=m1.width(); 03342 int m=m1.length(); 03343 int l=m2.width(); 03344 #ifdef BOUNDCHECK 03345 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 03346 PLERROR("transposeProduct2Acc(Mat,Mat), incompatible arguments " 03347 "%dx%d' times %dx%d into %dx%d", 03348 m1.length(),m1.width(), m2.length(),m2.width(), mat.length(), mat.width()); 03349 #endif 03350 for (int i=0;i<n;i++) 03351 { 03352 T* m1ki = m1.data()+i; 03353 T* mi = mat[i]; 03354 for (int k=0;k<m;k++,m1ki+=m1.mod()) 03355 { 03356 const T* m2k = m2[k]; 03357 T m1_ki = *m1ki; 03358 for (int j=0;j<l;j++) 03359 { 03360 T m2kj = m2k[j]; 03361 mi[j] += m1_ki * m2kj * m2kj; 03362 } 03363 } 03364 } 03365 } 03366 03367 // result[i,j] = sum_k m1[k,i] * m2[j,k] 03368 template<class T> 03369 void transposeTransposeProduct(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03370 { 03371 #ifdef BOUNDCHECK 03372 if (m1.length()!=m2.width()) 03373 PLERROR("transposeTransposeProduct(Mat,Mat), incompatible arguments %dx%d' times %dx%d'", 03374 m1.length(),m1.width(), m2.length(),m2.width()); 03375 #endif 03376 int n=m1.width(); 03377 int m=m1.length(); 03378 int l=m2.length(); 03379 for (int i=0;i<n;i++) 03380 { 03381 T* m1ki0 = m1.data()+i; 03382 T* mi = mat[i]; 03383 for (int j=0;j<l;j++) 03384 { 03385 T s=0; 03386 const T* m2j = m2[j]; 03387 T* m1ki = m1ki0; 03388 for (int k=0;k<m;k++,m1ki+=m1.mod()) 03389 s += (*m1ki) * m2j[k]; 03390 mi[j] = s; 03391 } 03392 } 03393 } 03394 03395 // result[i,j] += sum_k m1[k,i] * m2[j,k] 03396 template<class T> 03397 void transposeTransposeProductAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03398 { 03399 if (m1.length()!=m2.width()) 03400 PLERROR("transposeTransposeProductAcc(Mat,Mat), incompatible arguments %dx%d' times %dx%d", 03401 m1.length(),m1.width(), m2.length(),m2.width()); 03402 int n=m1.width(); 03403 int m=m1.length(); 03404 int l=m2.width(); 03405 for (int i=0;i<n;i++) 03406 { 03407 T* m1ki0 = m1.data()+i; 03408 T* mi = mat[i]; 03409 for (int j=0;j<l;j++) 03410 { 03411 T s=0; 03412 const T* m2j = m2[j]; 03413 T* m1ki = m1ki0; 03414 for (int k=0;k<m;k++,m1ki+=m1.mod()) 03415 s += (*m1ki) * m2j[k]; 03416 mi[j] += s; 03417 } 03418 } 03419 } 03420 03421 template<class T> 03422 T trace(const TMat<T>& mat) 03423 { 03424 if (!mat.isSquare()) 03425 PLERROR( "In trace()\nThe matrix must be square." ); 03426 T tr = mat.firstElement(); 03427 for ( int i = 1; i < mat.length(); i++ ) 03428 tr += mat(i,i); 03429 return tr; 03430 } 03431 03433 template<class T> 03434 void regularizeMatrix(const TMat<T>& mat, T tolerance) 03435 { 03436 static T reg; 03437 static T* k; 03438 static int shift; 03439 reg = tolerance * trace(mat); 03440 k = mat.data(); 03441 shift = mat.mod() + 1; 03442 for (int i = 0; i < mat.length(); i++) { 03443 *k += reg; 03444 k += shift; 03445 } 03446 } 03447 03448 03449 template<class T> 03450 void makeRowsSumTo1(const TMat<T>& mat) 03451 { 03452 for (int i = 0; i < mat.length(); i++) 03453 { 03454 TVec<T> row_i = mat(i); 03455 divide(row_i, sum(row_i), row_i); 03456 } 03457 } 03458 03459 03460 // result[i,j] = x[i,j]*scale; 03461 template<class T> 03462 void multiply(const TMat<T>& result, const TMat<T>& x, T scale) 03463 { 03464 #ifdef BOUNDCHECK 03465 if (result.length()!=x.length() || result.width()!=x.width()) 03466 PLERROR("multiply incompatible dimensions: %dx%d <- %dx%d", 03467 result.length(),result.width(),x.length(),x.width()); 03468 #endif 03469 if(result.isCompact() && x.isCompact()) 03470 { 03471 typename TMat<T>::compact_iterator itm = result.compact_begin(); 03472 typename TMat<T>::compact_iterator itmend = result.compact_end(); 03473 typename TMat<T>::compact_iterator itx = x.compact_begin(); 03474 for(; itm!=itmend; ++itm, ++itx) 03475 *itm = *itx * scale; 03476 } 03477 else // use non-compact iterators 03478 { 03479 typename TMat<T>::iterator itm = result.begin(); 03480 typename TMat<T>::iterator itmend = result.end(); 03481 typename TMat<T>::iterator itx = x.begin(); 03482 for(; itm!=itmend; ++itm, ++itx) 03483 *itm = *itx * scale; 03484 } 03485 } 03486 03487 template<class T> 03488 inline TMat<T> operator*(const TMat<T>& m, const T& scalar) 03489 { 03490 TMat<T> result(m.length(),m.width()); 03491 multiply(result, m, scalar); 03492 return result; 03493 } 03494 03495 template<class T> 03496 inline TMat<T> operator*(const T& scalar, const TMat<T>& m) 03497 { return m * scalar;} 03498 03499 // Will not work properly for integers... 03500 template<class T> 03501 inline TMat<T> operator/(const TMat<T>& m, const T& scalar) 03502 { return m * (T(1)/scalar); } 03503 03504 // result[i,j] += x[i,j]*scale; 03505 template<class T> 03506 void multiplyAcc(const TMat<T>& mat, const TMat<T>& x, T scale) 03507 { 03508 #ifdef BOUNDCHECK 03509 if (mat.length()!=x.length() || mat.width()!=x.width()) 03510 PLERROR("multiplyAcc incompatible dimensions: %dx%d <- %dx%d", 03511 mat.length(),mat.width(),x.length(),x.width()); 03512 #endif 03513 if(mat.isCompact() && x.isCompact()) 03514 { 03515 typename TMat<T>::compact_iterator itm = mat.compact_begin(); 03516 typename TMat<T>::compact_iterator itmend = mat.compact_end(); 03517 typename TMat<T>::compact_iterator itx = x.compact_begin(); 03518 for(; itm!=itmend; ++itm, ++itx) 03519 *itm += *itx * scale; 03520 } 03521 else // use non-compact iterators 03522 { 03523 typename TMat<T>::iterator itm = mat.begin(); 03524 typename TMat<T>::iterator itmend = mat.end(); 03525 typename TMat<T>::iterator itx = x.begin(); 03526 for(; itm!=itmend; ++itm, ++itx) 03527 *itm += *itx * scale; 03528 } 03529 } 03530 03531 // result[i,j] += x[i,j]*y[i,j]; 03532 template<class T> 03533 void multiplyAcc(const TMat<T>& mat, const TMat<T>& x, const TMat<T>& y) 03534 { 03535 int n=mat.length()*mat.width(); 03536 if (mat.length()!=x.length() || mat.width()!=x.width() || y.length()!=mat.length() || y.width()!=mat.width()) 03537 PLERROR("multiplyAcc this has size=%dx%d, x is %dx%d, y is %dx%d", 03538 mat.length(),mat.width(),x.length(),x.width(),y.length(),y.width()); 03539 T* p=mat.data(); 03540 T* xp=x.data(); 03541 T* yp=y.data(); 03542 for (int i=0;i<n;i++) 03543 p[i] += xp[i] * yp[i]; 03544 } 03545 03546 // result[i,j] += x[i,j]*x[i,j]*scale; 03547 template<class T> 03548 void squareMultiplyAcc(const TMat<T>& mat, const TMat<T>& x, T scale) 03549 { 03550 int n=x.length()*x.width(); 03551 if (mat.length()*mat.width()!=n) 03552 PLERROR("squareMultiplyAcc this has size=%d and x has size=%d", 03553 mat.width()*mat.length(),n); 03554 T* p=mat.data(); 03555 T* xp=x.data(); 03556 for (int i=0;i<n;i++) 03557 { 03558 T xpi = xp[i]; 03559 p[i] += scale * xpi * xpi; 03560 } 03561 } 03562 03563 // result[i,j] += (x[i,j]-y[i,j])^2*scale; 03564 template<class T> 03565 void diffSquareMultiplyAcc(const TMat<T>& mat, const TMat<T>& x, const TMat<T>& y, T scale) 03566 { 03567 int n=x.length()*x.width(); 03568 if (mat.length()*mat.width()!=n) 03569 PLERROR("diffSquareMultiplyAcc this has size=%d and x has size=%d", 03570 mat.width()*mat.length(),n); 03571 T* p=mat.data(); 03572 T* xp=x.data(); 03573 T* yp=y.data(); 03574 for (int i=0;i<n;i++) 03575 { 03576 T diff = (xp[i]-yp[i]); 03577 p[i] += scale * diff * diff; 03578 } 03579 } 03580 03581 03582 03583 template<class T> 03584 TVec<T> selectAndOrder(const TMat<T>& mat, int pos, int colnum=0) 03585 { 03586 #ifdef BOUNDCHECK 03587 if (colnum<0 || colnum>=mat.width()) PLERROR("Bad column number (%d)", colnum); 03588 if (pos<0 || pos>=mat.length()) PLERROR("Bad position (%d)", pos); 03589 #endif 03590 03591 int l=0; 03592 int h=mat.length()-1; 03593 TMat<T> v = mat.column(colnum); 03594 03595 while (l<h) 03596 { 03597 T p = v((l+h)/2,0); 03598 int x = l; 03599 int y = h; 03600 03601 do 03602 { 03603 while (v(x,0)<p) x++; 03604 while (p<v(y,0)) y--; 03605 if (x<=y) 03606 { 03607 mat.swapRows(x,y); 03608 x++; 03609 y--; 03610 } 03611 } while (x<=y); 03612 03613 if (pos>=x) l=x; 03614 else h=x-1; 03615 } 03616 03617 return mat(l); 03618 } 03619 03620 03621 // result[i,i] += lambda 03622 template<class T> 03623 void addToDiagonal(const TMat<T>& mat, T lambda) 03624 { 03625 T *d = mat.data(); 03626 for (int i=0;i<mat.length();i++,d+=mat.mod()+1) *d+=lambda; 03627 } 03628 03629 03630 03631 // result[i,i] += lambda[i] 03632 03633 template<class T> 03634 void addToDiagonal(const TMat<T>& mat, const TVec<T>& lambda) 03635 { 03636 #ifdef BOUNDCHECK 03637 if (lambda.length()!=mat.length()) 03638 PLERROR("Mat(%d)::addToDiagonal(Vec(%d)) inconsistent lengths", 03639 mat.length(), lambda.length()); 03640 #endif 03641 T *l = lambda.data(); 03642 T *d = mat.data(); 03643 for (int i=0;i<mat.length();i++,d+=mat.mod()+1,l++) *d += *l; 03644 } 03645 03646 03647 template<class T> 03648 void diag(const TMat<T>& mat, const TVec<T>& d) 03649 { 03650 T* d_ = d.data(); 03651 for (int i=0;i<mat.length();i++) 03652 d_[i] = mat(i,i); 03653 } 03654 03655 template<class T> 03656 TVec<T> diag(const TMat<T>& mat) 03657 { 03658 TVec<T> d(mat.length()); 03659 diag(mat, d); 03660 return d; 03661 } 03662 03663 template<class T> 03664 void diagonalOfSquare(const TMat<T>& mat, const TVec<T>& d) 03665 { 03666 T* d_=d.data(); 03667 for (int i=0;i<mat.length();i++) 03668 d_[i]=pownorm(mat(i)); 03669 } 03670 03671 03672 template<class T> 03673 void projectOnOrthogonalSubspace(const TMat<T>& mat, TMat<T> orthonormal_subspace) 03674 { 03675 for (int i=0;i<mat.length();i++) 03676 { 03677 TVec<T> row_i = mat(i); 03678 projectOnOrthogonalSubspace(row_i, orthonormal_subspace); 03679 } 03680 } 03681 03682 03683 template<class T> 03684 void averageAcrossRowsAndColumns(const TMat<T>& mat, TVec<T>& avg_across_rows, TVec<T>& avg_across_columns, bool ignored) 03685 { 03686 avg_across_rows.resize(mat.width()); 03687 avg_across_columns.resize(mat.length()); 03688 avg_across_rows.clear(); 03689 avg_across_columns.clear(); 03690 T* row_i=mat.data(); 03691 for (int i=0;i<mat.length();i++) 03692 { 03693 T& avg_cols_i=avg_across_columns[i]; 03694 T* avg_rows = avg_across_rows.data(); 03695 for (int j=0;j<mat.width();j++) 03696 { 03697 T row_ij=row_i[j]; 03698 avg_cols_i += row_ij; 03699 avg_rows[j] += row_ij; 03700 } 03701 row_i+=mat.mod(); 03702 } 03703 avg_across_rows /= mat.length(); 03704 avg_across_columns /= mat.width(); 03705 } 03706 03707 03708 template<class T> 03709 void addToRows(const TMat<T>& mat, const TVec<T> row, bool ignored) 03710 { 03711 for (int i=0;i<mat.length();i++) 03712 { 03713 TVec<T> row_i = mat(i); 03714 row_i += row; 03715 } 03716 } 03717 03718 03719 template<class T> 03720 void addToColumns(const TMat<T>& mat, const TVec<T> col, bool ignored) 03721 { 03722 T* row_i=mat.data(); 03723 for (int i=0;i<mat.length();i++) 03724 { 03725 T col_i=col[i]; 03726 for (int j=0;j<mat.width();j++) 03727 row_i[j] += col_i; 03728 row_i+=mat.mod(); 03729 } 03730 } 03731 03732 template<class T> 03733 void substractFromRows(const TMat<T>& mat, const TVec<T> row, bool ignored) 03734 { 03735 for (int i=0;i<mat.length();i++) 03736 { 03737 TVec<T> row_i = mat(i); 03738 row_i -= row; 03739 } 03740 } 03741 03742 03743 03744 // Probably bugged!!! 03745 template<class T> 03746 void substractFromColumns(const TMat<T>& mat, const TVec<T> col, bool ignored) 03747 { 03748 T* row_i=mat.data(); 03749 for (int i=0;i<mat.length();i++) 03750 { 03751 T col_i=col[i]; 03752 for (int j=0;j<mat.width();j++) 03753 row_i[j] -= col_i; 03754 row_i+=mat.mod(); 03755 } 03756 } 03757 03758 03759 template<class T> 03760 void addToMat(const TMat<T>& mat, T scalar, bool ignored) 03761 { mat += scalar; } 03762 03763 03764 // -------------- taken and adapted from Mat_maths.cc ------------------ 03765 03766 template<class T> 03767 T sum(const TMat<T>& mat, bool ignore_missing=false) 03768 { 03769 double res = 0.0; 03770 T* m_i = mat.data(); 03771 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03772 { 03773 for(int j=0; j<mat.width(); j++) 03774 { 03775 if (!is_missing(m_i[j])) res += m_i[j]; 03776 else if (!ignore_missing) return MISSING_VALUE; 03777 } 03778 } 03779 return T(res); 03780 } 03781 03782 template<class T> 03783 T product(const TMat<T>& mat) 03784 { 03785 double res = 1.0; 03786 T* m_i = mat.data(); 03787 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03788 for(int j=0; j<mat.width(); j++) 03789 res *= m_i[j]; 03790 return T(res); 03791 } 03792 03793 template<class T> 03794 T sum_of_squares(const TMat<T>& mat) 03795 { 03796 double res = 0.0; 03797 T* m_i = mat.data(); 03798 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03799 for(int j=0; j<mat.width(); j++) 03800 { 03801 T v = m_i[j]; 03802 res += v*v; 03803 } 03804 return T(res); 03805 } 03806 03807 template<class T> 03808 T mean(const TMat<T>& mat) 03809 { 03810 #ifdef BOUNDCHECK 03811 if(mat.length()==0 || mat.width()==0) 03812 PLERROR("IN T mean(const TMat<T>& mat) mat has 0 size"); 03813 #endif 03814 double res = 0.0; 03815 T* m_i = mat.data(); 03816 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03817 for(int j=0; j<mat.width(); j++) 03818 res += m_i[j]; 03819 return T(res/(mat.length()*mat.width())); 03820 } 03821 03822 template<class T> 03823 T geometric_mean(const TMat<T>& mat) 03824 { 03825 #ifdef BOUNDCHECK 03826 if(mat.length()==0 || mat.width()==0) 03827 PLERROR("IN T geometric_mean(const TMat<T>& mat) mat has 0 size"); 03828 #endif 03829 double res = 0.0; 03830 T* m_i = mat.data(); 03831 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03832 for(int j=0; j<mat.width(); j++) 03833 { 03834 T mij = m_i[j]; 03835 if (mij<=0) 03836 PLERROR("geometric_mean(TMat<T>): argument %g <=0 at position (%d,%d)", 03837 mij,i,j); 03838 res += log(m_i[j]); 03839 } 03840 return T(exp(res/(mat.length()*mat.width()))); 03841 } 03842 03843 template<class T> 03844 T variance(const TMat<T>& mat, T meanval) 03845 { 03846 #ifdef BOUNDCHECK 03847 if(mat.length()==0 || mat.width()==0) 03848 PLERROR("IN T variance(const TMat<T>& mat, T meanval) mat has 0 size"); 03849 #endif 03850 double res = 0.0; 03851 T* m_i = mat.data(); 03852 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03853 for(int j=0; j<mat.width(); j++) 03854 { 03855 double diff = m_i[j]-meanval; 03856 res += diff*diff; 03857 } 03858 return res/(mat.length()*mat.width()-1); 03859 } 03860 03861 template<class T> 03862 T correlation(const TMat<T>& mat) 03863 { 03864 int n = mat.length(); 03865 #ifdef BOUNDCHECK 03866 if(n==0 || mat.width()==0) 03867 PLERROR("In T correlation(const TMat<T>& mat) mat has 0 size"); 03868 #endif 03869 if (mat.width() != 2) 03870 PLERROR("In T correlation(const TMat<T>& mat), mat width (%d) must be 2", mat.width()); 03871 03872 double s_x=0, s_y=0, s_xy=0, s_x2=0, s_y2=0; 03873 for (int i=0; i<n; i++) 03874 { 03875 T x = mat(i,0); 03876 T y = mat(i,1); 03877 s_x += x; 03878 s_x2 += x*x; 03879 s_y += y; 03880 s_y2 += y*y; 03881 s_xy += x*y; 03882 } 03883 03884 return (n*s_xy - s_x*s_y)/sqrt((n*s_x2 - s_x*s_x)*(n*s_y2 - s_y*s_y)); 03885 } 03886 03887 template<class T> 03888 T correlation(const TVec<T>& x, const TVec<T>& y) 03889 { 03890 int n = x.length(); 03891 #ifdef BOUNDCHECK 03892 if(n==0 || y.length()==0) 03893 PLERROR("In T correlation(const TVec<T>& x, const TVec<T>& y), one Vec has 0 size"); 03894 #endif 03895 if (n != y.length()) 03896 PLERROR("In T correlation(const TVec<T>& x, const TVec<T>& y), both Vec must have same length (%d != %d)", n, y.length()); 03897 03898 double s_x=0, s_y=0, s_xy=0, s_x2=0, s_y2=0; 03899 for (int i=0; i<n; i++) 03900 { 03901 T x_val = x[i]; 03902 T y_val = y[i]; 03903 s_x += x_val; 03904 s_x2 += x_val*x_val; 03905 s_y += y_val; 03906 s_y2 += y_val*y_val; 03907 s_xy += x_val*y_val; 03908 } 03909 03910 return (n*s_xy - s_x*s_y)/sqrt((n*s_x2 - s_x*s_x)*(n*s_y2 - s_y*s_y)); 03911 } 03912 03913 template<class T> 03914 T min(const TMat<T>& mat) 03915 { 03916 #ifdef BOUNDCHECK 03917 if(mat.length()==0 || mat.width()==0) 03918 PLERROR("IN T min(const TMat<T>& mat) mat has 0 size"); 03919 #endif 03920 T* m_i = mat.data(); 03921 double minval = m_i[0]; 03922 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03923 for(int j=0; j<mat.width(); j++) 03924 if(m_i[j]<minval) 03925 minval = m_i[j]; 03926 return minval; 03927 } 03928 03929 template<class T> 03930 T max(const TMat<T>& mat) 03931 { 03932 #ifdef BOUNDCHECK 03933 if(mat.length()==0 || mat.width()==0) 03934 PLERROR("IN T max(const TMat<T>& mat) mat has 0 size"); 03935 #endif 03936 T* m_i = mat.data(); 03937 double maxval = m_i[0]; 03938 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03939 for(int j=0; j<mat.width(); j++) 03940 if(m_i[j]>maxval) 03941 maxval = m_i[j]; 03942 return maxval; 03943 } 03944 03946 template<class T> 03947 void argmin(const TMat<T>& mat, int& mini, int& minj) 03948 { 03949 #ifdef BOUNDCHECK 03950 if(mat.length()==0 || mat.width()==0) 03951 PLERROR("IN void argmin(const TMat<T>& mat, int& mini, iny& minj) mat has 0 size"); 03952 #endif 03953 T* m_i = mat.data(); 03954 mini=0; 03955 minj=0; 03956 double minval = m_i[0]; 03957 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03958 for(int j=0; j<mat.width(); j++) 03959 if(m_i[j]<minval) 03960 { 03961 minval = m_i[j]; 03962 mini = i; 03963 minj = j; 03964 } 03965 } 03966 03967 // Same as above with the max. 03968 template<class T> 03969 void argmax(const TMat<T>& mat, int& maxi, int& maxj) 03970 { 03971 #ifdef BOUNDCHECK 03972 if(mat.length()==0 || mat.width()==0) 03973 PLERROR("IN void argmax(const TMat<T>& mat, int& maxi, iny& maxj) mat has 0 size"); 03974 #endif 03975 T* m_i = mat.data(); 03976 maxi=0; 03977 maxj=0; 03978 double maxval = m_i[0]; 03979 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 03980 for(int j=0; j<mat.width(); j++) 03981 if(m_i[j]>maxval) 03982 { 03983 maxval = m_i[j]; 03984 maxi = i; 03985 maxj = j; 03986 } 03987 } 03988 03990 template<class T> 03991 int argmin(const TMat<T>& m) 03992 { 03993 int imin, jmin; 03994 argmin(m,imin,jmin); 03995 return (imin*m.width()+jmin); 03996 } 03997 03999 template<class T> 04000 int argmax(const TMat<T>& m) 04001 { 04002 int imax, jmax; 04003 argmax(m,imax,jmax); 04004 return (imax*m.width()+jmax); 04005 } 04006 04012 template<class T> 04013 void rowSum(const TMat<T>& mat, const TMat<T>& singlecolumn) 04014 { 04015 #ifdef BOUNDCHECK 04016 if(singlecolumn.length()!=mat.length() || singlecolumn.width() != 1) 04017 PLERROR("IN void rowSum(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix"); 04018 #endif 04019 for(int i=0; i<mat.length(); i++) 04020 singlecolumn(i,0) = sum(mat(i)); 04021 } 04022 04023 04024 template<class T> 04025 void rowSum(const TMat<T>& mat, const TVec<T>& colvec) 04026 { 04027 #ifdef BOUNDCHECK 04028 if(colvec.length()!=mat.length()) 04029 PLERROR("IN void rowSum(const TMat<T>& mat, const TVec<T>& colvec) colvec must have same length as mat"); 04030 #endif 04031 for(int i=0; i<mat.length(); i++) 04032 colvec[i] = sum(mat(i)); 04033 } 04034 04035 template<class T> 04036 void rowMean(const TMat<T>& mat, const TMat<T>& singlecolumn) 04037 { 04038 #ifdef BOUNDCHECK 04039 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 04040 PLERROR("IN void rowMean(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix, and mat must have non-zero width"); 04041 #endif 04042 for(int i=0; i<mat.length(); i++) 04043 singlecolumn(i,0) = mean(mat(i)); 04044 } 04045 04046 template<class T> 04047 void rowVariance(const TMat<T>& mat, const TMat<T>& singlecolumn, const TMat<T>& rowmean) 04048 { 04049 #ifdef BOUNDCHECK 04050 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || rowmean.length()!=mat.length() || rowmean.width()!=1 || mat.width()==0) 04051 PLERROR("IN void rowVariance(const TMat<T>& mat, TMat<T>& singlecolumn, const TMat<T>& rowmean) singlecolumn and rowmean must be mat.length() x 1 matrices, mat must have non-zero width"); 04052 #endif 04053 for(int i=0; i<mat.length(); i++) 04054 singlecolumn(i,0) = variance(mat(i),rowmean(i,0)); 04055 } 04056 04057 template<class T> 04058 void rowSumOfSquares(const TMat<T>& mat, const TMat<T>& singlecolumn) 04059 { 04060 #ifdef BOUNDCHECK 04061 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1) 04062 PLERROR("IN void rowSumOfSquares(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix"); 04063 #endif 04064 for (int i=0;i<mat.length();i++) 04065 { 04066 T ss=0; 04067 T* mi=mat[i]; 04068 for (int j=0;j<mat.width();j++) { T mij=mi[j]; ss+=mij*mij; } 04069 singlecolumn(i,0)=ss; 04070 } 04071 } 04072 04073 template<class T> 04074 void rowMax(const TMat<T>& mat, const TMat<T>& singlecolumn) 04075 { 04076 #ifdef BOUNDCHECK 04077 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 04078 PLERROR("IN void rowMax(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix, and mat must have non-zero width"); 04079 #endif 04080 for(int i=0; i<mat.length(); i++) 04081 singlecolumn(i,0) = max(mat(i)); 04082 } 04083 04084 template<class T> 04085 void rowMax(const TMat<T>& mat, const TVec<T>& colvec) 04086 { 04087 #ifdef BOUNDCHECK 04088 if(colvec.length()!=mat.length()) 04089 PLERROR("IN void rowSum(const TMat<T>& mat, const TVec<T>& colvec) colvec must have same length as mat"); 04090 #endif 04091 for(int i=0; i<mat.length(); i++) 04092 colvec[i] = max(mat(i)); 04093 } 04094 04095 template<class T> 04096 void rowMin(const TMat<T>& mat, const TMat<T>& singlecolumn) 04097 { 04098 #ifdef BOUNDCHECK 04099 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 04100 PLERROR("IN void rowMin(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix, and mat must have non-zero width"); 04101 #endif 04102 for(int i=0; i<mat.length(); i++) 04103 singlecolumn(i,0) = min(mat(i)); 04104 } 04105 04106 04107 template<class T> 04108 void rowMin(const TMat<T>& mat, const TVec<T>& colvec) 04109 { 04110 #ifdef BOUNDCHECK 04111 if(colvec.length()!=mat.length()) 04112 PLERROR("IN void rowSum(const TMat<T>& mat, const TVec<T>& colvec) colvec must have same length as mat"); 04113 #endif 04114 for(int i=0; i<mat.length(); i++) 04115 colvec[i] = min(mat(i)); 04116 } 04117 04118 template<class T> 04119 void rowArgmax(const TMat<T>& mat, const TMat<T>& singlecolumn) 04120 { 04121 #ifdef BOUNDCHECK 04122 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 04123 PLERROR("IN void rowMax(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix, and mat must have non-zero width"); 04124 #endif 04125 for(int i=0; i<mat.length(); i++) 04126 singlecolumn(i,0) = argmax(mat(i)); 04127 } 04128 04129 template<class T> 04130 void rowArgmin(const TMat<T>& mat, const TMat<T>& singlecolumn) 04131 { 04132 #ifdef BOUNDCHECK 04133 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 04134 PLERROR("IN void rowMax(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix, and mat must have non-zero width"); 04135 #endif 04136 for(int i=0; i<mat.length(); i++) 04137 singlecolumn(i,0) = argmin(mat(i)); 04138 } 04139 04145 template<class T> 04146 void columnSum(const TMat<T>& mat, TVec<T>& result) 04147 { 04148 #ifdef BOUNDCHECK 04149 if(result.length()!=mat.width()) 04150 PLERROR("IN void columnSum(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat"); 04151 #endif 04152 int l = mat.length(); 04153 result << mat(0); 04154 for(int j=1; j<l; j++) 04155 result += mat(j); 04156 } 04157 04158 template<class T> 04159 void columnSumOfSquares(const TMat<T>& mat, TVec<T>& result) 04160 { 04161 #ifdef BOUNDCHECK 04162 if(result.length()!=mat.width()) 04163 PLERROR("IN void columnSumOfSquares(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat"); 04164 #endif 04165 for(int j=0; j<mat.width(); j++) 04166 result[j] = sum_of_squares(mat.column(j)); 04167 } 04168 04169 template<class T> 04170 void columnMean(const TMat<T>& mat, TVec<T>& result) 04171 { 04172 #ifdef BOUNDCHECK 04173 if(result.length()!=mat.width() || mat.length()==0) 04174 PLERROR("IN void columnMean(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat and mat must have non-zero length"); 04175 #endif 04176 columnSum(mat,result); 04177 result /= real(mat.length()); 04178 } 04179 04180 template<class T> 04181 void columnWeightedMean(const TMat<T>& mat, TVec<T>& result) 04182 { 04183 #ifdef BOUNDCHECK 04184 if(result.length()!=mat.width()-1 || mat.length()<=1) 04185 PLERROR("IN void columnWeightedMean(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width - 1 of mat and mat must have at least 1 length"); 04186 #endif 04187 TVec<T> column_j_vec(mat.length()), weights_vec(mat.length()); 04188 TMat<T> column_j_mat(mat.length(), 1), weights_mat(mat.length(), 1); 04189 for(int j=0; j<mat.width()-1; j++){ 04190 column_j_mat = mat.column(j); 04191 weights_mat = mat.column(mat.width()-1); 04192 column_j_vec = column_j_mat.toVecCopy(); 04193 weights_vec = weights_mat.toVecCopy(); 04194 result[j] = weighted_mean(column_j_vec, weights_vec); 04195 } 04196 } 04197 04198 template<class T> 04199 void columnVariance(const TMat<T>& mat, TVec<T>& result, const TVec<T>& columnmean) 04200 { 04201 #ifdef BOUNDCHECK 04202 if(result.length()!=mat.width() || columnmean.length()!=mat.width() || mat.length()==0) 04203 PLERROR("IN void columnVariance(const TMat<T>& mat, TVec<T>& result, const TVec<T>& columnmean) the length of result and columnmean must equal the width of mat and mat must have non-zero length"); 04204 #endif 04205 for(int j=0; j<mat.width(); j++) 04206 result[j] = variance(mat.column(j),columnmean[j]); 04207 } 04208 04209 template<class T> 04210 void columnWeightedVariance(const TMat<T>& mat, TVec<T>& result, const TVec<T>& column_weighted_mean) 04211 { 04212 #ifdef BOUNDCHECK 04213 if(result.length()!=mat.width()-1 || column_weighted_mean.length()!=mat.width()-1 || mat.length()<=1) 04214 PLERROR("IN void columnWeightedVariance(const TMat<T>& mat, TVec<T>& result, const TVec<T>& column_weighted_mean) the length of result and column_weighted_mean must equal the width - 1 of mat and mat must have at least 1 length"); 04215 #endif 04216 T column_no_weighted_mean_j; 04217 TVec<T> column_j_vec(mat.length()), weights_vec(mat.length()); 04218 TMat<T> column_j_mat(mat.length(), 1), weights_mat(mat.length(), 1); 04219 for(int j=0; j<mat.width()-1; j++){ 04220 column_j_mat = mat.column(j); 04221 weights_mat = mat.column(mat.width()-1); 04222 column_j_vec = column_j_mat.toVecCopy(); 04223 weights_vec = weights_mat.toVecCopy(); 04224 column_no_weighted_mean_j = mean(mat.column(j)); 04225 result[j] = weighted_variance(column_j_vec, weights_vec, column_no_weighted_mean_j, column_weighted_mean[j]); 04226 } 04227 } 04228 04229 template<class T> 04230 void columnMax(const TMat<T>& mat, TVec<T>& result) 04231 { 04232 #ifdef BOUNDCHECK 04233 if(result.length()!=mat.width() || mat.length()==0) 04234 PLERROR("IN void columnMax(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat and mat must have non-zero length"); 04235 #endif 04236 for(int j=0; j<mat.width(); j++) 04237 result[j] = max(mat.column(j)); 04238 } 04239 04240 template<class T> 04241 void columnMin(const TMat<T>& mat, TVec<T>& result) 04242 { 04243 #ifdef BOUNDCHECK 04244 if(result.length()!=mat.width() || mat.length()==0) 04245 PLERROR("IN void columnMax(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat and mat must have non-zero length"); 04246 #endif 04247 for(int j=0; j<mat.width(); j++) 04248 result[j] = min(mat.column(j)); 04249 } 04250 04251 template<class T> 04252 void columnArgmax(const TMat<T>& mat, TVec<T>& result) 04253 { 04254 #ifdef BOUNDCHECK 04255 if(result.length()!=mat.width() || mat.length()==0) 04256 PLERROR("IN void columnMax(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat and mat must have non-zero length"); 04257 #endif 04258 int imax, jmax; 04259 for(int j=0; j<mat.width(); j++) 04260 { 04261 argmax(mat.column(j), imax, jmax); 04262 result[j] = imax; 04263 } 04264 } 04265 04266 template<class T> 04267 void columnArgmin(const TMat<T>& mat, TVec<T>& result) 04268 { 04269 #ifdef BOUNDCHECK 04270 if(result.length()!=mat.width() || mat.length()==0) 04271 PLERROR("IN void columnMax(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat and mat must have non-zero length"); 04272 #endif 04273 int imin, jmin; 04274 for(int j=0; j<mat.width(); j++) 04275 { 04276 argmin(mat.column(j), imin, jmin); 04277 result[j] = imin; 04278 } 04279 } 04280 04281 template<class T> 04282 T mahalanobis_distance(const TVec<T>& input, const TVec<T>& meanvec, const 04283 TMat<T>& inversecovmat) 04284 { 04285 TVec<T> diff = input-meanvec; 04286 return dot(diff,product(inversecovmat,diff)); 04287 } 04288 04290 template<class T> 04291 inline void computeMean(const TMat<T>& m, TVec<T>& meanvec) { columnMean(m,meanvec); } 04292 04294 template<class T> 04295 void computeMeanAndVariance(const TMat<T>& m, TVec<T>& meanvec, TVec<T>& variancevec) 04296 { 04297 columnMean(m,meanvec); 04298 columnVariance(m,variancevec,meanvec); 04299 } 04300 04301 template<class T> 04302 void computeCovar(const TMat<T>& m, const TVec<T>& meanvec, TMat<T>& covarmat) 04303 { 04304 int n = m.width(); 04305 covarmat.resize(n,n); 04306 transposeProduct(covarmat,m,m); 04307 covarmat /= T(m.length()); 04308 externalProductScaleAcc(covarmat,meanvec,meanvec,T(-1)); 04309 } 04310 04311 template<class T> 04312 void computeMeanAndCovar(const TMat<T>& m, TVec<T>& meanvec, TMat<T>& covarmat) 04313 { 04314 int n = m.width(); 04315 meanvec.resize(n); 04316 covarmat.resize(n,n); 04317 columnMean(m,meanvec); 04318 04319 transposeProduct(covarmat,m,m); 04320 covarmat /= T(m.length()); 04321 externalProductScaleAcc(covarmat,meanvec,meanvec,T(-1)); 04322 04323 /* 04324 Mat mm = m.copy(); 04325 mm -= meanvec; 04326 transposeProduct(covarmat,mm,mm); 04327 covarmat /= T(m.length()); 04328 */ 04329 } 04330 04332 template<class T> 04333 void computeMeanAndStddev(const TMat<T>& m, TVec<T>& meanvec, TVec<T>& stddevvec) 04334 { 04335 columnMean(m,meanvec); 04336 columnVariance(m,stddevvec,meanvec); 04337 for(int i=0; i<stddevvec.length(); i++) 04338 stddevvec[i] = sqrt(stddevvec[i]); 04339 } 04340 04341 04344 template<class T> 04345 void computeColumnsMeanAndStddev(const TMat<T>& m, TMat<T>& meanvec, TMat<T>& stddevvec) 04346 { 04347 rowMean(m,meanvec); 04348 rowVariance(m,stddevvec,meanvec); 04349 for(int i=0; i<stddevvec.length(); i++) 04350 stddevvec[i][0] = sqrt(stddevvec[i][0]); 04351 } 04352 04354 template<class T> 04355 void normalize(TMat<T>& m) 04356 { 04357 TVec<T> meanvec(m.width()); 04358 TVec<T> stddevvec(m.width()); 04359 computeMeanAndStddev(m,meanvec,stddevvec); 04360 m -= meanvec; 04361 m /= stddevvec; 04362 } 04363 04365 template<class T> 04366 void normalizeRows(const TMat<T>& m) 04367 { 04368 int l = m.length(); 04369 for(int i=0; i<l; i++) 04370 { 04371 TVec<T> v = m(i); 04372 v /= sum(v); 04373 } 04374 } 04375 04377 template<class T> 04378 void normalizeColumns(const TMat<T>& m) 04379 { 04380 int w = m.width(); 04381 for(int j=0; j<w; j++) 04382 { 04383 TMat<T> v = m.column(j); 04384 v /= sum(v); 04385 } 04386 } 04387 04389 template<class T> 04390 void normalize(TMat<T>& m, double n) 04391 { 04392 for(int i=0; i<m.length(); i++) 04393 { 04394 TVec<T> m_i = m(i); 04395 normalize(m_i,n); 04396 } 04397 } 04398 04399 template<class T> 04400 void operator+=(const TMat<T>& m, T scalar) 04401 { 04402 T* m_i = m.data(); 04403 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 04404 for(int j=0; j<m.width(); j++) 04405 m_i[j] += scalar; 04406 } 04407 04408 template<class T> 04409 void operator*=(const TMat<T>& m, T scalar) 04410 { 04411 T* m_i = m.data(); 04412 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 04413 for(int j=0; j<m.width(); j++) 04414 m_i[j] *= scalar; 04415 } 04416 04417 template<class T> 04418 inline void operator-=(const TMat<T>& m, T scalar) { m += (-scalar); } 04419 04420 template<class T> 04421 inline void operator/=(const TMat<T>& m, T scalar) { m *= (T(1)/scalar); } 04422 04423 template<class T> 04424 inline void operator/=(const TMat<T>& m, int scalar) { m *= (T(1)/scalar); } 04425 04426 04428 template<class T> 04429 void operator+=(const TMat<T>& m, const TVec<T>& v) 04430 { 04431 #ifdef BOUNDCHECK 04432 if(m.width()!=v.length()) 04433 PLERROR("IN operator+=(const TMat<T>& m, const TVec<T>& v) v must be as long as m is wide"); 04434 #endif 04435 T* m_i = m.data(); 04436 T* vv = v.data(); 04437 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 04438 for(int j=0; j<m.width(); j++) 04439 m_i[j] += vv[j]; 04440 } 04441 04443 template<class T> 04444 void operator-=(const TMat<T>& m, const TVec<T>& v) 04445 { 04446 #ifdef BOUNDCHECK 04447 if(m.width()!=v.length()) 04448 PLERROR("IN operator-=(const TMat<T>& m, const TVec<T>& v) v must be as long as m is wide"); 04449 #endif 04450 T* m_i = m.data(); 04451 T* vv = v.data(); 04452 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 04453 for(int j=0; j<m.width(); j++) 04454 m_i[j] -= vv[j]; 04455 } 04456 04458 template<class T> 04459 void operator*=(const TMat<T>& m, const TVec<T>& v) 04460 { 04461 #ifdef BOUNDCHECK 04462 if(m.width()!=v.length()) 04463 PLERROR("IN operator*=(const TMat<T>& m, const TVec<T>& v) v must be as long as m is wide"); 04464 #endif 04465 T* m_i = m.data(); 04466 T* vv = v.data(); 04467 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 04468 for(int j=0; j<m.width(); j++) 04469 m_i[j] *= vv[j]; 04470 } 04471 04473 template<class T> 04474 void operator*=(const TMat<T>& m1, const TMat<T>& m2) 04475 { 04476 int n=m1.length(); 04477 int l=m1.width(); 04478 #ifdef BOUNDCHECK 04479 if(l!=m2.width() || n!=m2.length()) 04480 PLERROR("IN operator*=(const TMat<T>& m1(%d,%d), const TMat<T>& m2(%d,%d)) sizes differ", 04481 m1.length(),m1.width(),m2.length(),m2.width()); 04482 #endif 04483 T* m1_i = m1.data(); 04484 T* m2_i = m2.data(); 04485 for(int i=0; i<n; i++, m1_i+=m1.mod(),m2_i+=m2.mod()) 04486 for(int j=0; j<l; j++) 04487 m1_i[j] *= m2_i[j]; 04488 } 04489 04491 template<class T> 04492 void operator/=(const TMat<T>& m, const TVec<T>& v) 04493 { 04494 #ifdef BOUNDCHECK 04495 if(m.width()!=v.length()) 04496 PLERROR("IN operator/=(const TMat<T>& m, const TVec<T>& v) v must be as long as m is wide"); 04497 #endif 04498 T* m_i = m.data(); 04499 T* vv = v.data(); 04500 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 04501 for(int j=0; j<m.width(); j++) 04502 m_i[j] /= vv[j]; 04503 } 04504 04506 template<class T> 04507 void operator/=(const TMat<T>& m1, const TMat<T>& m2) 04508 { 04509 int n=m1.length(); 04510 int l=m1.width(); 04511 #ifdef BOUNDCHECK 04512 if(l!=m2.width() || n!=m2.length()) 04513 PLERROR("IN operator/=(const TMat<T>& m1(%d,%d), const TMat<T>& m2(%d,%d)) sizes differ", 04514 m1.length(),m1.width(),m2.length(),m2.width()); 04515 #endif 04516 T* m1_i = m1.data(); 04517 T* m2_i = m2.data(); 04518 for(int i=0; i<n; i++, m1_i+=m1.mod(),m2_i+=m2.mod()) 04519 for(int j=0; j<l; j++) 04520 m1_i[j] /= m2_i[j]; 04521 } 04522 04523 template<class T> 04524 void operator+=(const TMat<T>& m1, const TMat<T>& m2) 04525 { 04526 int n=m1.length(); 04527 int l=m1.width(); 04528 #ifdef BOUNDCHECK 04529 if(m1.width()!=m2.width() || m1.length()!=m2.length()) 04530 PLERROR("IN operator+=(const TMat<T>& m1(%d,%d), const TMat<T>& m2(%d,%d)): m1 and m2 must have same dimensions", 04531 m1.length(),m1.width(),m2.length(),m2.width()); 04532 #endif 04533 T* m1_i = m1.data(); 04534 T* m2_i = m2.data(); 04535 for(int i=0; i<n; i++, m1_i+=m1.mod(),m2_i+=m2.mod()) 04536 for(int j=0; j<l; j++) 04537 m1_i[j] += m2_i[j]; 04538 } 04539 04540 template<class T> 04541 void operator-=(const TMat<T>& m1, const TMat<T>& m2) 04542 { 04543 int n=m1.length(); 04544 int l=m1.width(); 04545 #ifdef BOUNDCHECK 04546 if(m1.width()!=m2.width() || m1.length()!=m2.length()) 04547 PLERROR("IN operator+=(const TMat<T>& m1(%d,%d), const TMat<T>& m2(%d,%d)): m1 and m2 must have same dimensions", 04548 m1.length(),m1.width(),m2.length(),m2.width()); 04549 #endif 04550 T* m1_i = m1.data(); 04551 T* m2_i = m2.data(); 04552 for(int i=0; i<n; i++, m1_i+=m1.mod(),m2_i+=m2.mod()) 04553 for(int j=0; j<l; j++) 04554 m1_i[j] -= m2_i[j]; 04555 } 04556 04557 template<class T> 04558 TMat<T> operator-(const TMat<T>& m1, const TMat<T>& m2) 04559 { 04560 TMat<T> result(m1.length(), m1.width()); 04561 substract(m1,m2,result); 04562 return result; 04563 } 04564 04565 template<class T> 04566 TMat<T> operator+(const TMat<T>& m1, const TMat<T>& m2) 04567 { 04568 TMat<T> result(m1.length(), m1.width()); 04569 add(m1,m2,result); 04570 return result; 04571 } 04572 04573 template<class T> 04574 void substract(const TMat<T>& m1, const TMat<T>& m2, TMat<T>& destination) 04575 { 04576 #ifdef BOUNDCHECK 04577 if(m1.width()!=m2.width() || m1.length()!=m2.length() 04578 || m1.width()!=destination.width() || m1.length()!=destination.length()) 04579 PLERROR("IN substract(m1(%d,%d), m2(%d,%d), dest(%d,%d)): args must have same dimensions", 04580 m1.length(),m1.width(),m2.length(),m2.width(),destination.length(), 04581 destination.width()); 04582 #endif 04583 T* m1_i = m1.data(); 04584 T* m2_i = m2.data(); 04585 T* d_i = destination.data(); 04586 int m1_mod = m1.mod(); 04587 int m2_mod = m2.mod(); 04588 int d_mod = destination.mod(); 04589 int w = m1.width(); 04590 for (int i=0;i<m1.length();i++,m1_i+=m1_mod,m2_i+=m2_mod,d_i+=d_mod) 04591 for (int j=0;j<w;j++) 04592 d_i[j] = m1_i[j] - m2_i[j]; 04593 } 04594 04595 template<class T> 04596 void add(const TMat<T>& m1, const TMat<T>& m2, TMat<T>& destination) 04597 { 04598 #ifdef BOUNDCHECK 04599 if(m1.width()!=m2.width() || m1.length()!=m2.length() 04600 || m1.width()!=destination.width() || m1.length()!=destination.length()) 04601 PLERROR("IN substract(m1(%d,%d), m2(%d,%d), dest(%d,%d)): args must have same dimensions", 04602 m1.length(),m1.width(),m2.length(),m2.width(),destination.length(), 04603 destination.width()); 04604 #endif 04605 T* m1_i = m1.data(); 04606 T* m2_i = m2.data(); 04607 T* d_i = destination.data(); 04608 int m1_mod = m1.mod(); 04609 int m2_mod = m2.mod(); 04610 int d_mod = destination.mod(); 04611 int w = m1.width(); 04612 for (int i=0;i<m1.length();i++,m1_i+=m1_mod,m2_i+=m2_mod,d_i+=d_mod) 04613 for (int j=0;j<w;j++) 04614 d_i[j] = m1_i[j] + m2_i[j]; 04615 } 04616 04618 template<class T> 04619 TMat<T> operator-(const TMat<T>& m) 04620 { 04621 TMat<T> opposite(m.length(),m.width()); 04622 T *m_i=m.data(); 04623 T *o_i=opposite.data(); 04624 for (int i=0;i<m.length();i++,m_i+=m.mod(),o_i+=opposite.mod()) 04625 for (int j=0;j<m.width();j++) 04626 o_i[j] = - m_i[j]; 04627 return opposite; 04628 } 04629 04631 template<class T> 04632 void negateElements(const TMat<T>& m) 04633 { 04634 T* m_i = m.data(); 04635 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 04636 for(int j=0; j<m.width(); j++) 04637 m_i[j] = -m_i[j]; 04638 } 04639 04641 template<class T> 04642 void invertElements(const TMat<T>& m) 04643 { 04644 T* m_i = m.data(); 04645 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 04646 for(int j=0; j<m.width(); j++) 04647 m_i[j] = 1.0/m_i[j]; 04648 } 04649 04650 // result * m = identity 04651 // (works only if m.length()>=m.width()) 04652 template<class T> 04653 TMat<T> leftPseudoInverse(TMat<T>& m) 04654 { 04655 TMat<T> inv(m.width(), m.length()); 04656 leftPseudoInverse(m,inv); 04657 return inv; 04658 } 04659 04660 // result * m = identity 04661 // (works only if m.length()>=m.width()) 04662 template<class T> 04663 void leftPseudoInverse(const TMat<T>& m, TMat<T>& inv) 04664 { 04665 if (m.length()==m.width()) 04666 inverse(m,inv); 04667 if (m.length()<m.width()) 04668 PLERROR("leftPseudoInverse: matrix length(%d) must be >= width(%d)", 04669 m.length(), m.width()); 04670 PLERROR("SVD not implemented yet"); 04671 } 04672 04673 // m * result = identity 04674 // (works only if m.length()<=m.width()) 04675 template<class T> 04676 TMat<T> rightPseudoInverse(TMat<T>& m) 04677 { 04678 TMat<T> inv(m.width(), m.length()); 04679 rightPseudoInverse(m,inv); 04680 return inv; 04681 } 04682 04683 // m * result = identity 04684 // (works only if m.length()<=m.width()) 04685 template<class T> 04686 void rightPseudoInverse(const TMat<T>& m, TMat<T>& inv) 04687 { 04688 if (m.length()==m.width()) 04689 inverse(m,inv); 04690 if (m.length()>m.width()) 04691 PLERROR("rightPseudoInverse: matrix length(%d) must be <= width(%d)", 04692 m.length(), m.width()); 04693 PLERROR("SVD not implemented yet"); 04694 } 04695 04696 // find and return inv s.t. m * inv = inv * m = I = identity 04697 // (m must be square) 04698 template<class T> 04699 TMat<T> inverse(TMat<T>& m) 04700 { 04701 TMat<T> inv(m.length(),m.length()); 04702 inverse(m,inv); 04703 return inv; 04704 } 04705 04706 // find inv s.t. m * inv = inv * m = I = identity 04707 // (m must be square) 04708 template<class T> 04709 void inverse(const TMat<T>& m, TMat<T>& inv) 04710 { 04711 int n=m.length(); 04712 if (m.width()!=n) 04713 PLERROR("inverse(TMat<T>,TMat<T>): argument(%d,%d) must be square matrix", 04714 m.width(), n); 04715 inv.resize(n,n); 04716 if (n==1) 04717 inv.data()[0]=1.0/m.data()[0]; 04718 else 04719 PLERROR("matrix inverse not implemented yet"); 04720 } 04721 04722 // for square positive definite symmetric matrices A, 04723 // find X(n,m) s.t. A(n,n) X(n,m) = B(n,m). 04724 // This is obtained by doing a Cholesky decomposition 04725 // A = L L', with L lower diagonal, thus to solve 04726 // L L' X = B. 04727 // We use the CholeskySolve function which solves for x_i in L L' x_i = b_i 04728 // (on the columns x_i and b_i of X and B respectively). 04729 // Optionally provide pointers to the temporary matrix L(n,n) and vector y(n) 04730 // to avoid memory allocations. 04731 template<class T> 04732 void solveLinearSystemByCholesky(const TMat<T>& A, const TMat<T>& B, TMat<T>& X, TMat<T>* pL=0, TVec<T>* py=0) 04733 { 04734 int n=A.length(); 04735 int m=X.width(); 04736 if (X.length()!=n || A.width()!=n || B.length()!=n || B.width()!=m) 04737 PLERROR("solveLinearSystemByCholesky: A(%d,%d) * X(%d,%d) == B(%d,%d), incompatible", 04738 n,A.width(),X.length(),m,B.length(),B.width()); 04739 TMat<T>* L; 04740 TVec<T>* y; 04741 if (pL) L=pL; else L = new TMat<T>(n,n); 04742 if (py) y=py; else y = new TVec<T>(n); 04743 choleskyDecomposition(A,*L); 04744 choleskySolve(*L,B,X,*y); 04745 if (!pL) delete L; 04746 if (!py) delete y; 04747 } 04748 04749 // for square positive definite symmetric matrices A, 04750 // find X(n,m) s.t. X(n,m) A(m,m) = B(n,m). 04751 // This is obtained by doing a Cholesky decomposition 04752 // A = L L', with L lower diagonal, thus to solve 04753 // X L L' = B. 04754 // We use the CholeskySolve function which solves for x_i in L L' x_i = b_i: 04755 // L L' X' = B' 04756 // is solved on the rows of X (x_i) and the columns of B (b_i). 04757 // Optionally provide pointers to the temporary matrices L and y 04758 // to avoid memory allocations. 04759 template<class T> 04760 void solveTransposeLinearSystemByCholesky(const TMat<T>& A, const TMat<T>& B, TMat<T>& X,TMat<T>* pL=0, TVec<T>* py=0) 04761 { 04762 int n=X.length(); 04763 int m=X.width(); 04764 if (A.length()!=m || A.width()!=m || B.length()!=n || B.width()!=m) 04765 PLERROR("solveTransposeLinearSystemByCholesky: X(%d,%d) * A(%d,%d) == B(%d,%d), incompatible", 04766 n,m,A.length(),A.width(),B.length(),B.width()); 04767 TMat<T>* L; 04768 TVec<T>* y; 04769 if (pL) L=pL; else L = new TMat<T>(m,m); 04770 if (py) y=py; else y = new TVec<T>(m); 04771 choleskyDecomposition(A,*L); 04772 for (int i=0;i<n;i++) 04773 choleskySolve(*L,B(i),X(i),*y); 04774 if (!pL) delete L; 04775 if (!py) delete y; 04776 } 04777 04778 /* Perform a Cholesky decomposition of nxn symmetric positive definite 04779 matrix A, i.e., decompose it into 04780 A = L L' 04781 where L is a lower diagonal matrix (with zeros above the diagonal). 04782 L be used to solve a linear system A x = b, i.e., LL'x=b, with choleskySolve(L,b,x). 04783 See choleskySolve(TMat<T>*,TVec<T>*) for an example of use. 04784 04785 From the above equation, one obtains 04786 04787 for i=0..n-1 04788 L[i][i] = sqrt(A[i][i] - sum_{k=0}^{i-1} L[i][k]^2) 04789 for j=i+1... n-1 04790 L[j][i] = (1/L[i][i]) ( A[i][j] - sum_{k=0}^{i-1} L[i][k] L[j][k] ) 04791 04792 */ 04793 template<class T> 04794 void choleskyDecomposition(const TMat<T>& A, TMat<T>& L) 04795 { 04796 int n = A.length(); 04797 if (n!=A.width()) 04798 PLERROR("choleskyDecomposition: non-square matrix %dx%d\n",n,A.width()); 04799 L.resize(n,n); 04800 int i,j,k; 04801 T sum; 04802 bool restart=false; 04803 do 04804 { 04805 restart=false; 04806 for (i=0;i<n;i++) 04807 { 04808 const T* Ai = A[i]; 04809 T* Li = L[i]; 04810 T Lii=0; 04811 for (j=i;j<n;j++) 04812 { 04813 T* Lj = L[j]; 04814 for (sum=Ai[j],k=i-1;k>=0;k--) sum -= Li[k] * Lj[k]; 04815 if (i==j) 04816 { 04817 if (sum <= 0.0) 04818 { 04819 T eps = -1.1*sum; 04820 if (sum==0) eps=1e-8; 04821 PLWARNING("Cholesky decomposition would fail: add %g to diagonal",eps); 04822 T* Aii=A.data(); 04823 int addm=A.mod()+1; 04824 for (int ii=0;ii<n;ii++,Aii+=addm) *Aii += eps; 04825 restart=true; 04826 break; 04827 } 04828 Lii = sqrt(sum); 04829 } 04830 else Lj[i] = sum/Lii; 04831 } 04832 if (restart) break; 04833 Li[i] = Lii; 04834 } 04835 } 04836 while (restart); 04837 04838 } 04839 04840 /* Back-propagate through the call to choleskyDecomposition(A,L). 04841 The argument A holds the original symmetric positive definite 04842 matrix while is the lower diagonal matrix with L L' = A. 04843 Given the derivative of C wrt L, fill the derivative 04844 of C wrt A. dC_dA must have been cleared beforehand. 04845 We are given A, L, dC_dL, and write into dC_dA. 04846 Note that dC_dL is modified after the call 04847 because of the internal dependencies between the L's. 04848 04849 for i=n-1..0 04850 for j=n-1..i+1 04851 dC_dL[i][i] -= dC_dL[j][i] L[j][i] / L[i][i] 04852 dC_dA[i][j] += dC_dL[j][i] / L[i][i] 04853 for k=0..i-1 04854 dC_dL[i][k] -= dC_dL[j][i] L[j][k] / L[i][i] 04855 dC_dL[j][k] -= dC_dL[j][i] L[i][k] / L[i][i] 04856 dC_dA[i][i] += 0.5 * dC_dL[i][i] / L[i][i] 04857 for k=0..i-1 04858 dC_dL[i][k] -= dC_dL[i][i] L[i][k] / L[i][i] 04859 04860 */ 04861 template<class T> 04862 void bpropCholeskyDecomposition(const TMat<T>& A, const TMat<T>& L, 04863 TMat<T>& dC_dA, TMat<T>& dC_dL) 04864 { 04865 int n = A.length(); 04866 if (dC_dA) 04867 dC_dA.resize(n,n); 04868 int i,j,k; 04869 for (i=n-1;i>=0;i--) 04870 { 04871 const T* Li = L[i]; 04872 T* dC_dLi = dC_dL[i]; 04873 T* dC_dAi = dC_dA[i]; 04874 T invLii = 1.0/Li[i]; 04875 for (j=n-1;j>i;j--) 04876 { 04877 const T* Lj = L[j]; 04878 T* dC_dLj = dC_dL[j]; 04879 T dC_dLji = dC_dLj[i]; 04880 dC_dLi[i] -= dC_dLji * Lj[i] * invLii; 04881 dC_dAi[j] += dC_dLji * invLii; 04882 for (k=0;k<i;k++) 04883 { 04884 dC_dLi[k] -= dC_dLji * Lj[k] * invLii; 04885 dC_dLj[k] -= dC_dLji * Li[k] * invLii; 04886 } 04887 } 04888 T dC_dLii = dC_dLi[i]; 04889 dC_dAi[i] += 0.5 * dC_dLii * invLii; 04890 for (k=0;k<i;k++) 04891 dC_dLi[k] -= dC_dLii * Li[k] * invLii; 04892 } 04893 } 04894 04895 /* Solve the linear system A x = L L' x = b using a Cholesky decomposition 04896 of A into L L' performed with a prior call to choleskyDecomposition(A,L) 04897 (which on return has the matrix L, that is lower diagonal, and A = L L'). 04898 The solution of the linear system L L' x = b will be in x. 04899 See choleskySolve(TMat<T>*,TVec<T>*) for an example of use. 04900 The algorithm is first to solve L y = b, and then L' x = y. 04901 The argument y is optional and can be used to hold the intermediate 04902 solution to L y = b. 04903 04904 The solution to L L' x = b is obtained as follows: 04905 04906 * Solve L y = b by iterating once through the rows of L 04907 (store result in x): 04908 y[i] = (b[i] - sum_{k<i} L[i][k] y[k])/L[i][i] 04909 04910 * Solve L' x = y by iterating once (backwards) through the rows of L. 04911 x[i] = (y[i] - sum_{k>i} L[k][i] x[k])/L[i][i] 04912 04913 */ 04914 template<class T> 04915 void choleskySolve(const TMat<T>& L, TVec<T> b, TVec<T> x, TVec<T>& y) 04916 { 04917 int i,k; 04918 T sum; 04919 int n = L.length(); 04920 if (L.width()!=n) 04921 PLERROR("choleskySolve: matrix L (%d x %d) is not square!", 04922 n, L.width()); 04923 if (b.length()!=n || x.length()!=n) 04924 PLERROR("choleskySolve: RHS vector b(%d) or unknown x(%d) incompatiable with L(%d,%d)", 04925 b.length(),x.length(),n,n); 04926 04927 T* bp = b.data(); 04928 T* xp = x.data(); 04929 T* yp = y.data(); 04930 04931 // solve L y = b (in variable x if y=0): 04932 // for i=0..n-1 04933 // y[i] = (b[i] - sum_{k<i} L[i][k] y[k])/L[i][i] 04934 for (i=0;i<n;i++) 04935 { 04936 const T* Li = L[i]; 04937 for (sum=bp[i],k=i-1;k>=0;k--) sum -= Li[k] * yp[k]; 04938 yp[i] = sum / Li[i]; 04939 } 04940 // solve L' x = y 04941 // for i=n-1..0 04942 // x[i] = (y[i] - sum_{k>i} L[k][i] x[k])/L[i][i] 04943 for (i=n-1;i>=0;i--) 04944 { 04945 for (sum=yp[i],k=i+1;k<n;k++) sum -= L[k][i] * xp[k]; 04946 xp[i] = sum / L[i][i]; 04947 } 04948 } 04949 04950 // same as the previous choleskySolve but do it m times on the columns 04951 // of nxm matrices X and B. 04952 template<class T> 04953 void choleskySolve(const TMat<T>& L, const TMat<T>& B, TMat<T>& X, TVec<T>& y) 04954 { 04955 int i,k; 04956 T sum; 04957 int n = L.length(); 04958 int m = X.width(); 04959 if (L.width()!=n) 04960 PLERROR("choleskySolve: matrix L (%d x %d) is not square!", 04961 n, L.width()); 04962 if (B.length()!=n || B.width() !=m) 04963 PLERROR("choleskySolve: RHS matrix B(%d,%d) instead of (%d,%d) like X", 04964 B.length(),B.width(), n, m); 04965 if (X.length()!=n) 04966 PLERROR("choleskySolve: X(%d,%d) not compatible with L(%d,%d)", 04967 X.length(),m,n,n); 04968 if (y.length()!=n) 04969 PLERROR("choleskySolve: y(%d) not compatible with L(%d,%d)", 04970 y.length(),n,n); 04971 int bmod = B.mod(); 04972 int xmod = X.mod(); 04973 // loop over columns b and x of B and X 04974 for (int j=0;j<m;j++) 04975 { 04976 T* bp = B.data()+j; 04977 T* yp = y.data(); 04978 // solve L y = b (in variable x if y=0): 04979 // for i=0..n-1 04980 // y[i] = (b[i] - sum_{k<i} L[i][k] y[k])/L[i][i] 04981 for (i=0;i<n;i++,bp+=bmod) 04982 { 04983 const T* Li = L[i]; 04984 for (sum = *bp,k=i-1;k>=0;k--) sum -= Li[k] * yp[k]; 04985 yp[i] = sum / Li[i]; 04986 } 04987 // solve L' x = y 04988 // for i=n-1..0 04989 // x[i] = (y[i] - sum_{k>i} L[k][i] x[k])/L[i][i] 04990 for (i=n-1;i>=0;i--) 04991 { 04992 sum=yp[i]; 04993 if (i+1<n) 04994 { 04995 T* xp = &X(i+1,j); 04996 for (k=i+1;k<n;k++,xp+=xmod) sum -= L[k][i] * *xp; 04997 } 04998 X(i,j) = sum / L[i][i]; 04999 } 05000 } 05001 } 05002 05003 /* 05004 Back-propagate through the CholeskySolve(L,b,x,y) operation 05005 (the optional argument y of this call must have been provided). 05006 05007 dC_dL and dC_db must have been cleared beforehand. 05008 dC_dx will be modified (because of the dependencies between the x's. 05009 05010 (1) back-prop through step L' x = y: 05011 for i=0..n-1 05012 dC_dy[i] = dC_dx[i] / L[i][i] 05013 dC_dL[i][i] -= dC_dx[i] x[i] / L[i][i] 05014 for k=i+1..n 05015 dC_dx[k] -= dC_dx[i] L[k][i] / L[i][i] 05016 dC_dL[k][i] -= dC_dx[i] x[k] / L[i][i] 05017 05018 (2) back-prop through step L y = b: 05019 for i=n-1..0 05020 dC_db[i] = dC_dy[i] / L[i][i] 05021 dC_dL[i][i] -= dC_dy[i] y[i] / L[i][i] 05022 for k=0..i-1 05023 dC_dy[k] -= dC_dy[i] L[i][k] / L[i][i] 05024 dC_dL[i][k] -= dC_dy[i] * y[k] / L[i][i] 05025 */ 05026 template<class T> 05027 void bpropCholeskySolve(const TMat<T>& L, const TVec<T>& x, const TVec<T>& y, 05028 TMat<T>& dC_dL, TVec<T>& dC_db, TVec<T>& dC_dx) 05029 { 05030 int n = L.length(); 05031 int i,k; 05032 TVec<T> dC_dy(n); 05033 const T *xp = x.data(); 05034 const T *yp = y.data(); 05035 T *dC_dbp = dC_db.data(); 05036 T *dC_dxp = dC_dx.data(); 05037 T* dC_dyp = dC_dy.data(); 05038 05039 // (1) back-prop through step L' x = y: 05040 for (i=0;i<n;i++) 05041 { 05042 const T* Li = L[i]; 05043 T invLii = 1.0 / Li[i]; 05044 dC_dyp[i] = dC_dxp[i] * invLii; 05045 T dC_dxi = dC_dxp[i]; 05046 dC_dL[i][i] -= dC_dxp[i] * xp[i] * invLii; 05047 for (k=i+1;k<n;k++) 05048 { 05049 dC_dxp[k] -= dC_dxi * L[k][i] * invLii; 05050 dC_dL[k][i] -= dC_dxi * xp[k] * invLii; 05051 } 05052 } 05053 05054 // (2) back-prop through step L y = b: 05055 for (i=n-1;i>=0;i--) 05056 { 05057 const T* Li = L[i]; 05058 T* dC_dLi = dC_dL[i]; 05059 T invLii = 1.0 / Li[i]; 05060 T dC_dyi = dC_dyp[i]; 05061 T dC_dyi_over_Lii = dC_dyi * invLii; 05062 dC_dbp[i] += dC_dyi_over_Lii; 05063 dC_dLi[i] -= dC_dyi_over_Lii * yp[i]; 05064 for (k=0;k<i;k++) 05065 { 05066 dC_dyp[k] -= dC_dyi_over_Lii * Li[k]; 05067 dC_dLi[k] -= dC_dyi_over_Lii * yp[k]; 05068 } 05069 }; 05070 } 05071 05072 /* Use Cholesky decomposition to invert symmetric 05073 positive definite matrix A. 05074 Also returns the log of the determinant of A 05075 05076 We have L L' = A, and we want to solve L L' Ainv = I. 05077 05078 1) solve L Linv = I, i.e., invert L 05079 05080 for j=0..n-1 05081 Linv[j][j] = 1 / L[j][j] 05082 for i=j+1..n-1 05083 Linv[i][j] = - sum_{j<=k<i} L[i][k] Linv[k][j] / L[i][i] 05084 and 0 elsewhere (Linv is lower diagonal) 05085 05086 2) solve L' Ainv = Linv 05087 05088 for j=0..n-1 05089 for i=n-1..0 05090 Ainv[i][j] = (Linv[i][j] - sum_{k>i} L[k][i] Ainv[k][j])/L[i][i] 05091 05092 */ 05093 template<class T> 05094 real choleskyInvert(const TMat<T>& A, TMat<T>& Ainv) 05095 { 05096 int n= A.length(); 05097 TMat<T> L(n,n); 05098 Ainv.resize(n,n); 05099 05100 choleskyDecomposition(A,L); 05101 // now L L' = A 05102 05103 real logdet = log(fabs(L(0,0))); 05104 for(int i=1; i<n; i++) 05105 logdet += log(fabs(L(i,i))); 05106 logdet *= 2; 05107 05108 // Compute Linv and put its transpose above L's diagonal. 05109 // and put Linv[i][i] = 1 / L[i][i] in L's diagonal. 05110 int i,j; 05111 T *Lii = L.data(); 05112 for (i=0;i<n;i++,Lii+=1+n) 05113 *Lii = 1.0 / *Lii; 05114 05115 for (j=0;j<n;j++) 05116 { 05117 T *Linv_xj = L[j]; // Linv' in L's upper triangle 05118 for (i=j+1;i<n;i++) 05119 { 05120 T sum=0.0; 05121 T* Li = L[i]; 05122 int k; 05123 for (k=j;k<i;k++) sum -= Li[k] * Linv_xj[k]; 05124 Linv_xj[i] = sum * Li[i]; // * not / because inverse already done above 05125 } 05126 } 05127 // recall: now Linv above and on diagonal of L, L below it, 05128 05129 // compute A's inverse 05130 for (j=0;j<n;j++) 05131 { 05132 T* Linv_xj = L[j]; 05133 for (i=n-1;i>=j;i--) 05134 { 05135 T sum = Linv_xj[i]; // this is Linv[i][j] 05136 int k; 05137 for (k=i+1;k<n;k++) 05138 sum -= L[k][i] * Ainv[k][j]; 05139 Ainv[i][j] = sum * L[i][i]; 05140 } 05141 for (i=j-1;i>=0;i--) // symmetric part 05142 Ainv[i][j] = Ainv[j][i]; 05143 }; 05144 05145 return logdet; 05146 } 05147 05148 /* Solve a linear system of equations A x = b, when A is 05149 symmetric positive definite. Return x. */ 05150 template<class T> 05151 TVec<T> choleskySolve(const TMat<T>& A, const TVec<T>& b) 05152 { 05153 int n = A.length(); 05154 TMat<T> L(n,n); 05155 TVec<T> x(n); 05156 choleskyDecomposition(A,L); 05157 choleskySolve(L,b,x); 05158 return x; 05159 } 05160 05161 /* return inverse of positive definite matrix A 05162 using Cholesky decomposition. No side-effect on A. */ 05163 template<class T> 05164 TMat<T> choleskyInvert(const TMat<T>& A) 05165 { 05166 int n=A.length(); 05167 TMat<T> Ainv(n,n); 05168 choleskyInvert(A,Ainv); 05169 return Ainv; 05170 } 05171 05172 template<class T> 05173 void LU_decomposition(TMat<T>& A, TVec<T>& Trow, int& detsign, TVec<T>* p=0) 05174 { 05175 int n=A.length(); 05176 if (n!=A.width()) 05177 PLERROR("LU_decomposition: matrix A(%d,%d) should be square", n,A.width()); 05178 TVec<T>* pivot = (p==0)?new TVec<T>(n):p; 05179 T* pv = pivot->data(); 05180 detsign = 1; 05181 for (int i=0;i<n;i++) 05182 { 05183 T max_abs = maxabs(A(i)); 05184 if (max_abs==0) 05185 PLERROR("LU_decomposition: %d-th row has only zeros",i); 05186 pv[i] = 1.0 / max_abs; 05187 } 05188 int mod = A.mod(); 05189 for (int j=0;j<n;j++) 05190 { 05191 for (int i=0;i<j;i++) 05192 { 05193 T* Ai = A[i]; 05194 T* Akj = A.data()+j; 05195 T Uij = Ai[j]; 05196 for (int k=0;k<i;k++,Akj+=mod) 05197 Uij -= Ai[k] * *Akj; 05198 Ai[j] = Uij; 05199 } 05200 T max_abs = 0; 05201 int maxi = 0; 05202 for (int i=j;i<n;i++) 05203 { 05204 T* Ai = A[i]; 05205 T* Akj = A.data()+j; 05206 T Lij = Ai[j]; 05207 for (int k=0;k<j;k++,Akj+=mod) 05208 Lij -= Ai[k] * *Akj; 05209 Ai[j] = Lij; 05210 T piv = fabs(Lij) * pv[i]; 05211 if (piv >= max_abs) 05212 { 05213 maxi = i; 05214 max_abs = piv; 05215 } 05216 } 05217 if (j!=maxi) 05218 // swap row j and row maxi 05219 { 05220 A.swapRows(j,maxi); 05221 pv[maxi]=pv[j]; 05222 detsign = -detsign; 05223 } 05224 Trow[j] = maxi; 05225 T& Ajj = A(j,j); 05226 if (Ajj==0) Ajj=1e-20; // some regularization of singular matrices 05227 if (j<n-1) 05228 { 05229 T denom = 1.0/Ajj; 05230 T* Aij = &A(j+1,j); 05231 for (int i=j+1;i<n;i++, Aij+=mod) 05232 *Aij *= denom; 05233 } 05234 } 05235 if (p!=0) delete pivot; 05236 } 05237 05238 // return the determinant of A, using LU_decomposition 05239 template<class T> 05240 T det(const TMat<T>& A) 05241 { 05242 int n = A.length(); 05243 if (n!=A.width()) 05244 PLERROR("det(const TMat<T>& A): A(%d,%d) is not square!",n,A.width()); 05245 for (int i=0;i<n;i++) 05246 { 05247 T max_abs = maxabs(A(i)); 05248 if (max_abs==0) 05249 return 0.0; 05250 } 05251 TMat<T> LU(A.copy()); 05252 TVec<T> Trow(n); 05253 TVec<T> p(n); 05254 int detsign; 05255 LU_decomposition(LU, Trow, detsign); 05256 return det(LU,detsign); 05257 } 05258 05259 // return the determinant of A, whose LU decomposition is given 05260 // (detsign is as set by LU_decomposition) 05261 template<class T> 05262 T det(const TMat<T>& LU, int detsign) 05263 { 05264 T determinant = detsign; 05265 int mod = LU.mod(); 05266 int n = LU.width(); 05267 if (n!=LU.width()) 05268 PLERROR("det(const TMat<T>& LU, int detsign): LU(%d,%d) is not square!",n,LU.width()); 05269 T* LUii = LU.data(); 05270 for (int i=0;i<n;i++, LUii+=1+mod) 05271 determinant *= *LUii; 05272 return determinant; 05273 } 05274 05275 // dest[i,j] = 1 if src[i,j]==v, 0 otherwise 05276 template<class T> 05277 void equals(const TMat<T>& src, T v, TMat<T>& dest) 05278 { 05279 int l=src.length(); 05280 int w=src.width(); 05281 #ifdef BOUNDCHECK 05282 if (l!=dest.length() || w!=dest.width()) 05283 PLERROR("equals(TMat<T>(%d,%d),T,TMat<T>(%d,%d)) args of unequal dimensions", 05284 src.length(),src.width(),dest.length(),dest.width()); 05285 #endif 05286 for (int i=0;i<l;i++) 05287 { 05288 const T* s=src[i]; 05289 T* d=dest[i]; 05290 for (int j=0;j<w;j++) 05291 if (s[j]==v) d[j]=1.0; else d[j]=0.0; 05292 } 05293 } 05294 05295 // dest[i,j] = src[i,j] 05296 template<class T> 05297 void transpose(const TMat<T> src, TMat<T> dest) 05298 { 05299 int l=src.length(); 05300 int w=src.width(); 05301 #ifdef BOUNDCHECK 05302 if (w!=dest.length() || l!=dest.width()) 05303 PLERROR("tranpose(TMat<T>(%d,%d),T,TMat<T>(%d,%d)) args of unequal dimensions", 05304 src.length(),src.width(),dest.length(),dest.width()); 05305 #endif 05306 int dmod=dest.mod(); 05307 for (int i=0;i<l;i++) 05308 { 05309 const T* si=src[i]; 05310 T* dji= &dest(0,i); 05311 for (int j=0;j<w;j++,dji+=dmod) 05312 *dji = si[j]; 05313 } 05314 } 05315 05316 // res[i,j] = src[i,j] 05317 template<class T> 05318 TMat<T> transpose(const TMat<T>& src) 05319 { 05320 TMat<T> res(src.width(),src.length()); 05321 transpose(src,res); 05322 return res; 05323 } 05324 05325 // Apply a vector operation to each row of matrices, result in rows of a matrix 05326 template<class T> 05327 void apply(T (*func)(const TVec<T>&), const TMat<T>& m, TMat<T>& dest) 05328 { 05329 if (dest.length()!=m.length()) 05330 PLERROR("apply: m.length_=%d, dest.length_=%d", 05331 m.length(),dest.length()); 05332 for (int i=0;i<m.length();i++) 05333 dest(i,0)=func(m(i)); 05334 } 05335 05336 template<class T> 05337 void apply(T (*func)(const TVec<T>&,const TVec<T>&), const TMat<T>& m1, const TMat<T>& m2, 05338 TMat<T>& dest) 05339 { 05340 if (dest.length()!=m1.length() || m1.length()!=m2.length()) 05341 PLERROR("apply: m1.length_=%d, m2.length_=%d, dest.length_=%d", 05342 m1.length(),m2.length(),dest.length()); 05343 for (int i=0;i<m1.length();i++) 05344 dest(i,0)=func(m1(i),m2(i)); 05345 } 05346 05347 // Perform a traditional linear regression (but with weight decay), 05348 // without bias term. i.e. find weights such that: 05349 // 05350 // norm(weights*inputs - outputs) + weight_decay*norm(weights) 05351 // 05352 // is minimized, 05353 // 05354 // This is achieved by solving the following linear system: 05355 // 05356 // (X' X + weight_decay I) * weights = X' outputs 05357 05358 template<class T> 05359 void linearRegressionNoBias(TMat<T> inputs, TMat<T> outputs, T weight_decay, 05360 TMat<T> weights) 05361 { 05362 int inputsize = inputs.width(); 05363 int outputsize = outputs.width(); 05364 int l = inputs.length(); 05365 if(outputs.length()!=l) 05366 PLERROR("In linearRegressionNoBias: inputs and outputs should have the same length"); 05367 if(weights.length()!=inputsize || weights.width()!=outputsize) 05368 PLERROR("In linearRegressionNoBias: weights should be a (inputsize x outputsize) matrix (%d x %d)",inputsize,outputsize); 05369 TMat<T> XtX(inputsize,inputsize); 05370 transposeProduct(XtX, inputs,inputs); 05371 TMat<T> XtY(inputsize,outputsize); 05372 transposeProduct(XtY, inputs,outputs); 05373 for(int i=0; i<inputsize; i++) 05374 XtX(i,i) += weight_decay; 05375 solveLinearSystemByCholesky(XtX,XtY,weights); 05376 } 05377 05378 05379 // Perform a traditional linear regression (but with weight decay), 05380 // i.e. find bias and weights such that 05381 // 05382 // norm(bias + weights*inputs - outputs) + weight_decay*norm(weights) 05383 // 05384 // is minimized, where theta'=(biases;weights') {biases in first row} 05385 // 05386 // This is achieved by solving the following linear system: 05387 // 05388 // (X' X + weight_decay I) * theta = X' outputs 05389 // 05390 // where X = augmented inputs, i.e. X(t) = (1,inputs(t)) 05391 // 05392 template<class T> 05393 void linearRegression(TMat<T> inputs, TMat<T> outputs, T weight_decay, 05394 TMat<T> theta_t) 05395 { 05396 int l = inputs.length(); 05397 int n_inputs = inputs.width(); 05398 int n_outputs = outputs.width(); 05399 if (outputs.length()!=l) 05400 PLERROR("linearRegression: inputs.length_=%d while outputs.length_=%d", 05401 l,outputs.length()); 05402 if (theta_t.length()!=n_inputs+1 || theta_t.width()!=n_outputs) 05403 PLERROR("linearRegression: theta_t(%d,%d) should be (n_inputs(%d)+1)xn_outputs(%d)", 05404 theta_t.length(),theta_t.width(),n_inputs,n_outputs); 05405 05406 int n=n_inputs+1; 05407 05408 TMat<T> XtX(n,n); 05409 TMat<T> XtY(n,n_outputs); 05410 // compute X' X and X'Y: 05411 // XtX(i,j) = sum_t X[t,i]*X[t,j] (with X[t,0]=1, X[t,i+1]=inputs[t,i]) 05412 // YtY(i,j) = sum_t X[t,i]*Y[t,j] 05413 // 05414 int xmod=inputs.mod(); 05415 int ymod=outputs.mod(); 05416 T *xt = inputs.data(); 05417 T *yt = outputs.data(); 05418 XtX(0,0) = l; // we know the answer ahead of time for element (0,0) 05419 for (int t=0;t<l;t++,xt+=xmod,yt+=ymod) 05420 { 05421 T* xx0 = XtX.data(); 05422 T* xy0 = XtY.data(); 05423 for (int j=0;j<n_outputs;j++) 05424 xy0[j] += yt[j]; 05425 T *xxi = xx0+n; // start the inner matrix at (1,0) 05426 T *xyi = xy0+n_outputs; // start xy at (1,0) 05427 for (int i=0;i<n_inputs;i++,xxi+=n,xyi+=n_outputs) 05428 { 05429 T xti = xt[i]; 05430 xxi[0]+=xti; 05431 T *xxip=xxi+1; 05432 for (int j=0;j<i;j++) 05433 xxip[j] += xti*xt[j]; 05434 xxip[i]+=xti*xti; 05435 for (int j=0;j<n_outputs;j++) 05436 xyi[j] += xti * yt[j]; 05437 } 05438 } 05439 // now do the symmetric part of XtX 05440 T* xx = XtX.data(); 05441 T* xxi = xx+n; 05442 for (int i=1;i<n;i++,xxi+=n) 05443 { 05444 T *xx_i=xx+i; 05445 for (int j=0;j<i;j++,xx_i+=n) 05446 *xx_i = xxi[j]; 05447 } 05448 05449 // add weight_decay on the diagonal of XX' (except for the bias) 05450 T* xxii = &XtX(1,1); 05451 for (int i=0;i<n_inputs;i++,xxii+=1+n) 05452 *xxii += weight_decay; 05453 05454 // now solve by Cholesky decomposition 05455 solveLinearSystemByCholesky(XtX,XtY,theta_t); 05456 } 05457 05458 // Compute a linear fitting of 2 dimensional data resulting 05459 // in parameters m et b for y = mx + b 05460 // 1 1 05461 // Cost function used: C = - Sum[t] { (m * x_t + b - y_t)^2 } + - weight_decay * m^2 05462 // 2 2 05463 05464 template<class T> 05465 void linearRegression(TVec<T> inputs, TVec<T> outputs, T weight_decay, TVec<T> theta_t) 05466 { 05467 int npts = inputs.length(); 05468 05469 if (outputs.length()!=npts) 05470 PLERROR("linearRegression: inputs.length_=%d while outputs.length_=%d", 05471 inputs.length(),outputs.length()); 05472 if (theta_t.length()!=2) 05473 PLERROR("linearRegression: theta_t(%d) should be 2", theta_t.length()); 05474 05475 T sum_x = 0, sum_y = 0, sum_xy = 0, sum_x2 = 0, sum2_x = 0, sum2_y = 0; 05476 05477 for (int i = 0; i < npts; ++i) { 05478 sum_x += inputs[i]; 05479 sum_y += outputs[i]; 05480 sum_xy += inputs[i] * outputs[i]; 05481 sum_x2 += inputs[i] * inputs[i]; 05482 } 05483 sum2_x = sum_x * sum_x; 05484 sum2_y = sum_y * sum_y; 05485 05486 // m 05487 theta_t[1] = (sum_xy - (sum_x * sum_y) / npts) / (sum_x2 + weight_decay - sum2_x / npts); 05488 // b 05489 theta_t[0] = (sum_y - theta_t[1] * sum_x) / npts; 05490 } 05491 05492 05493 template<class T> 05494 TMat<T> smooth(TMat<T> data, int windowsize) 05495 { 05496 TVec<T> sumvec(data.width()); 05497 TMat<T> result(data.length(), data.width()); 05498 int currentwindowsize = windowsize/2; 05499 for(int k=0; k<currentwindowsize; k++) 05500 sumvec += data(k); 05501 result(0) << sumvec; 05502 //result(0) /= (T)currentwindowsize; 05503 TVec<T> res0 = result(0); 05504 res0 /= (T)currentwindowsize; 05505 result(0) << res0; 05506 05507 for(int i=0; i<data.length(); i++) 05508 { 05509 int lowi = i-(windowsize-1)/2; // lowest index of window rows (inclusive) 05510 int highi = i+windowsize/2; // highest index of window rows (inclusive) 05511 if(lowi-1>=0) // remove row lowi-1 if it exists 05512 { 05513 sumvec -= data(lowi-1); 05514 currentwindowsize--; 05515 } 05516 if(highi<data.length()) // add row highi if it exists 05517 { 05518 sumvec += data(highi); 05519 currentwindowsize++; 05520 } 05521 result(i) << sumvec; 05522 //result(i) /= (T)currentwindowsize; 05523 TVec<T> resi = result(i); 05524 resi /= (T)currentwindowsize; 05525 result(i) << resi; 05526 } 05527 05528 05529 return result; 05530 } 05531 05532 05533 template<class T> 05534 TMat<T> square(const TMat<T>& m) 05535 { 05536 TMat<T> res(m.length(), m.width()); 05537 for(int i=0; i<m.length(); i++) 05538 for(int j=0; j<m.width(); j++) 05539 res(i,j) = square(m(i,j)); 05540 return res; 05541 } 05542 05543 template<class T> 05544 TMat<T> sqrt(const TMat<T>& m) 05545 { 05546 TMat<T> res(m.length(), m.width()); 05547 for(int i=0; i<m.length(); i++) 05548 for(int j=0; j<m.width(); j++) 05549 res(i,j) = sqrt(m(i,j)); 05550 return res; 05551 } 05552 05553 template<class T> 05554 inline void affineMatrixInitialize(TMat<T> W, bool output_on_columns=true, real scale=1.0) 05555 { 05556 int n_inputs = output_on_columns?W.width():W.length(); 05557 real delta = scale/n_inputs; 05558 fill_random_uniform(W,-delta,delta); 05559 W(0).clear(); 05560 } 05561 05562 template<class T> 05563 TMat<T> grep(TMat<T> data, int col, TVec<T> values, bool exclude=false) 05564 { 05565 TMat<T> result(data.length(),data.width()); 05566 int length=0; 05567 05568 for(int i=0; i<data.length(); i++) 05569 { 05570 bool contains = values.contains(data(i,col)); 05571 if( (!exclude && contains) || (exclude && !contains) ) 05572 result(length++) << data(i); 05573 } 05574 result.resize(length,result.width()); 05575 result.compact(); // use less memory 05576 return result; 05577 } 05578 05579 05580 template<class T> 05581 void convolve(TMat<T> m, TMat<T> mask, TMat<T> result) 05582 { 05583 if(result.length() != m.length()-mask.length()+1 || result.width() != m.width()-mask.width()+1) 05584 PLERROR("In convolve(TMat<T> m, TMat<T> mask, TMat<T> result), result does not have the appropriate dimensions"); 05585 T sum; 05586 for(int i=0; i<result.length(); i++) 05587 for(int j=0; j<result.width(); j++) 05588 { 05589 T* maskptr = mask.data(); 05590 T* mptr = m[i]+j; 05591 sum = 0.0; 05592 for(int l=0; l<mask.length(); l++, maskptr += mask.mod(), mptr += m.mod()) 05593 for(int c=0; c<mask.width(); c++) 05594 sum += maskptr[c] * mptr[c]; 05595 result(i,j) = sum; 05596 } 05597 } 05598 05599 template<class T> 05600 void subsample(TMat<T> m, int thesubsamplefactor, TMat<T> result) 05601 { 05602 T sum; 05603 int norm = thesubsamplefactor * thesubsamplefactor; 05604 for(int i=0; i<result.length(); i++) 05605 for(int j=0; j<result.width(); j++) 05606 { 05607 T* mptr = m[thesubsamplefactor*i]+thesubsamplefactor*j; 05608 sum = 0.0; 05609 for(int l=0; l<thesubsamplefactor; l++, mptr += m.mod()) 05610 for(int c=0; c<thesubsamplefactor; c++) 05611 sum += mptr[c]; 05612 result(i,j) = sum/norm; 05613 } 05614 } 05615 05616 05617 template<class T> 05618 void classification_confusion_matrix(TMat<T> outputs, TMat<T> target_classes, TMat<T> confusion_matrix) 05619 { 05620 int argmax, target; 05621 T v_max, tmp; 05622 05623 for (int i=0; i<outputs.length(); i++) { 05624 // Find argmax(outputs) 05625 v_max = outputs(i,0); 05626 argmax = 0; 05627 for (int j=1; j<outputs.width(); ++j) { 05628 tmp = outputs(i,j); 05629 if (tmp > v_max) { 05630 argmax = j; 05631 v_max = tmp; 05632 } 05633 } 05634 // Update confusion matrix 05635 target = (int) target_classes(i,0); 05636 confusion_matrix(argmax, target) ++; 05637 } 05638 } 05639 05653 template<class T> 05654 int GramSchmidtOrthogonalization(TMat<T> A, T tolerance=1e-6) 05655 { 05656 int n_basis = 0; 05657 for (int i=0;i<A.length();i++) 05658 { 05659 TVec<T> Ai=A(i); 05660 if (n_basis!=i) 05661 { 05662 TVec<T> Ab = A(n_basis); 05663 Ab << Ai; 05664 Ai=Ab; 05665 } 05666 if (i>0) 05667 projectOnOrthogonalSubspace(Ai, A.subMatRows(0,n_basis)); 05668 T normAi = norm(Ai); 05669 if (normAi>1e-6) 05670 { 05671 if (normAi!=1) 05672 Ai/=normAi; 05673 n_basis++; 05674 } 05675 // else ignore row i 05676 } 05677 return n_basis; 05678 } 05679 05681 05683 template<class T> 05684 inline TVec<T> product(const TMat<T>& m, const TVec<T>& v) 05685 { TVec<T> res(m.length()); product(res, m,v); return res; } 05686 05688 template<class T> 05689 inline TVec<T> transposeProduct(const TMat<T>& m, const TVec<T>& v) 05690 { TVec<T> res(m.width()); transposeProduct(res, m,v); return res; } 05691 05693 template<class T> 05694 inline TMat<T> product(const TMat<T>& m1, const TMat<T>& m2) 05695 { TMat<T> res(m1.length(),m2.width()); product(res, m1,m2); return res; } 05696 05698 template<class T> 05699 inline TMat<T> transposeProduct(const TMat<T>& m1, const TMat<T>& m2) 05700 { TMat<T> res(m1.width(),m2.width()); transposeProduct(res, m1,m2); return res; } 05701 05703 template<class T> 05704 inline TMat<T> productTranspose(const TMat<T>& m1, const TMat<T>& m2) 05705 { TMat<T> res(m1.length(),m2.length()); productTranspose(res, m1,m2); return res; } 05706 05708 template<class T> 05709 inline TMat<T> operator+(const TMat<T>& m, const TVec<T>& v) 05710 { TMat<T> res = m.copy(); res+=v; return res; } 05711 05713 template<class T> 05714 inline TMat<T> operator-(const TMat<T>& m, const TVec<T>& v) 05715 { TMat<T> res = m.copy(); res-=v; return res; } 05716 05718 template<class T> 05719 inline TMat<T> operator*(const TMat<T>& m, const TVec<T>& v) 05720 { TMat<T> res = m.copy(); res*=v; return res; } 05721 05723 template<class T> 05724 inline TMat<T> operator/(const TMat<T>& m, const TVec<T>& v) 05725 { TMat<T> res = m.copy(); res/=v; return res; } 05726 05728 template<class T> 05729 inline TMat<T> operator/(const TMat<T>& m1, const TMat<T>& m2) 05730 { TMat<T> res = m1.copy(); res/=m2; return res; } 05731 05732 template<class T> 05733 inline void choleskySolve(const TMat<T>& L, TVec<T> b, TVec<T> x) 05734 { TVec<T> y; choleskySolve(L,b,x,y); } 05735 05737 template<class T> 05738 inline TMat<T> grep(TMat<T> data, int col, T value, bool exclude=false) 05739 { return grep(data,col,TVec<T>(1,value),exclude); } 05740 05741 template<class T> 05742 void addIfNonMissing(const TVec<T>& source, const TVec<int>& nnonmissing, TVec<T> destination) 05743 { 05744 #ifdef BOUNDCHECK 05745 if (source.length()!=nnonmissing.length() || source.length()!=destination.length()) 05746 PLERROR("addIfNonMissing: all arguments should have the same length, got %d,%d,%d\n", 05747 source.length(),nnonmissing.length(),destination.length()); 05748 #endif 05749 T* s=source.data(); 05750 T* d=destination.data(); 05751 int* n=nnonmissing.data(); 05752 int size=source.length(); 05753 for (int i=0;i<size;i++) 05754 if (finite(s[i])) 05755 { 05756 d[i] += s[i]; 05757 n[i]++; 05758 } 05759 } 05760 05761 template<class T> 05762 void addXandX2IfNonMissing(const TVec<T>& source, const TVec<int>& nnonmissing, TVec<T> somme, TVec<T> somme2) 05763 { 05764 #ifdef BOUNDCHECK 05765 if (source.length()!=nnonmissing.length() || source.length()!=somme.length() || source.length()!=somme2.length()) 05766 PLERROR("addIfNonMissing: all arguments should have the same length, got %d,%d,%d,%d\n", 05767 source.length(),nnonmissing.length(),somme.length(),somme2.length()); 05768 #endif 05769 T* s=source.data(); 05770 T* s1=somme.data(); 05771 T* s2=somme.data(); 05772 int* n=nnonmissing.data(); 05773 int size=source.length(); 05774 for (int i=0;i<size;i++) 05775 if (finite(s[i])) 05776 { 05777 s1[i] += s[i]; 05778 s2[i] += s[i]*s[i]; 05779 n[i]++; 05780 } 05781 } 05782 05783 } // end of namespace PLearn 05784 05785 05786 // Norman: replaced the code below with this wrapper 05787 SET_HASH_FUNCTION(PLearn::TVec<T>, T, v, sumsquare(v)) 05788 SET_HASH_WITH_FUNCTION(PLearn::Vec, v, sumsquare(v)) 05789 05790 //#if __GNUC__==3 && __GNUC_MINOR__>0 05791 //namespace __gnu_cxx { 05792 //#else 05793 //namespace std { 05794 //#endif 05795 // 05796 //template<class T> 05797 //struct hash<PLearn::TVec<T> > 05798 //{ 05799 // size_t operator()(PLearn::TVec<T> v) const { return hash<T>()(sumsquare(v));} 05800 //}; 05801 05802 //} // end of namespace std 05803 05804 05805 #endif // TMat_maths_impl_H