PLearn 0.1
|
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 // Copyright (C) 2005 Yoshua Bengio, Mantas Lukosevicius 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 10125 2009-04-16 19:22:03Z nouiz $ 00041 * AUTHORS: Pascal Vincent & Yoshua Bengio & Rejean Ducharme & Mantas Lukosevicius 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 #include <limits> 00052 #include <plearn/sys/Profiler.h> 00053 00054 namespace PLearn { 00055 using namespace std; 00056 00057 template <class T> 00058 TVec<T> sign(const TVec<T>& vec) 00059 { 00060 int len = vec.length(); 00061 00062 TVec<T> sign_( len ); 00063 if (len > 0) { 00064 T* v = vec.data(); 00065 T* s = sign_.data(); 00066 00067 while(--len >= 0) 00068 { 00069 *s = sign( *v ); 00070 v++; s++; 00071 } 00072 } 00073 return sign_; 00074 } 00075 00076 template <class T> 00077 void compute_sign(const TVec<T>& vec, const TVec<T>& dest) 00078 { 00079 int len = vec.length(); 00080 if (len > 0) { 00081 T* v = vec.data(); 00082 T* s = dest.data(); 00083 while(--len >= 0) 00084 { 00085 *s = sign( *v ); 00086 v++; s++; 00087 } 00088 } 00089 } 00090 00094 template <class T> 00095 bool sortedVectorsIntersect(const TVec<T>& v1, const TVec<T>& v2) 00096 { 00097 int i1=0,i2=0; 00098 do 00099 { 00100 T v1i = v1[i1]; 00101 T v2i = v2[i2]; 00102 if (v1i==v2i) return true; 00103 if (v1i<v2i) i1++; 00104 else i2++; 00105 } 00106 while (i1<v1.size() && i2<v2.size()); 00107 return false; 00108 } 00109 00110 // target is an integer between 0 and N-1 00111 // output is a vector of N discriminant functions 00112 // (each of which tries to separate class i from the others) 00113 template <class T> 00114 real one_against_all_hinge_loss(const TVec<T>& output, 00115 const int target) 00116 { 00117 int N = output.length(); 00118 T total_hinge_loss = 0; 00119 if (N > 0) { 00120 T* o = output.data(); 00121 while(--N >= 0) 00122 { 00123 if (N==target) 00124 total_hinge_loss += hinge_loss(*o,1); 00125 else 00126 total_hinge_loss += hinge_loss(*o,-1); 00127 o++; 00128 } 00129 } 00130 return total_hinge_loss; 00131 } 00132 00133 // target is an integer between 0 and N-1 00134 // output is a vector of N discriminant functions 00135 // (each of which tries to separate class i from the others) 00136 // compute derivative of hinge loss wrt each output, in d_output 00137 template <class T> 00138 void one_against_all_hinge_loss_bprop(const TVec<T>& output, 00139 const int target, 00140 TVec<T> d_output) 00141 { 00142 int N = output.length(); 00143 d_output.resize(N); 00144 if (N > 0) { 00145 T* o = output.data(); 00146 T* d_o = d_output.data(); 00147 //MNT old buggy code (opposite numbering of outputs): 00148 /*while(--N >= 0) 00149 { 00150 if (N==target) 00151 *d_o = d_hinge_loss(*o,1); 00152 else 00153 *d_o = d_hinge_loss(*o,-1); 00154 o++; d_o++; 00155 } 00156 */ 00157 for( int i = 0; i < N; i++ ) { 00158 if ( i == target ) 00159 *d_o = d_hinge_loss( *o, 1 ); 00160 else 00161 *d_o = d_hinge_loss( *o, -1 ); 00162 o++; 00163 d_o++; 00164 } 00165 } 00166 } 00167 00169 template <class T> 00170 void softmax(const TVec<T>& x, const TVec<T>& y) 00171 { 00172 int n = x.length(); 00173 if (n>0) 00174 { 00175 T* yp = y.data(); 00176 T* xp = x.data(); 00177 T maxx = max(x); 00178 real s = 0; 00179 for (int i=0;i<n;i++) 00180 s += (yp[i] = safeexp(xp[i]-maxx)); 00181 if (s == 0) PLERROR("trying to divide by 0 in softmax"); 00182 s = 1.0 / s; 00183 for (int i=0;i<n;i++) 00184 yp[i] *= s; 00185 } 00186 } 00187 00189 template <class T> 00190 void softmaxMinus(const TVec<T>& x, const TVec<T>& y) 00191 { 00192 int n = x.length(); 00193 if (n>0) 00194 { 00195 T* yp = y.data(); 00196 T* xp = x.data(); 00197 T minx = min(x); 00198 real s = 0; 00199 for (int i=0;i<n;i++) 00200 s += (yp[i] = safeexp(-xp[i]+minx)); 00201 if (s == 0) PLERROR("trying to divide by 0 in softmax"); 00202 s = 1.0 / s; 00203 for (int i=0;i<n;i++) 00204 yp[i] *= s; 00205 } 00206 } 00207 00208 // returns y = log(sofmax(x)) 00209 template <class T> 00210 void log_softmax(const TVec<T> &x, TVec<T> &y) 00211 { 00212 if (x.length() > 0) { 00213 y << x; 00214 y -= max(x); 00215 y -= logadd(y); 00216 } 00217 } 00218 00220 template <class T> 00221 void exp(const TVec<T>& x, TVec<T>& y) 00222 { 00223 y.resize(x.length()); 00224 int n = x.length(); 00225 if (!n) 00226 return; 00227 T* xp = x.data(); 00228 T* yp = y.data(); 00229 while(n--) 00230 *yp++ = exp(*xp++); 00231 } 00232 00234 template<class T> 00235 T sumsquare(const TVec<T>& x) 00236 { 00237 if (x.length() == 0) 00238 return T(0); 00239 T* v = x.data(); 00240 T res = square(v[0]); 00241 int l = x.length(); 00242 for(int i=1; i<l; i++) 00243 res += square(v[i]); 00244 return res; 00245 } 00246 00248 template<class T> 00249 T sumabs(const TVec<T>& x) 00250 { 00251 if (x.length() == 0) 00252 return T(0); 00253 T* v = x.data(); 00254 T res = (T)(fabs((real)v[0])); 00255 int l = x.length(); 00256 for(int i=1; i<l; i++) 00257 res += (T)(fabs((real)v[i])); 00258 return res; 00259 } 00260 00262 template<class T> 00263 void squareElements(const TVec<T>& x) 00264 { 00265 if (x.length() == 0) 00266 return; 00267 T* ptr = x.data(); 00268 int l = x.length(); 00269 while(l--) 00270 { 00271 *ptr *= *ptr; 00272 ++ptr; 00273 } 00274 } 00275 00277 template<class T> 00278 void squareElements(const TMat<T>& m) 00279 { 00280 if (m.size()==0) 00281 return; 00282 if(m.isCompact()) { 00283 typename TMat<T>::compact_iterator it = m.compact_begin(); 00284 typename TMat<T>::compact_iterator itend = m.compact_end(); 00285 for(; it != itend; ++it) 00286 *it = square(*it); 00287 } else { 00288 typename TMat<T>::iterator it = m.begin(); 00289 typename TMat<T>::iterator itend = m.end(); 00290 for(; it != itend; ++it) 00291 *it = square(*it); 00292 } 00293 } 00294 00296 template<class T> 00297 T sumsquare(const TMat<T>& m) 00298 { 00299 if (m.size()==0) 00300 return T(0); 00301 if(m.isCompact()) 00302 { 00303 typename TMat<T>::compact_iterator it = m.compact_begin(); 00304 typename TMat<T>::compact_iterator itend = m.compact_end(); 00305 T res = square(*it); 00306 ++it; 00307 for(; it!=itend; ++it) 00308 res += square(*it); 00309 return res; 00310 } 00311 else 00312 { 00313 typename TMat<T>::iterator it = m.begin(); 00314 typename TMat<T>::iterator itend = m.end(); 00315 T res = square(*it); 00316 ++it; 00317 for(; it!=itend; ++it) 00318 res += square(*it); 00319 return res; 00320 } 00321 } 00322 00323 00325 template<class T> 00326 T sumabs(const TMat<T>& m) 00327 { 00328 if (m.size()==0) 00329 return T(0); 00330 if(m.isCompact()) 00331 { 00332 typename TMat<T>::compact_iterator it = m.compact_begin(); 00333 typename TMat<T>::compact_iterator itend = m.compact_end(); 00334 T res = fabs(*it); 00335 ++it; 00336 for(; it!=itend; ++it) 00337 res += fabs(*it); 00338 return res; 00339 } 00340 else 00341 { 00342 typename TMat<T>::iterator it = m.begin(); 00343 typename TMat<T>::iterator itend = m.end(); 00344 T res = fabs(*it); 00345 ++it; 00346 for(; it!=itend; ++it) 00347 res += fabs(*it); 00348 return res; 00349 } 00350 } 00351 00352 // res[i,j] = scale*(mat[i,j] - avg[i] - avg[j] + mean(avg)) 00353 template<class T> 00354 void doubleCentering(const TMat<T>& mat, TVec<T>& avg, TMat<T>& res, T scale=T(1)) 00355 { 00356 T moy = mean(avg); 00357 int n=avg.length(); 00358 if (!n) 00359 return; 00360 T* a = avg.data(); 00361 if (scale==T(1)) 00362 for (int i=0;i<n;i++) 00363 { 00364 T* Mi = mat[i]; 00365 T* Ri = res[i]; 00366 T term = moy-a[i]; 00367 for (int j=0;j<n;j++) 00368 Ri[j] = Mi[j] - a[j] + term; 00369 } 00370 else 00371 for (int i=0;i<n;i++) 00372 { 00373 T* Mi = mat[i]; 00374 T* Ri = res[i]; 00375 T term = moy-a[i]; 00376 for (int j=0;j<n;j++) 00377 Ri[j] = scale*(Mi[j] - a[j] + term); 00378 } 00379 } 00380 00381 00383 template <class T> 00384 inline void multiply(const TVec<T>& source1, T source2, TVec<T>& destination) 00385 { 00386 int n=source1.length(); 00387 if (n!=destination.length()) 00388 destination.resize(n); 00389 if (!n) 00390 return; 00391 T* s1=source1.data(); 00392 T* d=destination.data(); 00393 for (int i=0;i<n;i++) 00394 d[i] = s1[i]*source2; 00395 } 00396 00397 00398 //------- These were previously in Vec_maths 00399 00402 template<class T> 00403 T sum(const TVec<T>& vec, bool ignore_missing) 00404 { 00405 double res = 0.0; 00406 if (vec.size() == 0) 00407 return res; 00408 T* v = vec.data(); 00409 for(int i=0; i<vec.length(); i++) 00410 { 00411 if (!is_missing(v[i])) res += v[i]; 00412 else if (!ignore_missing) return MISSING_VALUE; 00413 } 00414 return T(res); 00415 } 00416 00419 template<class T> 00420 T sum(const TVec<T>& vec) 00421 { 00422 T res = T(0); 00423 if (vec.size() == 0) 00424 return res; 00425 T* v = vec.data(); 00426 for(int i=0; i<vec.length(); i++) 00427 res += v[i]; 00428 return res; 00429 } 00430 00434 template<class T> 00435 T sum_of_log(const TVec<T>& vec) 00436 { 00437 double res = 0.0; 00438 if (vec.size() == 0) 00439 return res; 00440 T* v = vec.data(); 00441 for(int i=0; i<vec.length(); i++) 00442 res += pl_log(v[i]); 00443 return T(res); 00444 } 00445 00446 template<class T> 00447 T product(const TVec<T>& vec) 00448 { 00449 T res(static_cast<T>(1.0)); 00450 if (vec.size() == 0) 00451 return res; 00452 T* v = vec.data(); 00453 for(int i=0; i<vec.length(); i++) 00454 res *= v[i]; 00455 return res; 00456 } 00457 00458 /* 00459 template<class T> 00460 T mean(const TVec<T>& vec) 00461 { 00462 #ifdef BOUNDCHECK 00463 if(vec.length()==0) 00464 PLERROR("IN T mean(const TVec<T>& vec) vec has zero length"); 00465 #endif 00466 double res = 0.0; 00467 T* v = vec.data(); 00468 for(int i=0; i<vec.length(); i++) 00469 res += v[i]; 00470 return T(res/vec.length()); 00471 } 00472 */ 00473 00478 template<class T> 00479 T mean(const TVec<T>& vec, bool ignore_missing=false) 00480 { 00481 #ifdef BOUNDCHECK 00482 if(vec.length()==0) 00483 PLERROR("IN T mean(const TVec<T>& vec) vec has zero length"); 00484 #endif 00485 if (vec.size() == 0) 00486 return MISSING_VALUE; 00487 double res = 0.0; 00488 int n = 0; 00489 T* v = vec.data(); 00490 for(int i=0; i<vec.length(); i++) 00491 { 00492 if (!is_missing(v[i])) 00493 { 00494 res += v[i]; 00495 n++; 00496 } 00497 else if (!ignore_missing) 00498 return MISSING_VALUE; 00499 } 00500 00501 if (n == 0) 00502 return MISSING_VALUE; 00503 return T(res/double(n)); 00504 } 00505 00506 template<class T> 00507 T harmonic_mean(const TVec<T>& vec, bool ignore_missing=false) 00508 { 00509 #ifdef BOUNDCHECK 00510 if(vec.length()==0) 00511 PLERROR("IN T mean(const TVec<T>& vec) vec has zero length"); 00512 #endif 00513 if (vec.size() == 0) 00514 return MISSING_VALUE; 00515 double res = 0.0; 00516 int n = 0; 00517 T* v = vec.data(); 00518 for(int i=0; i<vec.length(); i++) 00519 { 00520 if (!is_missing(v[i])) 00521 { 00522 res += 1.0/v[i]; 00523 n++; 00524 } 00525 else if (!ignore_missing) 00526 return MISSING_VALUE; 00527 } 00528 00529 if (n == 0) 00530 return MISSING_VALUE; 00531 return T(double(n)/res); 00532 } 00533 00534 // This one won't really work if you have missing values in the vector 00535 // template<class T> 00536 // T avgdev(const TVec<T>& vec, T meanval) 00537 // { 00538 // #ifdef BOUNDCHECK 00539 // if(vec.length()==0) 00540 // PLERROR("IN T avgdev(const TVec<T>& vec, T meanval) vec has zero length"); 00541 // #endif 00542 // double res = 0.0; 00543 // T* v = vec.data(); 00544 // for(int i=0; i<vec.length(); i++) 00545 // res += fabs(v[i]-meanval); 00546 // return res/vec.length(); 00547 // } 00548 00549 // Does avgdev with/without missing values. ignore_missing=true ignores the missing values 00550 // and computes the avgdev without'em 00551 template<class T> 00552 T avgdev(const TVec<T>& vec, T meanval, bool ignore_missing = false) 00553 { 00554 #ifdef BOUNDCHECK 00555 if(vec.length()==0) 00556 PLERROR("IN T avgdev(const TVec<T>& vec, T meanval) vec has zero length"); 00557 #endif 00558 double res = 0.0; 00559 int n = 0; 00560 if (vec.size() == 0) 00561 return MISSING_VALUE; 00562 T* v = vec.data(); 00563 for(int i=0; i<vec.length(); i++) 00564 if (!is_missing(v[i])) 00565 { 00566 res += fabs(v[i]-meanval); 00567 n++; 00568 } 00569 else if (!ignore_missing) 00570 return MISSING_VALUE; 00571 if (n == 0) 00572 return MISSING_VALUE; 00573 else 00574 return T(res/n); 00575 } 00576 00577 template<class T> 00578 T geometric_mean(const TVec<T>& vec) 00579 { 00580 #ifdef BOUNDCHECK 00581 if(vec.length()==0) 00582 PLERROR("IN T geometric_mean(const TVec<T>& vec) vec has zero length"); 00583 #endif 00584 if (vec.size() == 0) 00585 return MISSING_VALUE; 00586 double res = 0.0; 00587 T* v = vec.data(); 00588 for(int i=0; i<vec.length(); i++) 00589 { 00590 T vi = v[i]; 00591 if (vi<=0) 00592 PLERROR("geometric_mean(TVec<T>): argument %g <=0 at position [%d]", 00593 vi,i); 00594 res += v[i]; 00595 } 00596 return T(exp(res/vec.length())); 00597 } 00598 00599 template<class T> 00600 T weighted_mean(const TVec<T>& vec, const TVec<T>& weights, bool ignore_missing=false) 00601 { 00602 #ifdef BOUNDCHECK 00603 if(vec.length()!=weights.length() || vec.length() == 0) 00604 PLERROR("IN T weighted_mean(const TVec<T>& vec, const TVec<T>& weights) vec and weights must have equal (non-zero) lengths"); 00605 #endif 00606 if (vec.size() == 0) 00607 return MISSING_VALUE; 00608 double res = 0.0; 00609 T sum_weights = 0.0; 00610 T* v = vec.data(); 00611 T* w = weights.data(); 00612 for(int i=0; i<vec.length(); i++) 00613 { 00614 if (!is_missing(v[i]) && !is_missing(w[i])) 00615 { 00616 res += v[i] * w[i]; 00617 sum_weights += w[i]; 00618 } 00619 else if (!ignore_missing) return MISSING_VALUE; 00620 } 00621 if (sum_weights == 0) 00622 PLERROR("IN T weighted_mean: sum(weights) == 0"); 00623 return T(res/sum_weights); 00624 } 00625 00626 // ignore_missing = true means that it computes the variance ignoring 00627 // the missing values 00628 template<class T> 00629 T variance(const TVec<T>& vec, T meanval, bool ignore_missing=false) 00630 { 00631 #ifdef BOUNDCHECK 00632 if(vec.length()<=1) 00633 PLERROR("IN T variance(const TVec<T>& vec, T meanval) vec length must be more than one"); 00634 #endif 00635 if (vec.size() == 0) 00636 return MISSING_VALUE; 00637 double res = 0.0; 00638 T* v = vec.data(); 00639 int n = 0; 00640 for(int i=0; i<vec.length(); i++) 00641 { 00642 if (!is_missing(v[i])) 00643 { 00644 double diff = v[i]-meanval; 00645 res += diff*diff; 00646 n++; 00647 } 00648 else if (!ignore_missing) 00649 return MISSING_VALUE; 00650 } 00651 if (n == 0) 00652 return MISSING_VALUE; 00653 else 00654 return T(res/n); 00655 } 00656 00657 template<class T> 00658 T covariance(const TVec<T>& vec1, const TVec<T>& vec2, T mean1, T mean2) 00659 { 00660 #ifdef BOUNDCHECK 00661 if(vec1.length()<=1) 00662 PLERROR("IN T covariance(const TVec<T>& vec1, const TVec<T>& vec2, T mean1, T mean2) vec1's length must be more than one"); 00663 if(vec2.length()<=1) 00664 PLERROR("IN T covariance(const TVec<T>& vec1, const TVec<T>& vec2, T mean1, T mean2) vec2's length must be more than one"); 00665 if(vec1.length() != vec2.length()) 00666 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"); 00667 #endif 00668 if (vec1.size() == 0 || vec2.size() == 0) 00669 return MISSING_VALUE; 00670 int length = vec1.length(); 00671 double res = 0.0; 00672 T* v1 = vec1.data(); 00673 T* v2 = vec2.data(); 00674 for(int i=0; i<length; i++) 00675 { 00676 double temp = (v1[i]-mean1)*(v2[i]-mean2); 00677 res += temp; 00678 } 00679 return res/(length - 1); 00680 } 00681 00682 template<class T> 00683 T weighted_variance(const TVec<T>& vec, const TVec<T>& weights, T no_weighted_mean, T weighted_mean) 00684 { 00685 #ifdef BOUNDCHECK 00686 if(vec.length()!=weights.length() || vec.length()==0) 00687 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"); 00688 #endif 00689 if (vec.size() == 0) 00690 return MISSING_VALUE; 00691 double res = 0.0; 00692 T* v = vec.data(); 00693 T* w = weights.data(); 00694 for(int i=0; i<vec.length(); i++) 00695 res += v[i] * v[i] * w[i]; 00696 T sum_weights = sum(weights, false); 00697 if (sum_weights == 0) 00698 PLERROR("IN T weighted_variance(const TVec<T>& vec, const TVec<T>& weights, T no_weighted_mean, T weighted_mean) sum(weights) == 0"); 00699 return (res/sum_weights - no_weighted_mean * (2*weighted_mean - no_weighted_mean))*vec.length()/(vec.length()-1); 00700 } 00701 00702 template<class T> 00703 TVec<T> histogram(const TVec<T>& vec, T minval, T maxval, int nbins) 00704 { 00705 TVec<T> histo(nbins); 00706 T deltaval = maxval-minval + 1e-6; 00707 for(int i=0; i<vec.length(); i++) 00708 { 00709 T val = vec[i]; 00710 int binpos = int((val-minval)/deltaval*nbins); 00711 if(binpos>=0 && binpos<nbins) 00712 histo[binpos]++; 00713 } 00714 return histo; 00715 } 00716 00717 00719 template <class T> 00720 T max(const TVec<T>& vec) 00721 { 00722 #ifdef BOUNDCHECK 00723 if(vec.length()==0) 00724 PLERROR("IN max(const NumericVec& vec) TVec has zero length()"); 00725 #endif 00726 int n = vec.length(); 00727 if (!n) 00728 return std::numeric_limits<T>::min(); 00729 T* pv = vec.data(); 00730 T maxval = *pv++; 00731 while(--n) 00732 { 00733 if(*pv>maxval) 00734 maxval = *pv; 00735 ++pv; 00736 } 00737 return maxval; 00738 } 00739 00741 template <class T> 00742 T max(const TVec<T>& vec, int& argmax) 00743 { 00744 PLASSERT(vec.length() != 0); 00745 00746 int n = vec.length(); 00747 if (n == 0) 00748 { 00749 argmax = -1; 00750 return std::numeric_limits<T>::min(); 00751 } 00752 T* pv = vec.data(); 00753 T maxval = *pv++; 00754 argmax = 0; 00755 for (int i=1; i<vec.length(); i++,pv++) 00756 if (*pv>maxval) 00757 { 00758 maxval = *pv; 00759 argmax = i; 00760 } 00761 return maxval; 00762 } 00763 00765 template<class T> 00766 T min(const TVec<T>& vec) 00767 { 00768 #ifdef BOUNDCHECK 00769 if(vec.length()==0) 00770 PLERROR("IN T min(const TVec<T>& vec) vec has zero length"); 00771 #endif 00772 if (vec.size() == 0) 00773 return std::numeric_limits<T>::max(); 00774 T* v = vec.data(); 00775 T minval = v[0]; 00776 for(int i=1; i<vec.length(); i++) 00777 if(v[i]<minval) 00778 minval = v[i]; 00779 return minval; 00780 } 00781 00783 template <class T> 00784 T min(const TVec<T>& vec, int& argmin) 00785 { 00786 PLASSERT(vec.length() != 0); 00787 00788 int n = vec.length(); 00789 if (n == 0) 00790 { 00791 argmin = -1; 00792 return std::numeric_limits<T>::max(); 00793 } 00794 T* pv = vec.data(); 00795 T minval = *pv++; 00796 argmin = 0; 00797 for (int i=1; i<vec.length(); i++,pv++) 00798 if (*pv<minval) 00799 { 00800 minval = *pv; 00801 argmin = i; 00802 } 00803 return minval; 00804 } 00805 00807 template<class T> 00808 T maxabs(const TVec<T>& vec) 00809 { 00810 #ifdef BOUNDCHECK 00811 if(vec.length()==0) 00812 PLERROR("IN T maxabs(const TVec<T>& vec) vec has zero length"); 00813 #endif 00814 if (vec.size() == 0) 00815 return std::numeric_limits<T>::min(); 00816 T* v = vec.data(); 00817 T maxval = fabs(v[0]); 00818 for(int i=1; i<vec.length(); i++) 00819 { 00820 T a=fabs(v[i]); 00821 if(a>maxval) 00822 maxval = a; 00823 } 00824 return maxval; 00825 } 00826 00828 template <class T> 00829 T maxabs(const TVec<T>& vec, int& argmax) 00830 { 00831 PLASSERT(vec.length() != 0); 00832 00833 int n = vec.length(); 00834 if (n == 0) 00835 { 00836 argmax = -1; 00837 return std::numeric_limits<T>::min(); 00838 } 00839 T* pv = vec.data(); 00840 T maxval = fabs(*pv++); 00841 argmax = 0; 00842 for (int i=1; i<vec.length(); i++,pv++) 00843 { 00844 T a = fabs(*pv); 00845 if (a>maxval) 00846 { 00847 maxval = a; 00848 argmax = i; 00849 } 00850 } 00851 return maxval; 00852 } 00853 00855 template<class T> 00856 T minabs(const TVec<T>& vec) 00857 { 00858 #ifdef BOUNDCHECK 00859 if(vec.length()==0) 00860 PLERROR("IN T minabs(const TVec<T>& vec) vec has zero length"); 00861 #endif 00862 int n = vec.length(); 00863 PLASSERT( n >= 1 ); 00864 T* v = vec.data(); 00865 T minval = fabs(v[0]); 00866 for(int i=1; i<n; i++) 00867 { 00868 T a=fabs(v[i]); 00869 if(a<minval) 00870 minval = a; 00871 } 00872 00873 return minval; 00874 } 00875 00877 template <class T> 00878 T minabs(const TVec<T>& vec, int& argmin) 00879 { 00880 PLASSERT(vec.length() != 0); 00881 00882 int n = vec.length(); 00883 if (n == 0) 00884 { 00885 argmin = -1; 00886 return std::numeric_limits<T>::max(); 00887 } 00888 T* pv = vec.data(); 00889 T minval = fabs(*pv++); 00890 argmin = 0; 00891 for (int i=1; i<n; i++,pv++) 00892 { 00893 T a = fabs(*pv); 00894 if (a<minval) 00895 { 00896 minval = a; 00897 argmin = i; 00898 } 00899 } 00900 return minval; 00901 } 00902 00903 template<class T> 00904 int argmax(const TVec<T>& vec) 00905 { 00906 #ifdef BOUNDCHECK 00907 if(vec.length()==0) 00908 PLERROR("IN int argmax(const TVec<T>& vec) vec has zero length"); 00909 #endif 00910 T* v = vec.data(); 00911 int indexmax = 0; 00912 T maxval = v[0]; 00913 for(int i=1; i<vec.length(); i++) 00914 if(v[i]>maxval) 00915 { 00916 maxval = v[i]; 00917 indexmax = i; 00918 } 00919 return indexmax; 00920 } 00921 00922 template<class T> 00923 int argmax(const TVec<T>& vec, bool ignore_missing) 00924 { 00925 #ifdef BOUNDCHECK 00926 if(vec.length()==0) 00927 PLERROR("IN int argmax(const TVec<T>& vec) vec has zero length"); 00928 #endif 00929 T* v = vec.data(); 00930 int indexmax = -1; 00931 T maxval = MISSING_VALUE; 00932 00933 for(int i=0; i<vec.length(); i++) 00934 { 00935 if( is_missing(v[i]) ) 00936 { 00937 if(ignore_missing) continue; 00938 else PLERROR("argmax(const TVec<T>& vec, bool ignore_missing) encountered a MISSING_VALUE\n" 00939 "at index %d and ignore_missing is false.", i); 00940 } 00941 00942 if( indexmax == -1 || 00943 v[i] > maxval ) 00944 { 00945 maxval = v[i]; 00946 indexmax = i; 00947 } 00948 } 00949 return indexmax; 00950 } 00951 00952 00953 template<class T> 00954 int argmin(const TVec<T>& vec) 00955 { 00956 #ifdef BOUNDCHECK 00957 if(vec.length()==0) 00958 PLERROR("IN int argmin(const TVec<T>& vec) vec has zero length"); 00959 #endif 00960 T* v = vec.data(); 00961 int indexmin = 0; 00962 T minval = v[0]; 00963 for(int i=1; i<vec.length(); i++) 00964 if(v[i]<minval) 00965 { 00966 minval = v[i]; 00967 indexmin = i; 00968 } 00969 return indexmin; 00970 } 00971 00972 template<class T> 00973 int argmin(const TVec<T>& vec, bool ignore_missing) 00974 { 00975 #ifdef BOUNDCHECK 00976 if(vec.length()==0) 00977 PLERROR("IN int argmin(const TVec<T>& vec) vec has zero length"); 00978 #endif 00979 T* v = vec.data(); 00980 int indexmin = -1; 00981 T minval = MISSING_VALUE; 00982 00983 for(int i=0; i<vec.length(); i++) 00984 { 00985 if( is_missing(v[i]) ) 00986 { 00987 if(ignore_missing) continue; 00988 else PLERROR("argmin(const TVec<T>& vec, bool ignore_missing) encountered a MISSING_VALUE\n" 00989 "at index %d and ignore_missing is false.", i); 00990 } 00991 00992 if( indexmin == -1 || 00993 v[i] < minval ) 00994 { 00995 minval = v[i]; 00996 indexmin = i; 00997 } 00998 } 00999 return indexmin; 01000 } 01001 01002 01003 01004 template<class T> 01005 T pownorm(const TVec<T>& vec, double n) 01006 { 01007 double result = 0.0; 01008 if (vec.size() == 0) 01009 return result; 01010 T* v = vec.data(); 01011 if(n==1.0) 01012 { 01013 for(int i=0; i<vec.length(); i++) 01014 { 01015 T val = v[i]; 01016 if(val>=0) 01017 result += val; 01018 else 01019 result -= val; 01020 } 01021 } 01022 else if(n==2.0) 01023 { 01024 for(int i=0; i<vec.length(); i++) 01025 { 01026 T val = v[i]; 01027 result += val*val; 01028 } 01029 } 01030 else if(n==0) 01031 { result = vec.length(); } 01032 else 01033 { 01034 for(int i=0; i<vec.length(); i++) 01035 result += mypow(fabs(v[i]),n); 01036 } 01037 return result; 01038 } 01039 01040 template<class T> 01041 inline T pownorm(const TVec<T>& vec) { return pownorm(vec,T(2.0)); } 01042 01043 template<class T> 01044 T norm(const TVec<T>& vec, double n) 01045 { 01046 if(n==T(1.0)) 01047 return pownorm(vec, T(1.0)); 01048 else if(n==T(2.0)) 01049 return sqrt(pownorm(vec,T(2.0))); 01050 else 01051 return mypow(pownorm(vec,n), T(1.0)/n); 01052 } 01053 01054 template<class T> 01055 inline T norm(const TVec<T>& vec) { return norm(vec,T(2.0)); } 01056 01057 template<class T> 01058 void normalize(const TVec<T>& vec, double n=2) 01059 { vec /= norm(vec,n); } 01060 01065 template<class T> 01066 T powdistance(const TVec<T>& vec1, const TVec<T>& vec2, double n, 01067 bool ignore_missing = false) 01068 { 01069 #ifdef BOUNDCHECK 01070 if(vec1.length() != vec2.length()) 01071 PLERROR("In weighted_powdistance: vec1, vec2 should have the same length (%d!=%d)", 01072 vec1.length(), vec2.length()); 01073 #endif 01074 int length = vec1.length(); 01075 if (length == 0) 01076 return 0.0; 01077 T result = 0; 01078 T diff = 0; 01079 T* v1 = vec1.data(); 01080 T* v2 = vec2.data(); 01081 if(fast_exact_is_equal(n, 1.0)) // L1 distance 01082 { 01083 for(int i=0; i<length; i++, v1++, v2++) 01084 if (!ignore_missing || (!is_missing(*v1) && !is_missing(*v2))) { 01085 diff = *v1 - *v2; 01086 if(diff >= 0) 01087 result += diff; 01088 else 01089 result -= diff; 01090 } 01091 } 01092 else if(fast_exact_is_equal(n, 2.0)) 01093 { 01094 for(int i=0; i<length; i++, v1++, v2++) 01095 if (!ignore_missing || (!is_missing(*v1) && !is_missing(*v2))) { 01096 diff = *v1 - *v2; 01097 result += diff*diff; 01098 } 01099 } 01100 else 01101 { 01102 for(int i=0; i<length; i++, v1++, v2++) 01103 if (!ignore_missing || (!is_missing(*v1) && !is_missing(*v2))) { 01104 diff = *v1 - *v2; 01105 if(diff<0) 01106 diff = -diff; 01107 result += mypow(diff,n); 01108 } 01109 } 01110 return result; 01111 } 01112 01113 template<class T> 01114 inline T powdistance(const TVec<T>& vec1, const TVec<T>& vec2) 01115 { return powdistance(vec1, vec2, 2.0); } 01116 01117 template<class T> 01118 T dist(const TVec<T>& vec1, const TVec<T>& vec2, double n) 01119 { 01120 if(fast_exact_is_equal(n, T(1.0))) 01121 return powdistance(vec1, vec2, T(1.0)); 01122 else if(fast_exact_is_equal(n, T(2.0))) 01123 return sqrt(powdistance(vec1, vec2, T(2.0))); 01124 else 01125 return mypow(powdistance(vec1, vec2, n), T(1.0)/n); 01126 } 01127 01128 template<class T> 01129 inline T L2distance(const TVec<T>& vec1, const TVec<T>& vec2) 01130 { return dist(vec1, vec2, 2.0); } 01131 01132 template<class T> 01133 inline T L1distance(const TVec<T>& vec1, const TVec<T>& vec2) 01134 { return dist(vec1, vec2, 1.0); } 01135 01136 01137 template<class T> 01138 T weighted_powdistance(const TVec<T>& vec1, const TVec<T>& vec2, double n, const TVec<T>& weights) 01139 { 01140 #ifdef BOUNDCHECK 01141 if(vec1.length() != weights.length() || vec2.length()!=weights.length()) 01142 PLERROR("In weighted_powdistance: vec1, vec2 and weights vector should have the same length"); 01143 #endif 01144 T result = 0.0; 01145 if (vec1.size() > 0 && vec2.size() > 0 && weights.size() > 0) { 01146 T* v1 = vec1.data(); 01147 T* v2 = vec2.data(); 01148 T* w = weights.data(); 01149 int length = vec1.length(); 01150 if(n==1.0) // L1 distance 01151 { 01152 for(int i=0; i<length; i++) 01153 { 01154 T diff = w[i]*(v1[i]-v2[i]); 01155 if(diff>=0) 01156 result += diff; 01157 else 01158 result -= diff; 01159 } 01160 } 01161 else if(n==2.0) 01162 { 01163 for(int i=0; i<length; i++) 01164 { 01165 T diff = w[i]*(v1[i]-v2[i]); 01166 result += diff*diff; 01167 } 01168 } 01169 else 01170 { 01171 for(int i=0; i<length; i++) 01172 { 01173 T diff = w[i]*(v1[i]-v2[i]); 01174 if(diff<0) 01175 diff = -diff; 01176 result += mypow(diff,n); 01177 } 01178 } 01179 } 01180 return result; 01181 } 01182 01183 template<class T> 01184 T weighted_distance(const TVec<T>& vec1, const TVec<T>& vec2, double n, const TVec<T>& weights) 01185 { 01186 if(n==1.0) 01187 return weighted_powdistance(vec1, vec2, 1.0, weights); 01188 else if(n==2.0) 01189 return sqrt(weighted_powdistance(vec1, vec2, 2.0, weights)); 01190 else 01191 return mypow(weighted_powdistance(vec1, vec2, n, weights), 1.0/n); 01192 } 01193 01194 01196 template<class T> 01197 inline void operator+=(const TVec<T>& vec1, const TVec<T>& vec2) 01198 { 01199 #ifdef BOUNDCHECK 01200 if(vec1.size() != vec2.size()) 01201 PLERROR("In operator+=, vec1 and vec2 vectors must have the same length"); 01202 #endif 01203 if (vec1.size() > 0 && vec2.size() > 0) { 01204 T* v1 = vec1.data(); 01205 T* v2 = vec2.data(); 01206 int l = vec1.length(); 01207 for(int i=0; i<l; i++) 01208 *v1++ += *v2++; 01209 } 01210 } 01211 01212 template<class T> 01213 void operator+=(const TVec<T>& vec, T scalar) 01214 { 01215 if (vec.size() > 0) { 01216 T* v = vec.data(); 01217 for(int i=0; i<vec.length(); i++) 01218 v[i] += scalar; 01219 } 01220 } 01221 01222 template<class T> 01223 TVec<T> operator-(const TVec<T>& vec) 01224 { 01225 if (vec.size() > 0) { 01226 TVec<T> opposite(vec.length()); 01227 T *v=vec.data(); 01228 T *o=opposite.data(); 01229 for (int i=0;i<vec.length();i++) 01230 o[i] = - v[i]; 01231 return opposite; 01232 } 01233 return TVec<T>(); 01234 } 01235 01236 template<class T> 01237 void operator-=(const TVec<T>& vec1, const TVec<T>& vec2) 01238 { 01239 #ifdef BOUNDCHECK 01240 if(vec1.size() != vec2.size()) 01241 PLERROR("In operator-=, vec1 and vec2 vectors must have the same length"); 01242 #endif 01243 if (vec1.size() > 0 && vec2.size() > 0) { 01244 T* v1 = vec1.data(); 01245 T* v2 = vec2.data(); 01246 for(int i=0; i<vec1.length(); i++) 01247 v1[i] -= v2[i]; 01248 } 01249 } 01250 01251 template<class T> 01252 void operator-=(const TVec<T>& vec, T scalar) 01253 { vec += -scalar; } 01254 01255 template<class T> 01256 void operator*=(const TVec<T>& vec1, const TVec<T>& vec2) 01257 { 01258 #ifdef BOUNDCHECK 01259 if(vec1.size() != vec2.size()) 01260 PLERROR("In operator*=, vec1 and vec2 vectors must have the same length"); 01261 #endif 01262 if (vec1.size() > 0 && vec2.size() > 0) { 01263 T* v1 = vec1.data(); 01264 T* v2 = vec2.data(); 01265 for(int i=0; i<vec1.length(); i++) 01266 v1[i] *= v2[i]; 01267 } 01268 } 01269 01270 template<class T> 01271 void operator*=(const TVec<T>& vec, T factor) 01272 { 01273 if (vec.size() > 0) { 01274 T* p = vec.data(); 01275 int l = vec.length(); 01276 for (int i=0;i<l;i++) 01277 *p++ *= factor; 01278 } 01279 } 01280 01281 template<class T> 01282 void operator/=(const TVec<T>& vec1, const TVec<T>& vec2) 01283 { 01284 #ifdef BOUNDCHECK 01285 if(vec1.size() != vec2.size()) 01286 PLERROR("In operator/=, vec1 and vec2 vectors must have the same length"); 01287 #endif 01288 if (vec1.size() > 0 && vec2.size() > 0) { 01289 T* v1 = vec1.data(); 01290 T* v2 = vec2.data(); 01291 int l=vec1.length(); 01292 for(int i=0; i<l; i++) 01293 v1[i] /= v2[i]; 01294 } 01295 } 01296 01297 template<class T> 01298 inline void operator/=(const TVec<T>& vec, T scalar) 01299 { vec *= T(1.0)/scalar; } 01300 01301 template<class T> 01302 inline void operator/=(const TVec<T>& vec, int scalar) 01303 { vec /= T(scalar); } 01304 01305 template<class T> 01306 void compute_log(const TVec<T>& src, const TVec<T>& dest) 01307 { 01308 #ifdef BOUNDCHECK 01309 if(src.length()!=dest.length()) 01310 PLERROR("In log, src and dest vectors must have the same length"); 01311 #endif 01312 if (src.size() > 0 && dest.size() > 0) { 01313 T* ps = src.data(); 01314 T* pd = dest.data(); 01315 int n = src.length(); 01316 for(int i=0; i<n; i++) 01317 *pd++ = pl_log(*ps++); 01318 } 01319 } 01320 01321 template<class T> 01322 inline TVec<T> log(const TVec<T>& src) 01323 { TVec<T> dest(src.length()); compute_log(src,dest); return dest; } 01324 01325 template<class T> 01326 void compute_sqrt(const TVec<T>& src, const TVec<T>& dest) 01327 { 01328 #ifdef BOUNDCHECK 01329 if(src.length()!=dest.length()) 01330 PLERROR("In sqrt, src and dest vectors must have the same length"); 01331 #endif 01332 if (src.size() > 0 && dest.size() > 0) { 01333 T* ps = src.data(); 01334 T* pd = dest.data(); 01335 int n = src.length(); 01336 for(int i=0; i<n; i++) 01337 *pd++ = sqrt(*ps++); 01338 } 01339 } 01340 01341 template<class T> 01342 inline TVec<T> sqrt(const TVec<T>& src) 01343 { TVec<T> dest(src.length()); compute_sqrt(src,dest); return dest; } 01344 01345 template<class T> 01346 void compute_safelog(const TVec<T>& src, const TVec<T>& dest) 01347 { 01348 #ifdef BOUNDCHECK 01349 if(src.length()!=dest.length()) 01350 PLERROR("In safelog, src and dest vectors must have the same length"); 01351 #endif 01352 if (src.size() > 0 && dest.size() > 0) { 01353 T* ps = src.data(); 01354 T* pd = dest.data(); 01355 int n = src.length(); 01356 for(int i=0; i<n; i++) 01357 *pd++ = safelog(*ps++); 01358 } 01359 } 01360 01361 template<class T> 01362 inline TVec<T> safelog(const TVec<T>& src) 01363 { TVec<T> dest(src.length()); compute_safelog(src,dest); return dest; } 01364 01365 template<class T> 01366 void compute_tanh(const TVec<T>& src, const TVec<T>& dest) 01367 { 01368 #ifdef BOUNDCHECK 01369 if(src.length()!=dest.length()) 01370 PLERROR("In tanh, src and dest vectors must have the same length"); 01371 #endif 01372 if (src.size() > 0 && dest.size() > 0) { 01373 T* ps = src.data(); 01374 T* pd = dest.data(); 01375 int n = src.length(); 01376 for(int i=0; i<n; i++) 01377 *pd++ = tanh(*ps++); 01378 } 01379 } 01380 01381 template<class T> 01382 void bprop_tanh(const TVec<T>& tanh_x, const TVec<T>& d_tanh_x, TVec<T>& d_x) 01383 { 01384 #ifdef BOUNDCHECK 01385 if(tanh_x.length()!=d_tanh_x.length()) 01386 PLERROR("In bprop_tanh, src and dest vectors must have the same length"); 01387 #endif 01388 if (tanh_x.size() > 0 && d_tanh_x.size() > 0 && d_x.size() > 0) { 01389 int n = tanh_x.length(); 01390 if (n != d_x.length()) d_x.resize(n); 01391 T* y = tanh_x.data(); 01392 T* dy = d_tanh_x.data(); 01393 T* dx = d_x.data(); 01394 for(int i=0; i<n; i++) 01395 { 01396 real yi = *y++; 01397 *dx++ = *dy++ * (1 - yi*yi); 01398 } 01399 } 01400 } 01401 01402 template<class T> 01403 inline TVec<T> tanh(const TVec<T>& src) 01404 { TVec<T> dest(src.length()); compute_tanh(src,dest); return dest; } 01405 01406 01407 template<class T> 01408 void compute_fasttanh(const TVec<T>& src, const TVec<T>& dest) 01409 { 01410 #ifdef BOUNDCHECK 01411 if(src.length()!=dest.length()) 01412 PLERROR("In fasttanh, src and dest vectors must have the same length"); 01413 #endif 01414 if (src.size() > 0 && dest.size() > 0) { 01415 T* ps = src.data(); 01416 T* pd = dest.data(); 01417 int n = src.length(); 01418 for(int i=0; i<n; i++) 01419 *pd++ = fasttanh(*ps++); 01420 } 01421 } 01422 01423 template<class T> 01424 inline TVec<T> fasttanh(const TVec<T>& src) 01425 { TVec<T> dest(src.length()); compute_fasttanh(src,dest); return dest; } 01426 01427 template<class T> 01428 void compute_sigmoid(const TVec<T>& src, const TVec<T>& dest) 01429 { 01430 #ifdef BOUNDCHECK 01431 if(src.length()!=dest.length()) 01432 PLERROR("In sigmoid, src and dest vectors must have the same length"); 01433 #endif 01434 if (src.size() > 0 && dest.size() > 0) { 01435 T* ps = src.data(); 01436 T* pd = dest.data(); 01437 int n = src.length(); 01438 for(int i=0; i<n; i++) 01439 *pd++ = sigmoid(*ps++); 01440 } 01441 } 01442 01443 template<class T> 01444 void log_sigmoid(const TVec<T>& src, const TVec<T>& dest) 01445 { 01446 #ifdef BOUNDCHECK 01447 if(src.length()!=dest.length()) 01448 PLERROR("In sigmoid, src and dest vectors must have the same length"); 01449 #endif 01450 if (src.size() > 0 && dest.size() > 0) { 01451 T* ps = src.data(); 01452 T* pd = dest.data(); 01453 int n = src.length(); 01454 for(int i=0; i<n; i++) 01455 *pd++ = log_sigmoid(*ps++); 01456 } 01457 } 01458 01459 template<class T> 01460 inline TVec<T> sigmoid(const TVec<T>& src) 01461 { TVec<T> dest(src.length()); compute_sigmoid(src,dest); return dest; } 01462 01463 01464 template<class T> 01465 void compute_fastsigmoid(const TVec<T>& src, const TVec<T>& dest) 01466 { 01467 #ifdef BOUNDCHECK 01468 if(src.length()!=dest.length()) 01469 PLERROR("In fastsigmoid, src and dest vectors must have the same length"); 01470 #endif 01471 if (src.size() > 0 && dest.size() > 0) { 01472 T* ps = src.data(); 01473 T* pd = dest.data(); 01474 int n = src.length(); 01475 for(int i=0; i<n; i++) 01476 *pd++ = fastsigmoid(*ps++); 01477 } 01478 } 01479 01480 template<class T> 01481 inline TVec<T> fastsigmoid(const TVec<T>& src) 01482 { TVec<T> dest(src.length()); compute_fastsigmoid(src,dest); return dest; } 01483 01484 template<class T> 01485 void compute_inverse_sigmoid(const TVec<T>& src, const TVec<T>& dest) 01486 { 01487 #ifdef BOUNDCHECK 01488 if(src.length()!=dest.length()) 01489 PLERROR("In inverse_sigmoid, src and dest vectors must have the same length"); 01490 #endif 01491 if (src.size() > 0 && dest.size() > 0) { 01492 T* ps = src.data(); 01493 T* pd = dest.data(); 01494 int n = src.length(); 01495 for(int i=0; i<n; i++) 01496 *pd++ = inverse_sigmoid(*ps++); 01497 } 01498 } 01499 01500 template<class T> 01501 inline TVec<T> inverse_sigmoid(const TVec<T>& src) 01502 { TVec<T> dest(src.length()); compute_inverse_sigmoid(src,dest); return dest; } 01503 01504 01505 template<class T> 01506 void negateElements(const TVec<T>& vec) 01507 { 01508 if (vec.size() > 0) { 01509 T* v = vec.data(); 01510 for(int i=0; i<vec.length(); i++) 01511 v[i] = -v[i]; 01512 } 01513 } 01514 01515 template<class T> 01516 void invertElements(const TVec<T>& vec) 01517 { 01518 if (vec.size() > 0) { 01519 T* v = vec.data(); 01520 for(int i=0; i<vec.length(); i++) 01521 v[i] = 1.0/v[i]; 01522 } 01523 } 01524 01525 template<class T> 01526 TVec<T> inverted(const TVec<T>& vec) 01527 { 01528 TVec<T> ret(vec.length()); 01529 if (vec.size() > 0) { 01530 T* v = vec.data(); 01531 for(int i=0; i<vec.length(); i++) 01532 ret[i] = 1.0/v[i]; 01533 } 01534 return ret; 01535 } 01536 01537 01538 template<class T> 01539 T dot(const TVec<T>& vec1, const TVec<T>& vec2) 01540 { 01541 #ifdef BOUNDCHECK 01542 if(vec1.length()!=vec2.length()) 01543 PLERROR("In T operator*(const TVec<T>& vec1, const TVec<T>& vec2) (dot product) the 2 vecs must have the same length."); 01544 #endif 01545 T res = 0; 01546 if (vec1.size() > 0 && vec2.size() > 0) { 01547 T* v1 = vec1.data(); 01548 T* v2 = vec2.data(); 01549 for(int i=0; i<vec1.length(); i++) 01550 res += v1[i]*v2[i]; 01551 } 01552 return res; 01553 } 01554 01560 template<class V, class T, class U> 01561 V dot(const TVec<T>& vec1, const TVec<U>& vec2) 01562 { 01563 #ifdef BOUNDCHECK 01564 if(vec1.length()!=vec2.length()) 01565 PLERROR("In T operator*(const TVec<T>& vec1, const TVec<T>& vec2) (dot product) the 2 vecs must have the same length."); 01566 #endif 01567 V res = 0; 01568 if (vec1.size() > 0 && vec2.size() > 0) { 01569 T* v1 = vec1.data(); 01570 U* v2 = vec2.data(); 01571 for(int i=0; i<vec1.length(); i++) 01572 res += v1[i]*v2[i]; 01573 } 01574 return res; 01575 } 01576 01577 template<class T> 01578 T dot(const TMat<T>& m1, const TMat<T>& m2) 01579 { 01580 #ifdef BOUNDCHECK 01581 if(m1.size()!=m2.size()) 01582 PLERROR("In T operator*(const TMat<T>& m1, const TVec<T>& vec2) (dot product) the 2 matrices must have the same number of elements."); 01583 #endif 01584 01585 T res = 0; 01586 if (m1.size() > 0 && m2.size() > 0) { 01587 T* v1 = m1.data(); 01588 T* v2 = m2.data(); 01589 if (m1.isCompact() && m2.isCompact()) 01590 for(int i=0; i<m1.size(); i++) 01591 res += v1[i]*v2[i]; 01592 else 01593 { 01594 TMatElementIterator<T> p1 = m1.begin(); 01595 TMatElementIterator<T> p2 = m2.begin(); 01596 for (int i=0; i<m1.size(); i++,++p1,++p2) 01597 res += *p1 * *p2; 01598 } 01599 } 01600 return res; 01601 } 01602 01603 template<class T> 01604 TVec<T> operator-(const TVec<T>& v1, const TVec<T>& v2) 01605 { 01606 if (v1.length() != v2.length()) 01607 PLERROR("TVec<T> - TVec<T>: different lengths %d and %d", 01608 v1.length(), v2.length()); 01609 TVec<T> v(v1.length()); 01610 v << v1; 01611 v-=v2; 01612 return v; 01613 } 01614 01615 template<class T> 01616 TVec<T> operator-(T v1, const TVec<T>& v2) 01617 { 01618 TVec<T> v(v2.length()); 01619 v = -v2; 01620 v += v1; 01621 return v; 01622 } 01623 01624 template<class T> 01625 TVec<T> operator-(const TVec<T>& v1, T v2) 01626 { 01627 TVec<T> v(v1.length()); 01628 substract(v1,v2,v); 01629 return v; 01630 } 01631 01632 template<class T> 01633 TVec<T> operator+(const TVec<T>& v1, const TVec<T>& v2) 01634 { 01635 if (v1.length() != v2.length()) 01636 PLERROR("TVec<T> + TVec<T>: different lengths %d and %d", 01637 v1.length(), v2.length()); 01638 TVec<T> v(v1.length()); 01639 v << v1; 01640 v+=v2; 01641 return v; 01642 } 01643 01644 template<class T> 01645 TVec<T> operator+(T v1, const TVec<T>& v2) 01646 { 01647 TVec<T> v(v2.length()); 01648 add(v2,v1,v); 01649 return v; 01650 } 01651 01652 template<class T> 01653 TVec<T> operator+(const TVec<T>& v1, T v2) 01654 { 01655 TVec<T> v(v1.length()); 01656 add(v1,v2,v); 01657 return v; 01658 } 01659 01660 template<class T> 01661 TVec<T> operator%(const TVec<T>& v1, const TVec<T>& v2) 01662 { 01663 if (v1.length() != v2.length()) 01664 PLERROR("TVec<T> + TVec<T>: different lengths %d and %d", 01665 v1.length(), v2.length()); 01666 TVec<T> v(v1.length()); 01667 v << v1; 01668 v*=v2; 01669 return v; 01670 } 01671 01672 template<class T> 01673 TVec<T> operator*(T scalar, const TVec<T>& v) 01674 { 01675 TVec<T> result(v.length()); 01676 multiply(v,scalar,result); 01677 return result; 01678 } 01679 01680 template<class T> 01681 TVec<T> operator*(const TVec<T>& v1, T v2) 01682 { 01683 TVec<T> v(v1.length()); 01684 multiply(v1,v2,v); 01685 return v; 01686 } 01687 01688 template<class T> 01689 TVec<T> operator/(const TVec<T>& v1, const TVec<T>& v2) 01690 { 01691 if (v1.length() != v2.length()) 01692 PLERROR("TVec<T> + TVec<T>: different lengths %d and %d", 01693 v1.length(), v2.length()); 01694 TVec<T> v(v1.length()); 01695 v << v1; 01696 v/=v2; 01697 return v; 01698 } 01699 01700 template<class T> 01701 TVec<T> operator/(T v1, const TVec<T>& v2) 01702 { 01703 int n=v2.length(); 01704 TVec<T> v(n); 01705 if (v2.size() > 0) { 01706 T* s2=v2.data(); 01707 T* d=v.data(); 01708 for (int i=0;i<n;i++) 01709 d[i] = v1/s2[i]; 01710 } 01711 return v; 01712 } 01713 01714 // norman: changed to unharmful declaration (see below old style) 01715 template<class T1, class T2> 01716 TVec<T1> operator/(const TVec<T1>& v1, T2 scalar) 01717 { 01718 TVec<T1> v(v1.length()); 01719 multiply(v1,T1(1.0)/(T1)scalar,v); 01720 return v; 01721 } 01722 01723 // norman: harmful declarations 01724 // Replaced with a better declaration above 01725 //template<class T> 01726 //TVec<T> operator/(const TVec<T>& v1, T scalar) 01727 //{ 01728 // TVec<T> v(v1.length()); 01729 // multiply(v1,T(1.0)/scalar,v); 01730 // return v; 01731 //} 01732 01733 // norman: This will cause problems if T = int (recursive declaration) 01734 // Replaced with a better declaration above 01735 //template<class T> 01736 //TVec<T> operator/(const TVec<T>& v1, int scalar) 01737 //{ return v1/T(scalar); } 01738 01739 template<class T> 01740 T logadd(const TVec<T>& vec) 01741 { 01742 int l = vec.length(); 01743 if(l==0) 01744 return LOG_INIT; 01745 01746 T *p_x = vec.data(); 01747 T sum = *p_x++; 01748 for (int i=1; i<l; i++, p_x++) 01749 sum = logadd(sum, *p_x); 01750 return sum; 01751 } 01752 01753 template<class T> 01754 T output_margin(const TVec<T>& class_scores, int correct_class) 01755 { 01756 T maxother = -FLT_MAX; 01757 for(int i=0; i<class_scores.length(); i++) 01758 { 01759 if(i!=correct_class && class_scores[i]>maxother) 01760 maxother = class_scores[i]; 01761 } 01762 return class_scores[correct_class]-maxother; 01763 } 01764 01765 template<class T> 01766 void fill_one_hot(const TVec<T>& vec, int hotpos, T coldvalue, T hotvalue) 01767 { 01768 #ifdef BOUNDCHECK 01769 if(!vec) 01770 PLERROR("In fill_one_hot given vec must have the correct size"); 01771 if(hotpos<0 || (vec.length()==1 && hotpos>1) || (vec.length()>1 && hotpos>=vec.length())) 01772 PLERROR("In fill_one_hot given hotpos out of vec range"); 01773 #endif 01774 if(vec.length()==1) 01775 vec[0] = (hotpos==0 ?coldvalue :hotvalue); 01776 else 01777 { 01778 vec.fill(coldvalue); 01779 vec[hotpos] = hotvalue; 01780 } 01781 } 01782 01783 template<class T> 01784 TVec<T> one_hot(int length, int hotpos, T coldvalue, T hotvalue) 01785 { 01786 TVec<T> result(length); 01787 fill_one_hot(result, hotpos, coldvalue, hotvalue); 01788 return result; 01789 } 01790 01791 template<class T> 01792 TVec<T> square(const TVec<T>& vec) 01793 { 01794 int n = vec.length(); 01795 TVec<T> result(n); 01796 square(result,vec); 01797 return result; 01798 } 01799 01800 template<class T> 01801 void square(TVec<T>& result, const TVec<T>& vec) 01802 { 01803 #ifdef BOUNDCHECK 01804 if (result.size() != vec.size()) 01805 PLERROR("In square, 'result' and 'vec' must have the same size"); 01806 #endif 01807 int n = vec.length(); 01808 if (n > 0) { 01809 T* v = vec.data(); 01810 T* r = result.data(); 01811 for(int i=0; i<n; i++) 01812 r[i] = v[i]*v[i]; 01813 } 01814 } 01815 01816 template<class T> 01817 TVec<T> squareroot(const TVec<T>& vec) 01818 { 01819 int n = vec.length(); 01820 TVec<T> result(n); 01821 if (n > 0) { 01822 T* v = vec.data(); 01823 T* r = result.data(); 01824 for(int i=0; i<n; i++) 01825 r[i] = sqrt(v[i]); 01826 } 01827 return result; 01828 } 01829 01832 template<class T> 01833 TVec<T> remove_missing(const TVec<T>& vec) 01834 { 01835 int n = vec.length(); 01836 int n_non_missing = 0; 01837 TVec<T> result(n); 01838 if (n > 0) { 01839 T* v = vec.data(); 01840 T* r = result.data(); 01841 for(int i=0; i<n; i++) { 01842 if (!is_missing(v[i])) 01843 r[n_non_missing++] = v[i]; 01844 } 01845 result.resize(n_non_missing); 01846 } 01847 return result; 01848 } 01849 01851 template<class T> 01852 void remove_missing_inplace(TVec<T>& v) 01853 { 01854 int n_non_missing=0; 01855 int next_non_missing=1; 01856 T* d = v.data(); 01857 for(;;) 01858 { 01859 while(n_non_missing<v.length()&&!is_missing(d[n_non_missing])) 01860 { 01861 n_non_missing++;next_non_missing++; 01862 } 01863 if(n_non_missing>=v.length()) 01864 return; 01865 while(next_non_missing<v.length()&&is_missing(d[next_non_missing])) 01866 next_non_missing++; 01867 if(next_non_missing>=v.length()) 01868 { 01869 v.resize(n_non_missing); 01870 return; 01871 } 01872 else 01873 { 01874 pl_swap(d[n_non_missing],d[next_non_missing]); 01875 } 01876 } 01877 } 01878 01881 template<class T, class U, class V> 01882 TVec<U> apply(const TVec<T>& vec, U (*func)(V)) 01883 { 01884 TVec<U> destination(vec.length()); 01885 apply(vec,destination,func); 01886 return destination; 01887 } 01888 01890 template<class T, class U> 01891 void apply(const TVec<T>& source, TVec<U>& destination, U (*func)(T)) 01892 { 01893 int n=source.length(); 01894 if (n!=destination.length()) 01895 PLERROR("apply: source(%d) and destination(%d) TVec<T>'s must have same length", 01896 n,destination.length()); 01897 if (n > 0) { 01898 T* s = source.data(); 01899 U* d = destination.data(); 01900 for(int i=0; i<n; i++) 01901 d[i]=func(s[i]); 01902 } 01903 } 01904 01907 template<class T, class U, class V> 01908 void apply(const TVec<T>& src1,const TVec<U>& src2, TVec<V>& dest, 01909 V (*func)(T,U)) 01910 { 01911 int n=src1.length(); 01912 if (n!=dest.length() || n!=src2.length()) 01913 PLERROR("apply: src1, src2 and destination TVec<T>'s must have same length"); 01914 if (n > 0) { 01915 T* s1 = src1.data(); 01916 U* s2 = src2.data(); 01917 V* d = dest.data(); 01918 for(int i=0; i<n; i++) 01919 d[i]=func(s1[i],s2[i]); 01920 } 01921 } 01922 01923 01924 // Efficient mathematical operations (without memory allocation) 01925 01926 // destination[i] = source1[i]*source2[i] 01927 template<class T> 01928 void multiply(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 01929 { 01930 int n=source1.length(); 01931 if (n!=source2.length()) 01932 PLERROR("multiply: two sources (l=%d and %d) must have same length", 01933 n,source2.length()); 01934 if (n!=destination.length()) 01935 destination.resize(n); 01936 if (n > 0) { 01937 T* s1=source1.data(); 01938 T* s2=source2.data(); 01939 T* d=destination.data(); 01940 for (int i=0;i<n;i++) 01941 d[i] = s1[i]*s2[i]; 01942 } 01943 } 01944 01945 // destination[i] = source1[i] + source2[i]*source3 01946 template<class T> 01947 void multiplyAdd(const TVec<T>& source1, const TVec<T>& source2, 01948 T source3, TVec<T>& destination) 01949 { 01950 int n=source1.length(); 01951 if (n!=source2.length()) 01952 PLERROR("multiply: two sources (l=%d and %d) must have same length", 01953 n,source2.length()); 01954 if (n!=destination.length()) 01955 destination.resize(n); 01956 if (n > 0) { 01957 T* s1=source1.data(); 01958 T* s2=source2.data(); 01959 T* d=destination.data(); 01960 for (int i=0;i<n;i++) 01961 d[i] = s1[i]+s2[i]*source3; 01962 } 01963 } 01964 01965 // destination[i] = a*destination[i] + b*source[i] 01966 template<class T> 01967 void multiplyScaledAdd(const TVec<T>& source, T a, T b, const TVec<T>& destination) 01968 { 01969 int n=source.length(); 01970 if (n!=destination.length()) 01971 PLERROR("multiply: source and destination (l=%d and %d) must have same length", 01972 n,destination.length()); 01973 if (n > 0) { 01974 T* s=source.data(); 01975 T* d=destination.data(); 01976 for (int i=0;i<n;i++) 01977 d[i] = a*d[i] + b*s[i]; 01978 } 01979 } 01980 01981 // destination[i,j] = a*destination[i,j] + b*source[i,j] 01982 template<class T> 01983 void multiplyScaledAdd(const TMat<T>& source, T a, T b, const TMat<T>& destination) 01984 { 01985 int n=source.length(); 01986 int m=source.width(); 01987 if (n!=destination.length() || m!=destination.width()) 01988 PLERROR("multiply: source and destination must have same dimensions"); 01989 if (n > 0) { 01990 int sm=source.mod(); 01991 int dm=destination.mod(); 01992 T* s=source.data(); 01993 T* d=destination.data(); 01994 for (int i=0;i<n;i++,s+=sm,d+=dm) 01995 for (int j=0;j<m;j++) 01996 d[j] = a*d[j] + b*s[j]; 01997 } 01998 } 01999 02000 // destination[i] = source1[i]+source2[i] 02001 template<class T> 02002 void add(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 02003 { 02004 int n=source1.length(); 02005 if (n!=source2.length()) 02006 PLERROR("add: two sources (l=%d and %d) must have same length", 02007 n,source2.length()); 02008 if (n!=destination.length()) 02009 destination.resize(n); 02010 if (n > 0) { 02011 T* s1=source1.data(); 02012 T* s2=source2.data(); 02013 T* d=destination.data(); 02014 for (int i=0;i<n;i++) 02015 d[i] = s1[i]+s2[i]; 02016 } 02017 } 02018 02019 // destination[i] = source1[i]+source2 02020 template<class T> 02021 void add(const TVec<T>& source1, T source2, TVec<T>& destination) 02022 { 02023 int n=source1.length(); 02024 if (n!=destination.length()) 02025 destination.resize(n); 02026 if (n > 0) { 02027 T* s1=source1.data(); 02028 T* d=destination.data(); 02029 for (int i=0;i<n;i++) 02030 d[i] = s1[i]+source2; 02031 } 02032 } 02033 02034 template<class T> 02035 inline void substract(const TVec<T>& source1, T source2, TVec<T>& destination) 02036 { add(source1,-source2,destination); } 02037 02038 // destination[i] = source1[i]-source2[i] 02039 template<class T> 02040 void substract(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 02041 { 02042 int n=source1.length(); 02043 if (n!=source2.length()) 02044 PLERROR("substract: two sources (l=%d and %d) must have same length", 02045 n,source2.length()); 02046 if (n!=destination.length()) 02047 destination.resize(n); 02048 if (n > 0) { 02049 T* s1=source1.data(); 02050 T* s2=source2.data(); 02051 T* d=destination.data(); 02052 for (int i=0;i<n;i++) 02053 d[i] = s1[i]-s2[i]; 02054 } 02055 } 02056 02057 // destination[i] += source1[i]-source2[i] 02058 template<class T> 02059 void substractAcc(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 02060 { 02061 int n=source1.length(); 02062 if (n!=source2.length()) 02063 PLERROR("substract: two sources (l=%d and %d) must have same length", 02064 n,source2.length()); 02065 if (n!=destination.length()) 02066 destination.resize(n); 02067 if (n > 0) { 02068 T* s1=source1.data(); 02069 T* s2=source2.data(); 02070 T* d=destination.data(); 02071 for (int i=0;i<n;i++) 02072 d[i] += s1[i]-s2[i]; 02073 } 02074 } 02075 02076 // destination[i] = source1-source2[i] 02077 template<class T> 02078 void substract(T source1, const TVec<T>& source2, TVec<T>& destination) 02079 { 02080 int n=source2.length(); 02081 if (n!=destination.length()) 02082 destination.resize(n); 02083 if (n > 0) { 02084 T* s2=source2.data(); 02085 T* d=destination.data(); 02086 for (int i=0;i<n;i++) 02087 d[i] = source1-s2[i]; 02088 } 02089 } 02090 02091 template<class T> 02092 inline void divide(const TVec<T>& source1, T source2, TVec<T>& destination) 02093 { multiply(source1,1.0/source2,destination); } 02094 02095 // destination[i] = source1[i]/source2[i] 02096 template<class T> 02097 void divide(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 02098 { 02099 int n=source1.length(); 02100 if (n!=source2.length()) 02101 PLERROR("divide: two sources (l=%d and %d) must have same length", 02102 n,source2.length()); 02103 if (n!=destination.length()) 02104 destination.resize(n); 02105 if (n > 0) { 02106 T* s1=source1.data(); 02107 T* s2=source2.data(); 02108 T* d=destination.data(); 02109 for (int i=0;i<n;i++) 02110 d[i] = s1[i]/s2[i]; 02111 } 02112 } 02113 02114 // destination[i] = source1/source2[i] 02115 template<class T> 02116 void divide(T source1, const TVec<T>& source2, TVec<T>& destination) 02117 { 02118 int n=source2.length(); 02119 if (n!=destination.length()) 02120 destination.resize(n); 02121 if (n > 0) { 02122 T* s2=source2.data(); 02123 T* d=destination.data(); 02124 for (int i=0;i<n;i++) 02125 d[i] = source1/s2[i]; 02126 } 02127 } 02128 02129 // destination[i] = max(source1[i],source2[i]) 02130 template<class T> 02131 void max(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 02132 { 02133 int n=source1.length(); 02134 if (n!=source2.length()) 02135 PLERROR("max: two sources (l=%d and %d) must have same length", 02136 n,source2.length()); 02137 if (n!=destination.length()) 02138 destination.resize(n); 02139 if (n > 0) { 02140 T* s1=source1.data(); 02141 T* s2=source2.data(); 02142 T* d=destination.data(); 02143 for (int i=0;i<n;i++) 02144 d[i] = MAX(s1[i],s2[i]); 02145 } 02146 } 02147 02148 // destination[i] = max(source1[i],source2) 02149 template<class T> 02150 void max(const TVec<T>& source1, T source2, TVec<T>& destination) 02151 { 02152 int n=source1.length(); 02153 if (n!=destination.length()) 02154 destination.resize(n); 02155 if (n > 0) { 02156 T* s1=source1.data(); 02157 T* d=destination.data(); 02158 for (int i=0;i<n;i++) 02159 d[i] = MAX(s1[i],source2); 02160 } 02161 } 02162 02163 02164 // destination[i] = min(source1[i],source2[i]) 02165 template<class T> 02166 void min(const TVec<T>& source1, const TVec<T>& source2, TVec<T>& destination) 02167 { 02168 int n=source1.length(); 02169 if (n!=source2.length()) 02170 PLERROR("min: two sources (l=%d and %d) must have same length", 02171 n,source2.length()); 02172 if (n!=destination.length()) 02173 destination.resize(n); 02174 if (n > 0) { 02175 T* s1=source1.data(); 02176 T* s2=source2.data(); 02177 T* d=destination.data(); 02178 for (int i=0;i<n;i++) 02179 d[i] = MIN(s1[i],s2[i]); 02180 } 02181 } 02182 02183 // destination[i] = min(source1[i],source2) 02184 template<class T> 02185 void min(const TVec<T>& source1, T source2, TVec<T>& destination) 02186 { 02187 int n=source1.length(); 02188 if (n!=destination.length()) 02189 destination.resize(n); 02190 if (n > 0) { 02191 T* s1=source1.data(); 02192 T* d=destination.data(); 02193 for (int i=0;i<n;i++) 02194 d[i] = MIN(s1[i],source2); 02195 } 02196 } 02197 02198 02199 template<class T> 02200 TVec<T> softmax(const TVec<T>& x) 02201 { 02202 TVec<T> y(x.length()); 02203 softmax(x,y); 02204 return y; 02205 } 02206 02207 template<class T> 02208 void tanh(const TVec<T>& x, TVec<T>& y) 02209 { 02210 int n = x.length(); 02211 #ifdef BOUNDCHECK 02212 if (y.length()!=n) 02213 PLERROR("tanh(TVec<T>,TVec<T>), second argument of length %d, first of length %d, should be =", 02214 n,y.length()); 02215 #endif 02216 if (n>0) 02217 { 02218 T* yp = y.data(); 02219 T* xp = x.data(); 02220 for (int i=0;i<n;i++) 02221 yp[i] = tanh(xp[i]); 02222 } 02223 } 02224 02225 template<class T> 02226 TVec<T> 02227 exp(const TVec<T>& vec) 02228 { 02229 TVec<T> res( vec.length() ); 02230 exp( vec, res ); 02231 return res; 02232 } 02233 02234 // return indices of non-zero elements 02235 template<class T> 02236 TVec<T> nonZeroIndices(TVec<T> v) 02237 { 02238 int n=v.length(); 02239 if (!n) 02240 return TVec<T>(); 02241 TVec<T> indices(n); 02242 int ni=0; 02243 T* val = v.data(); 02244 T* indx= indices.data(); 02245 for (int i=0;i<n;i++) 02246 if (val[i]!=0) 02247 indx[ni++]=i; 02248 indices.resize(ni); 02249 return indices; 02250 } 02251 02252 // return indices of non-zero elements 02253 template<class T> 02254 TVec<T> nonZeroIndices(TVec<bool> v) 02255 { 02256 int n=v.length(); 02257 if (!n) 02258 return TVec<T>(); 02259 TVec<T> indices(n); 02260 int ni=0; 02261 bool* val = v.data(); 02262 T* indx= indices.data(); 02263 for (int i=0;i<n;i++) 02264 if (val[i]) 02265 indx[ni++]=i; 02266 indices.resize(ni); 02267 return indices; 02268 } 02269 02270 // Set the complement indices, i.e. if 0<=i<n is not an element 02271 // of the indices vector it is put in the complement_indices vector. 02272 template<class T> 02273 void complement_indices(TVec<T>& indices, int n, 02274 TVec<T>& complement_indices, 02275 TVec<T>& buffer) 02276 { 02277 int ni=indices.length(); 02278 T* ind = indices.data(); 02279 T* cind = complement_indices.data(); 02280 buffer.resize(n); 02281 buffer.fill(0); 02282 T* buf=buffer.data(); 02283 for (int i=0;i<ni;i++) 02284 buf[(int)ind[i]]=1.0; 02285 for (int i=0,j=0;i<n;i++) 02286 if (buf[i]==0.0) 02287 cind[j++]=i; 02288 } 02289 02290 // dest[i] = 1 if src[i]==v, 0 otherwise 02291 template<class T> 02292 void equals(const TVec<T>& src, T v, TVec<T>& dest) 02293 { 02294 int n=src.length(); 02295 #ifdef BOUNDCHECK 02296 if (n!=dest.length()) 02297 PLERROR("equals(TVec<T>(%d),T,TVec<T>(%d)) args of unequal lengths", 02298 n,dest.length()); 02299 #endif 02300 if (n > 0) { 02301 T* s=src.data(); 02302 T* d=dest.data(); 02303 for (int i=0;i<n;i++) 02304 if (s[i]==v) d[i]=1.0; else d[i]=0.0; 02305 } 02306 } 02307 02308 // dest[i] = 1 if first[i] > second[i], 0 otherwise 02309 template<class T> 02310 void isLargerThan(const TVec<T>& first, const TVec<T>& second, TVec<T>& dest) 02311 { 02312 int n=first.length(); 02313 if(n!=second.length() || n!=dest.length()) 02314 PLERROR("isLargerThan(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal length", 02315 n, second.length(), dest.length()); 02316 if (n > 0) { 02317 T* f=first.data(); 02318 T* s=second.data(); 02319 T* d=dest.data(); 02320 for (int i=0; i<n; i++) 02321 d[i] = f[i] > s[i]; 02322 } 02323 } 02324 02325 // dest[i] = 1 if first[i] >= second[i], 0 otherwise 02326 template<class T> 02327 void isLargerThanOrEqualTo(const TVec<T>& first, const TVec<T>& second, TVec<T>& dest) 02328 { 02329 int n=first.length(); 02330 if(n!=second.length() || n!=dest.length()) 02331 PLERROR("isLargerThan(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal length", 02332 n, second.length(), dest.length()); 02333 if (n > 0) { 02334 T* f=first.data(); 02335 T* s=second.data(); 02336 T* d=dest.data(); 02337 for (int i=0; i<n; i++) 02338 d[i] = f[i] >= s[i]; 02339 } 02340 } 02341 02342 // dest[i] = 1 if first[i] < second[i], 0 otherwise 02343 template<class T> 02344 void isSmallerThan(const TVec<T>& first, const TVec<T>& second, TVec<T>& dest) 02345 { 02346 int n=first.length(); 02347 if(n!=second.length() || n!=dest.length()) 02348 PLERROR("isLargerThan(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal length", 02349 n, second.length(), dest.length()); 02350 if (n > 0) { 02351 T* f=first.data(); 02352 T* s=second.data(); 02353 T* d=dest.data(); 02354 for (int i=0; i<n; i++) 02355 d[i] = f[i] < s[i]; 02356 } 02357 } 02358 02359 // dest[i] = 1 if first[i] <= second[i], 0 otherwise 02360 template<class T> 02361 void isSmallerThanOrEqualTo(const TVec<T>& first, const TVec<T>& second, TVec<T>& dest) 02362 { 02363 int n=first.length(); 02364 if(n!=second.length() || n!=dest.length()) 02365 PLERROR("isLargerThan(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal length", 02366 n, second.length(), dest.length()); 02367 if (n > 0) { 02368 T* f=first.data(); 02369 T* s=second.data(); 02370 T* d=dest.data(); 02371 for (int i=0; i<n; i++) 02372 d[i] = f[i] <= s[i]; 02373 } 02374 } 02375 02376 // dest[i] = if_vec[i] ? then_vec[i] : else_vec[i]; 02377 template<class T> 02378 void ifThenElse(const TVec<T>& if_vec, const TVec<T>& then_vec, 02379 const TVec<T>& else_vec, TVec<T>& dest) 02380 { 02381 int n=if_vec.length(); 02382 if (n!=then_vec.length() || n!=else_vec.length() || n!=dest.length()) 02383 PLERROR("ifThenElse(TVec<T>(%d), TVec<T>(%d), TVec<T>(%d), TVec<T>(%d)) args of unequal lengths", 02384 n, then_vec.length(), else_vec.length(), dest.length()); 02385 if (n > 0) { 02386 T* i_=if_vec.data(); 02387 T* t_=then_vec.data(); 02388 T* e_=else_vec.data(); 02389 T* d_=dest.data(); 02390 for (int i=0;i<n;i++) 02391 d_[i] = i_[i] ? t_[i] : e_[i]; 02392 } 02393 } 02394 02395 // returns the number of times that src[i] == value 02396 template<class T> 02397 int vec_counts(const TVec<T>& src, T value) 02398 { 02399 int len = src.length(); 02400 int n = 0; 02401 if (len > 0) { 02402 T *p = src.data(); 02403 for (int i=0; i<len; i++, p++) 02404 if (*p == value) 02405 n++; 02406 } 02407 return n; 02408 } 02409 02410 // returns the position of f in src (-1 if f is not found) 02411 template<class T> 02412 int vec_find(const TVec<T>& src, T f) 02413 { 02414 int len = src.length(); 02415 if (len > 0) { 02416 T *p = src.data(); 02417 for (int i=0; i<len; i++, p++) 02418 if (*p == f) 02419 return(i); 02420 } 02421 return -1; 02422 } 02423 02424 02425 template<class T> 02426 T estimatedCumProb(T x, TVec<T> bins) 02427 { 02428 const int nbins = bins.length()-1; 02429 if (nbins<1) PLERROR("estimatedCumProb:: there should be at least two elements in the bins vector"); 02430 // +0.5 because we allocate mass 0.25 at the left and 0.25 at the right of the interval (bins(0),bins(nbins)) 02431 const T one_over_nbins = 1.0/(T)(nbins+0.5); 02432 02433 int k = binary_search(bins, x); 02434 02435 if (k == -1) 02436 return 0.25*one_over_nbins; 02437 else if (k == nbins-1) 02438 return 1.0 - 0.25*one_over_nbins; 02439 else if (bins[k] != bins[k+1]) 02440 return one_over_nbins*(0.25 + k + (x-bins[k])/(bins[k+1]-bins[k])); 02441 else 02442 return one_over_nbins*(0.75 + k); 02443 } 02444 02445 // returns the index of the kth ordered element of v 02446 // (dumb algorithm, takes time in k*n ) 02447 template<class T> 02448 int positionOfkthOrderedElement(const TVec<T>& vec, int k) 02449 { 02450 #ifdef BOUNDCHECK 02451 if(k<0 || k>=vec.length()) 02452 PLERROR("In positionOfkthOrderedElement, k out of bounds"); 02453 #endif 02454 02455 T* v = vec.data(); 02456 02457 T minval = -FLT_MAX; 02458 int pos = -1; 02459 int l=0; 02460 02461 while(l<=k) 02462 { 02463 int nelements_equal_to_newminval = 0; 02464 T newminval = FLT_MAX; 02465 for(int i=0; i<vec.length(); i++) 02466 { 02467 if(v[i]>minval) 02468 { 02469 if(v[i]<newminval) 02470 { 02471 newminval = v[i]; 02472 nelements_equal_to_newminval = 1; 02473 pos = i; 02474 } 02475 else if(v[i]==newminval) 02476 nelements_equal_to_newminval++; 02477 } 02478 } 02479 l += nelements_equal_to_newminval; 02480 minval = newminval; 02481 } 02482 02483 return pos; 02484 } 02485 02488 template<class T> 02489 inline T kthOrderedElement(const TVec<T>& vec, int k) 02490 { return vec[positionOfkthOrderedElement(vec,k)]; } 02491 02493 template<class T> 02494 inline T median(const TVec<T>& vec) 02495 { 02496 if (vec.isEmpty()) 02497 PLERROR("In median - Cannot compute median of an empty vector"); 02498 return kthOrderedElement(vec, (vec.length()-1)/2); 02499 } 02500 02501 //-------------- These were previouslty methods of TVec ---------------------------------- 02502 02503 02506 template<class T> 02507 T selectAndOrder(const TVec<T>& vec, int pos) 02508 { 02509 if (pos<0 || pos>=vec.length()) PLERROR("Bad position (%d)", pos); 02510 02511 int l=0; 02512 int h=vec.length()-1; 02513 T* v = vec.data(); 02514 02515 while (l<h) 02516 { 02517 T p = v[(l+h)/2]; 02518 int x = l; 02519 int y = h; 02520 02521 do 02522 { 02523 while (v[x]<p) x++; 02524 while (p<v[y]) y--; 02525 if (x<=y) 02526 { 02527 PLearn::swap(v[x],v[y]); 02528 x++; 02529 y--; 02530 } 02531 } while (x<=y); 02532 02533 if (pos>=x) l=x; 02534 else h=x-1; 02535 } 02536 02537 return v[l]; 02538 } 02539 02544 template<class T> 02545 TVec<T> getQuantiles(const TVec<T>& vec, int q) 02546 { 02547 int l = vec.length(); 02548 T* v = vec.data(); 02549 TVec<T> w(q+1); 02550 T linvq = T(l)/q; 02551 for(int i=0;i<q;i++) w[i] = v[int(linvq*i)]; 02552 w[q]=v[l-1]; 02553 return w; 02554 } 02555 02558 template<class T> 02559 TVec<T> nonZero(const TVec<T>& vec) 02560 { 02561 T *v =vec.data(); 02562 int n=0; 02563 for(int i=0;i<vec.length(); i++) if (v[i]!=0) n++; 02564 TVec<T> w(n); 02565 int j=0; 02566 for(int i=0;i<vec.length(); i++) if (v[i]!=0) w[j++]=v[i]; 02567 return(w); 02568 } 02569 02572 template<class T> 02573 TVec<T> positiveValues(const TVec<T>& vec) 02574 { 02575 T *v =vec.data(); 02576 int n=0; 02577 for(int i=0;i<vec.length(); i++) if (v[i]>0) n++; 02578 TVec<T> w(n); 02579 int j=0; 02580 for(int i=0;i<vec.length(); i++) if (v[i]>0) w[j++]=v[i]; 02581 return(w); 02582 } 02583 02588 template<class T> 02589 int positionOfClosestElement(const TVec<T>& vec, const T& value, bool is_sorted_vec=false) 02590 { 02591 T* v = vec.data(); 02592 if (is_sorted_vec) // dichotomy search 02593 { 02594 int pos = binary_search(vec, value); 02595 if (pos == -1) return 0; 02596 else if (pos == vec.length()-1) return pos; 02597 T dist1 = fabs(v[pos]-value); 02598 T dist2 = fabs(v[pos+1]-value); 02599 if (dist1 <= dist2) return pos; 02600 else return pos+1; 02601 } 02602 else // linear search 02603 { 02604 int pos_of_closest = 0; 02605 T dist_to_closest = fabs(v[0]-value); 02606 for(int i=1; i<vec.length(); i++) 02607 { 02608 T dist = fabs(v[i]-value); 02609 if(dist<dist_to_closest) 02610 { 02611 pos_of_closest = i; 02612 dist_to_closest = dist; 02613 } 02614 } 02615 return pos_of_closest; 02616 } 02617 } 02618 02619 02626 template <class T> 02627 void projectOnOrthogonalSubspace(const TVec<T>& vec, const TMat<T>& orthonormal_subspace) 02628 { 02629 for (int i=0;i<orthonormal_subspace.length();i++) 02630 { 02631 TVec<T> vi = orthonormal_subspace(i); 02632 T dp = dot(vec,vi); 02633 multiplyAcc(vec, vi,-dp); 02634 } 02635 } 02636 02637 02639 template<class T> 02640 void multiplyAcc(const TVec<T>& vec, const TVec<T>& x, T scale) 02641 { 02642 int n=x.length(); 02643 if (vec.length()!=n) 02644 PLERROR("TVec::multiplyAcc this has length_=%d and x has length_=%d", vec.length(),n); 02645 T* p=vec.data(); 02646 T* xp=x.data(); 02647 for (int i=0;i<n;i++) 02648 *p++ += scale * *xp++; 02649 } 02650 02652 template<class T> 02653 void exponentialMovingAverageUpdate(const TVec<T>& vec, const TVec<T>& x, T alpha) 02654 { 02655 int n=x.length(); 02656 if (vec.length()!=n) 02657 PLERROR("TVec::exponentialMovingAverageUpdate length_=%d and x has length_=%d", 02658 vec.length(),n); 02659 T* p=vec.data(); 02660 T* xp=x.data(); 02661 T one_minus_alpha = 1-alpha; 02662 for (int i=0;i<n;i++) 02663 p[i] = one_minus_alpha*p[i] + alpha*xp[i]; 02664 } 02665 02667 template<class T> 02668 void exponentialMovingVarianceUpdate(const TVec<T>& vec, const TVec<T>& x, const TVec<T>& mu, T alpha) 02669 { 02670 int n=x.length(); 02671 if (vec.length()!=n || vec.length()!=mu.length()) 02672 PLERROR("TVec::exponentialVarianceAverageUpdate length_=%d and" 02673 "x has length_=%d, mu has length() %d", 02674 vec.length(),n,mu.length()); 02675 T* p=vec.data(); 02676 T* xp=x.data(); 02677 T* mp=mu.data(); 02678 T one_minus_alpha = 1-alpha; 02679 for (int i=0;i<n;i++) 02680 { 02681 T dif = (xp[i]-mp[i]); 02682 p[i] = one_minus_alpha*p[i] + alpha*dif*dif; 02683 } 02684 } 02685 02686 02688 template<class T> 02689 void exponentialMovingSquareUpdate(const TVec<T>& vec, const TVec<T>& x, T alpha) 02690 { 02691 int n=x.length(); 02692 if (vec.length()!=n) 02693 PLERROR("TVec::exponentialMovingAverageUpdate length_=%d and x has length_=%d", 02694 vec.length(),n); 02695 T* p=vec.data(); 02696 T* xp=x.data(); 02697 T one_minus_alpha = 1-alpha; 02698 for (int i=0;i<n;i++) 02699 { 02700 T xpi = xp[i]; 02701 p[i] = one_minus_alpha*p[i] + alpha*xpi*xpi; 02702 } 02703 } 02704 02706 template<class T> 02707 void multiplyAcc(const TVec<T>& vec, const TVec<T>& x, const TVec<T>& y) 02708 { 02709 int n=x.length(); 02710 if (vec.length()!=n || y.length()!=n) 02711 PLERROR("TVec::multiplyAcc, this+=x*y: length_=%d, x.length_=%d, y.length_=%d", 02712 vec.length(),n,y.length()); 02713 T* p=vec.data(); 02714 T* xp=x.data(); 02715 T* yp=y.data(); 02716 for (int i=0;i<n;i++) 02717 p[i] += xp[i] * yp[i]; 02718 } 02719 02721 template<class T> 02722 void squareMultiplyAcc(const TVec<T>& vec, const TVec<T>& x, T scale) 02723 { 02724 int n=x.length(); 02725 if (vec.length()!=n) 02726 PLERROR("TVec::squareMultiplyAcc this has length_=%d and x has length_=%d", vec.length(),n); 02727 T* p=vec.data(); 02728 T* xp=x.data(); 02729 for (int i=0;i<n;i++) 02730 { 02731 T xpi = xp[i]; 02732 p[i] += scale * xpi * xpi; 02733 } 02734 } 02735 02737 template<class T> 02738 void squareAcc(const TVec<T>& vec, const TVec<T>& x) 02739 { 02740 int n=x.length(); 02741 if (vec.length()!=n) 02742 PLERROR("TVec::squareAcc this has length_=%d and x has length_=%d", vec.length(),n); 02743 T* p=vec.data(); 02744 T* xp=x.data(); 02745 for (int i=0;i<n;i++) 02746 { 02747 T xpi = xp[i]; 02748 p[i] += xpi * xpi; 02749 } 02750 } 02751 02753 template<class T> 02754 void squareSubtract(const TVec<T>& vec, const TVec<T>& x) 02755 { 02756 int n=x.length(); 02757 if (vec.length()!=n) 02758 PLERROR("TVec::squareDiff this has length_=%d and x has length_=%d", vec.length(),n); 02759 T* p=vec.data(); 02760 T* xp=x.data(); 02761 for (int i=0;i<n;i++) 02762 { 02763 T xpi = xp[i]; 02764 p[i] -= xpi * xpi; 02765 } 02766 } 02767 02769 template<class T> 02770 void diffSquareMultiplyAcc(const TVec<T>& vec, const TVec<T>& x, const TVec<T>& y, T scale) 02771 { 02772 int n=x.length(); 02773 if (vec.length()!=n || y.length()!=n) 02774 PLERROR("TVec::diffSquareMultiplyAcc this.length_=%d, x.length_=%d, y.length_=%d", 02775 vec.length(),n,y.length()); 02776 T* p=vec.data(); 02777 T* xp=x.data(); 02778 T* yp=y.data(); 02779 for (int i=0;i<n;i++) 02780 { 02781 T diff = xp[i]-yp[i]; 02782 p[i] += scale * diff * diff; 02783 } 02784 } 02785 02787 template<class T> 02788 void diffSquareMultiplyScaledAcc(const TVec<T>& vec, const TVec<T>& x, const TVec<T>& y, T fact1, T fact2) 02789 { 02790 int n=x.length(); 02791 if (vec.length()!=n || y.length()!=n) 02792 PLERROR("TVec::diffSquareMultiplyAcc this.length_=%d, x.length_=%d, y.length_=%d", 02793 vec.length(),n,y.length()); 02794 T* p=vec.data(); 02795 T* xp=x.data(); 02796 T* yp=y.data(); 02797 for (int i=0;i<n;i++) 02798 { 02799 T diff = xp[i]-yp[i]; 02800 p[i] = fact1 * p[i] + fact2 * diff * diff; 02801 } 02802 } 02803 02805 template <class T> 02806 void product(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v) 02807 { 02808 int l = m.length(); 02809 int w = m.width(); 02810 #ifdef BOUNDCHECK 02811 if (l!=result.length() || w!=v.length()) 02812 PLERROR("product(TVec, TMat, TVec), incompatible arguments:\n" 02813 "%d <- %dx%d times %d", 02814 result.length(), l, w, v.length()); 02815 #endif 02816 02817 if (m.isEmpty() || v.isEmpty() || result.isEmpty()) 02818 { 02819 // Size zero: no need to bother computing anything. 02820 // In such a case, the result of the matrix-vector multiplication, if 02821 // not empty, is necessarily zero, since R^0 = {0}. 02822 if (!result.isEmpty()) 02823 result.clear(); 02824 return; 02825 } 02826 02827 T *rp = result.data(); 02828 T *vp = v.data(); 02829 for (int i=0;i<l;i++) 02830 { 02831 const T* mi = m[i]; 02832 T s = 0; 02833 for (int j=0;j<w;j++) 02834 s += mi[j] * vp[j]; 02835 rp[i] = s; 02836 } 02837 } 02838 02840 template <class T> 02841 void productAcc(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v) 02842 { 02843 int l = m.length(); 02844 int w = m.width(); 02845 #ifdef BOUNDCHECK 02846 if (l!=result.length() || w!=v.length()) 02847 PLERROR("productAcc(TVec, TMat, TVec), incompatible arguments:\n" 02848 "%d <- %dx%d times %d", 02849 result.length(), l, w, v.length()); 02850 #endif 02851 02852 if (m.isEmpty() || v.isEmpty() || result.isEmpty()) 02853 { 02854 // Size zero: no need to bother computing anything. 02855 // In such a case, the result of the matrix-vector multiplication, if 02856 // not empty, is necessarily zero, since R^0 = {0}. 02857 return; 02858 } 02859 02860 T* rp = result.data(); 02861 T* mp = m.data(); 02862 T* vdata = v.data(); 02863 int deltam = m.mod()-m.width(); 02864 for (int i=0;i<l;i++) 02865 { 02866 T *vp = vdata; 02867 T s = *rp; 02868 for (int j=0;j<w;j++) 02869 s += *mp++ * *vp++; 02870 *rp++ = s; 02871 mp += deltam; 02872 } 02873 } 02874 02877 template <class T> 02878 void productScaleAcc(const TVec<T>& result, const TMat<T>& m, bool transpose_m, 02879 const TVec<T>& v, T alpha, T beta) 02880 { 02881 if (transpose_m) 02882 transposeProductScaleAcc(result, m, v, alpha, beta); 02883 else 02884 productScaleAcc(result, m, v, alpha, beta); 02885 } 02886 02888 template <class T> 02889 void productScaleAcc(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v, 02890 T alpha, T beta) 02891 { 02892 int l = m.length(); 02893 int w = m.width(); 02894 #ifdef BOUNDCHECK 02895 if (l!=result.length() || w!=v.length()) 02896 PLERROR("productScaleAcc(TVec, TMat, TVec), incompatible arguments:\n" 02897 "%d <- %dx%d times %d", 02898 result.length(), l, w, v.length()); 02899 #endif 02900 02901 if (m.isEmpty() || v.isEmpty() || result.isEmpty()) 02902 { 02903 // Size zero: no need to bother computing anything. 02904 // In such a case, the result of the matrix-vector multiplication, if 02905 // not empty, is necessarily zero, since R^0 = {0}. 02906 if (!result.isEmpty()) 02907 result *= beta; 02908 return; 02909 } 02910 02911 T* rp = result.data(); 02912 T* mp = m.data(); 02913 T* vdata = v.data(); 02914 int deltam = m.mod()-m.width(); 02915 for (int i=0;i<l;i++) 02916 { 02917 T *vp = vdata; 02918 T s = 0; 02919 for (int j=0;j<w;j++) 02920 s += *mp++ * *vp++; 02921 *rp = alpha * s + beta * (*rp); 02922 ++rp; 02923 mp += deltam; 02924 } 02925 } 02926 02930 template <class T> 02931 void transposeProduct(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v) 02932 { 02933 Profiler::pl_profile_start("transposeProduct T"); 02934 int l=m.length(); 02935 #ifdef BOUNDCHECK 02936 int w=m.width(); 02937 if (l!=v.length() || w!=result.length()) 02938 PLERROR("transposeProduct(TVec, TMat, TVec), incompatible arguments:\n" 02939 "%d <- %dx%d' times %d", 02940 result.length(), l, w, v.length()); 02941 #endif 02942 02943 if (m.isEmpty() || v.isEmpty() || result.isEmpty()) 02944 { 02945 // Size zero: no need to bother computing anything. 02946 // In such a case, the result of the matrix-vector multiplication, if 02947 // not empty, is necessarily zero, since R^0 = {0}. 02948 if (!result.isEmpty()) 02949 result.clear(); 02950 Profiler::pl_profile_end("transposeProduct T"); 02951 return; 02952 } 02953 02954 T *rp = result.data(); 02955 T *vp = v.data(); 02956 result.clear(); 02957 for (int j=0;j<l;j++) 02958 { 02959 const T* mj = m[j]; 02960 T vj = vp[j]; 02961 for (int i=0;i<result.length();i++) 02962 rp[i] += mj[i] * vj; 02963 } 02964 Profiler::pl_profile_end("transposeProduct T"); 02965 } 02966 02968 template <class T> 02969 void transposeProductAcc(const TVec<T>& result, const TMat<T>& m, 02970 const TVec<T>& v) 02971 { 02972 int l=m.length(); 02973 int w=m.width(); 02974 #ifdef BOUNDCHECK 02975 if (l!=v.length() || w!=result.length()) 02976 PLERROR("transposeProductAcc(TVec, TMat, TVec), incompatible arguments" 02977 ":\n" 02978 "%dx%d' times %d -> %d", 02979 result.length(), l, w, v.length()); 02980 #endif 02981 02982 if (m.isEmpty() || v.isEmpty() || result.isEmpty()) 02983 { 02984 // Size zero: no need to bother computing anything. 02985 // In such a case, the result of the matrix-vector multiplication, if 02986 // not empty, is necessarily zero, since R^0 = {0}. 02987 return; 02988 } 02989 02990 T* rdata = result.data(); 02991 T* vp = v.data(); 02992 T* mp = m.data(); 02993 int deltam = m.mod()-m.width(); 02994 for (int j=0;j<l;j++) 02995 { 02996 T vj = *vp++; 02997 02998 /* 02999 T* rp = rdata; 03000 for (int i=0;i<w;i++) 03001 *rp++ += vj * *mp++; 03002 mp += deltam; 03003 */ 03004 03005 if(vj!=0) 03006 { 03007 if(vj==1) 03008 { 03009 T* rp = rdata; 03010 for (int i=0;i<w;i++) 03011 *rp++ += *mp++; 03012 mp += deltam; 03013 } 03014 else 03015 { 03016 T* rp = rdata; 03017 for (int i=0;i<w;i++) 03018 *rp++ += vj * *mp++; 03019 mp += deltam; 03020 } 03021 } 03022 else mp += w + deltam; 03023 } 03024 } 03025 03027 template <class T> 03028 void transposeProductScaleAcc(const TVec<T>& result, const TMat<T>& m, 03029 const TVec<T>& v, T alpha, T beta) 03030 { 03031 int l=m.length(); 03032 int w=m.width(); 03033 #ifdef BOUNDCHECK 03034 if (l!=v.length() || w!=result.length()) 03035 PLERROR("transposeProductScaleAcc(TVec, TMat, TVec), incompatible" 03036 " arguments:\n" 03037 "%d <- %dx%d' times %d", 03038 result.length(), l, w, v.length()); 03039 #endif 03040 03041 if (m.isEmpty() || v.isEmpty() || result.isEmpty()) 03042 { 03043 // Size zero: no need to bother computing anything. 03044 // In such a case, the result of the matrix-vector multiplication, if 03045 // not empty, is necessarily zero, since R^0 = {0}. 03046 if (!result.isEmpty()) 03047 result *= beta; 03048 return; 03049 } 03050 03051 T* rdata = result.data(); 03052 T* vp = v.data(); 03053 T* mp = m.data(); 03054 int deltam = m.mod()-m.width(); 03055 03056 T* rp = rdata; 03057 // initial scaling 03058 for (int i=0;i<w;i++) 03059 *rp++ *= beta; 03060 03061 for (int j=0;j<l;j++) 03062 { 03063 T vj = *vp++; 03064 rp = rdata; 03065 for (int i=0;i<w;i++) 03066 *rp++ += alpha * vj * *mp++; 03067 mp += deltam; 03068 } 03069 } 03070 03071 /* Obsolete? Uncomment if needed 03073 // For compatibility 03074 template <class T> 03075 void transposeProductAcc(const TVec<T>& result, const TMat<T>& m, const TVec<T>& v, T alpha) 03076 { 03077 transposeProductAcc(result, m, v, alpha, 1.); 03078 } 03079 */ 03080 03081 /* Obsolete? Uncomment if needed 03082 template <class T> 03083 void compressedTransposeProductAcc(const TVec<T>& result, const TMat<T>& m, char* comprbufvec) 03084 { 03085 cout<<"using kasjdlkadja"<<endl; 03086 union { double d; char c[8]; } uni; 03087 int l=m.length(),n, idx=0; 03088 unsigned char mode; 03089 cout<<"l="<<l<<endl; 03090 for(int i=0;i<l;i++) 03091 cout<<i<<":"<<char(comprbufvec[i])<<endl; 03092 while(l>0) 03093 { 03094 read_compr_mode_and_size_ptr(comprbufvec, mode, n); 03095 if(mode==0 || mode==1) 03096 { 03097 idx+=n; 03098 cout<<"0x"<<n<<" "; 03099 l-=n; 03100 if(mode==1) 03101 { 03102 --l; 03103 result+=m(idx++); // !!!!!! 03104 cout<<"1 "; 03105 } 03106 } 03107 else if(mode==2) 03108 { 03109 while(n--) 03110 { 03111 cout<<double(*comprbufvec)<<" "<<endl; 03112 result+= m(idx++) * double(*comprbufvec++); // !!!!!! 03113 03114 --l; 03115 } 03116 } 03117 else if(mode==3) 03118 { 03119 while(n--) 03120 { 03121 memcpy(uni.c,comprbufvec,sizeof(double)); 03122 cout<<double(uni.d)<<" "<<endl; 03123 comprbufvec+=8; 03124 result+= m(idx++) * uni.d; // !!!!!!! 03125 --l; 03126 } 03127 } 03128 else 03129 PLERROR("BUG IN binread_compressed: mode is only 2 bits, so how can it be other than 0,1,2,3 ?"); 03130 } 03131 03132 if(l!=0) 03133 PLERROR("In compressed_dot_product : l is not 0 at exit of function, wrong data?"); 03134 } 03135 */ 03136 03138 template<class T> 03139 void diagonalizedFactorsProduct(TMat<T>& result, const TMat<T>& U, const TVec<T> d, const TMat<T> V, bool accumulate=false) 03140 { 03141 #ifdef BOUNDCHECK 03142 if (result.length()!=U.length() || result.width()!=V.width() || U.width()!=d.length() || V.length()!=d.length()) 03143 PLERROR("diagonalizedFactorsProduct: incompatible arguments: (%dx%d)*(%d)*(%dx%d) --> (%dx%d)", 03144 U.length(),U.width(),d.length(),V.length(),V.width(),result.length(),result.width()); 03145 #endif 03146 int n1=U.length(); 03147 int n2=U.width(); 03148 int n3=V.width(); 03149 T *r_ij = result.data(); 03150 if (accumulate) 03151 for (int i=0;i<n1;i++) 03152 { 03153 T *u_i = U[i]; 03154 for (int j=0;j<n3;j++,r_ij++) 03155 { 03156 T* d_k = d.data(); 03157 T res=0; 03158 for (int k=0;k<n2;k++,d_k++) 03159 res += *d_k * u_i[k] * V(k,j); 03160 *r_ij += res; 03161 } 03162 } 03163 else 03164 for (int i=0;i<n1;i++) 03165 { 03166 T *u_i = U[i]; 03167 for (int j=0;j<n3;j++,r_ij++) 03168 { 03169 T* d_k = d.data(); 03170 T res=0; 03171 for (int k=0;k<n2;k++,d_k++) 03172 res += *d_k * u_i[k] * V(k,j); 03173 *r_ij = res; 03174 } 03175 } 03176 } 03177 03183 template<class T> 03184 void diagonalizedFactorsProductBprop(const TMat<T>& dCdresult, const TMat<T>& U, const TVec<T> d, 03185 const TMat<T> V, TMat<T>& dCdU, TVec<T>& dCdd, TMat<T>& dCdV) 03186 { 03187 #ifdef BOUNDCHECK 03188 if (dCdU.length()!=U.length() || dCdU.width()!=U.width() || dCdd.length()!=d.length() 03189 || dCdV.length()!=V.length() || dCdV.width()!=V.width() || 03190 U.width()!=d.length() || V.length()!=d.length()) 03191 PLERROR("diagonalizedFactorsProductBprop: incompatible arguments"); 03192 #endif 03193 int n1=U.length(); 03194 int n2=U.width(); 03195 int n3=V.width(); 03196 T *dCdr_ij = dCdresult.data(); 03197 for (int i=0;i<n1;i++) 03198 { 03199 T *u_i = U[i]; 03200 T *dCdu_i = dCdU[i]; 03201 for (int j=0;j<n3;j++,dCdr_ij++) 03202 { 03203 T dcdr = *dCdr_ij; 03204 T* d_k = d.data(); 03205 T* dCdd_k = dCdd.data(); 03206 for (int k=0;k<n2;k++,d_k++,dCdd_k++) 03207 { 03208 T dk = *d_k; 03209 T u_ik = u_i[k]; 03210 T v_kj = V(k,j); 03211 dCdu_i[k] += dcdr * dk * v_kj; 03212 *dCdd_k += dcdr * u_ik * v_kj; 03213 dCdV(k,j) += dk * u_ik * dcdr; 03214 } 03215 } 03216 } 03217 } 03218 03220 template<class T> 03221 void diagonalizedFactorsProductTranspose(TMat<T>& result, const TMat<T>& U, const TVec<T> d, const TMat<T> V, bool accumulate=false) 03222 { 03223 #ifdef BOUNDCHECK 03224 if (result.length()!=U.length() || result.width()!=V.length() || U.width()!=d.length() || V.width()!=d.length()) 03225 PLERROR("diagonalizedFactorsProductTranspose: incompatible arguments: (%dx%d)*(%d)*(%dx%d)' --> (%dx%d)", 03226 U.length(),U.width(),d.length(),V.length(),V.width(),result.length(),result.width()); 03227 #endif 03228 int n1=U.length(); 03229 int n2=U.width(); 03230 int n3=V.length(); 03231 T *r_ij = result.data(); 03232 for (int i=0;i<n1;i++) 03233 { 03234 T *u_i = U[i]; 03235 for (int j=0;j<n3;j++,r_ij++) 03236 { 03237 T* d_k = d.data(); 03238 T* v_j = V[j]; 03239 T res=0; 03240 for (int k=0;k<n2;k++,d_k++) 03241 res += *d_k * u_i[k] * v_j[k]; 03242 if (accumulate) 03243 *r_ij += res; 03244 else 03245 *r_ij = res; 03246 } 03247 } 03248 } 03249 03250 // SINCE res[i,j] = sum_k U[i,k] d[k] V[j,k] ==> 03251 // gradients on dC/dU, dC/dd and dC/dV: 03252 // dC/dU[i,k] = sum_j dC/dres[i,j] d_k V[j,k] 03253 // dC/dd[k] = sum_{ij} dC/dres[i,j] U[i,k] V[j,k] 03254 // dC/dV[j,k] = sum_i dC/dres[i,j] d_k U[i,k] 03255 template<class T> 03256 void diagonalizedFactorsProductTransposeBprop(const TMat<T>& dCdresult, const TMat<T>& U, 03257 const TVec<T> d, const TMat<T> V, TMat<T>& dCdU, 03258 TVec<T>& dCdd, TMat<T>& dCdV) 03259 { 03260 #ifdef BOUNDCHECK 03261 if (dCdU.length()!=U.length() || dCdU.width()!=U.width() || dCdd.length()!=d.length() 03262 || dCdV.length()!=V.length() || dCdV.width()!=V.width() || 03263 U.width()!=d.length() || V.width()!=d.length()) 03264 PLERROR("diagonalizedFactorsProductTransposeBprop: incompatible arguments"); 03265 #endif 03266 int n1=U.length(); 03267 int n2=U.width(); 03268 int n3=V.length(); 03269 T *dCdr_ij = dCdresult.data(); 03270 for (int i=0;i<n1;i++) 03271 { 03272 T *u_i = U[i]; 03273 T *dCdu_i = dCdU[i]; 03274 for (int j=0;j<n3;j++,dCdr_ij++) 03275 { 03276 T* d_k = d.data(); 03277 T* dCdd_k = dCdd.data(); 03278 T* v_j = V[j]; 03279 T* dCdv_j = dCdV[j]; 03280 for (int k=0;k<n2;k++,d_k++,dCdd_k++) 03281 { 03282 T dcdr = *dCdr_ij; 03283 T dk = *d_k; 03284 T v_jk = v_j[k]; 03285 T u_ik = u_i[k]; 03286 dCdu_i[k] += dcdr * dk * v_jk; 03287 *dCdd_k += dcdr * u_ik * v_jk; 03288 dCdv_j[k] += dcdr * dk * u_ik; 03289 } 03290 } 03291 } 03292 } 03293 03295 template<class T> 03296 void diagonalizedFactorsTransposeProduct(TMat<T>& result, const TMat<T>& U, const TVec<T> d, const TMat<T> V, bool accumulate=false) 03297 { 03298 #ifdef BOUNDCHECK 03299 if (result.length()!=U.width() || result.width()!=V.width() || U.length()!=d.length() || V.length()!=d.length()) 03300 PLERROR("diagonalizedFactorsTransposeProduct: incompatible arguments: (%dx%d)'*(%d)*(%dx%d) --> (%dx%d)", 03301 U.length(),U.width(),d.length(),V.length(),V.width(),result.length(),result.width()); 03302 #endif 03303 int n1=U.width(); 03304 int n2=U.length(); 03305 int n3=V.width(); 03306 if (!accumulate) 03307 result.clear(); 03308 T* d_k = d.data(); 03309 for (int k=0;k<n2;k++,d_k++) 03310 { 03311 T *u_k = U[k]; 03312 T *v_k = V[k]; 03313 T *r_ij = result.data(); 03314 for (int i=0;i<n1;i++) 03315 { 03316 T u_ki = u_k[i]; 03317 for (int j=0;j<n3;j++,r_ij++) 03318 *r_ij += *d_k * u_ki * v_k[j]; 03319 } 03320 } 03321 } 03322 03323 // SINCE res[i,j] = sum_k U[k,i] d[k] V[k,j] ==> 03324 // gradients on dC/dU, dC/dd and dC/dV: 03325 // dC/dU[k,i] = d_k * sum_j dC/dres[i,j] V[k,j] 03326 // dC/dd[k] = sum_{ij} dC/dres[i,j] U[k,i] V[k,j] 03327 // dC/dV[k,j] = d_k sum_i dC/dres[i,j] U[k,i] 03328 template<class T> 03329 void diagonalizedFactorsTransposeProductBprop(const TMat<T>& dCdresult, const TMat<T>& U, 03330 const TVec<T> d, const TMat<T> V, TMat<T>& dCdU, 03331 TVec<T>& dCdd, TMat<T>& dCdV) 03332 { 03333 #ifdef BOUNDCHECK 03334 if (dCdU.length()!=U.length() || dCdU.width()!=U.width() || dCdd.length()!=d.length() 03335 || dCdV.length()!=V.length() || dCdV.width()!=V.width() || 03336 U.length()!=d.length() || V.length()!=d.length()) 03337 PLERROR("diagonalizedFactorsTransposeProductBprop: incompatible arguments"); 03338 #endif 03339 int n1=U.width(); 03340 int n2=U.length(); 03341 int n3=V.width(); 03342 T* d_k = d.data(); 03343 T* dCdd_k = dCdd.data(); 03344 for (int k=0;k<n2;k++,d_k++,dCdd_k++) 03345 { 03346 T dk = *d_k; 03347 T *u_k = U[k]; 03348 T *dCdu_k = dCdU[k]; 03349 T *v_k = V[k]; 03350 T *dCdv_k = dCdV[k]; 03351 T *dCdr_ij = dCdresult.data(); 03352 for (int i=0;i<n1;i++) 03353 { 03354 T u_ki = u_k[i]; 03355 T& dCdu_ki = dCdu_k[i]; 03356 for (int j=0;j<n3;j++,dCdr_ij++) 03357 { 03358 T dcdr = *dCdr_ij; 03359 T v_kj = v_k[j]; 03360 dCdu_ki += dcdr * dk * v_kj; 03361 *dCdd_k += dcdr * u_ki * v_kj; 03362 dCdv_k[j] += dcdr * dk * u_ki; 03363 } 03364 } 03365 } 03366 } 03367 03369 template<class T> 03370 void diagonalizedFactorsTransposeProductTranspose(TMat<T>& result, const TMat<T>& U, const TVec<T> d, const TMat<T> V, bool accumulate=false) 03371 { 03372 #ifdef BOUNDCHECK 03373 if (result.length()!=U.width() || result.width()!=V.length() || U.length()!=d.length() || V.width()!=d.length()) 03374 PLERROR("diagonalizedFactorsTransposeProductTranspose: incompatible arguments: (%dx%d)'*(%d)*(%dx%d)' --> (%dx%d)", 03375 U.length(),U.width(),d.length(),V.length(),V.width(),result.length(),result.width()); 03376 #endif 03377 int n1=U.width(); 03378 int n2=U.length(); 03379 int n3=V.length(); 03380 if (!accumulate) 03381 result.clear(); 03382 T* d_k = d.data(); 03383 for (int k=0;k<n2;k++,d_k++) 03384 { 03385 T *u_k = U[k]; 03386 T *r_ij = result.data(); 03387 for (int i=0;i<n1;i++) 03388 { 03389 T u_ki = u_k[i]; 03390 for (int j=0;j<n3;j++,r_ij++) 03391 *r_ij += *d_k * u_ki * V(j,k); 03392 } 03393 } 03394 } 03395 03396 // SINCE res[i,j] = sum_k U[k,i] d[k] V[j,k] ==> 03397 // gradients on dC/dU, dC/dd and dC/dV: 03398 // dC/dU[k,i] = d_k * sum_j dC/dres[i,j] V[j,k] 03399 // dC/dd[k] = sum_{ij} dC/dres[i,j] U[k,i] V[j,k] 03400 // dC/dV[j,k] = d_k * sum_i dC/dres[i,j] U[k,i] 03401 template<class T> 03402 void diagonalizedFactorsTransposeProductTransposeBprop(const TMat<T>& dCdresult, const TMat<T>& U, 03403 const TVec<T> d, const TMat<T> V, TMat<T>& dCdU, 03404 TVec<T>& dCdd, TMat<T>& dCdV) 03405 { 03406 #ifdef BOUNDCHECK 03407 if (dCdU.length()!=U.length() || dCdU.width()!=U.width() || dCdd.length()!=d.length() 03408 || dCdV.length()!=V.length() || dCdV.width()!=V.width() || 03409 U.length()!=d.length() || V.width()!=d.length()) 03410 PLERROR("diagonalizedFactorsTransposeProductTransposeBprop: incompatible arguments"); 03411 #endif 03412 int n1=U.width(); 03413 int n2=U.length(); 03414 int n3=V.length(); 03415 T* d_k = d.data(); 03416 T* dCdd_k = dCdd.data(); 03417 for (int k=0;k<n2;k++,d_k++,dCdd_k++) 03418 { 03419 T dk = *d_k; 03420 T *u_k = U[k]; 03421 T *dCdu_k = dCdU[k]; 03422 T *dCdr_ij = dCdresult.data(); 03423 for (int i=0;i<n1;i++) 03424 { 03425 T u_ki = u_k[i]; 03426 T& dCdu_ki = dCdu_k[i]; 03427 for (int j=0;j<n3;j++,dCdr_ij++) 03428 { 03429 T dcdr = *dCdr_ij; 03430 T v_jk = V(j,k); 03431 dCdu_ki += dcdr * dk * v_jk; 03432 *dCdd_k += dcdr * u_ki * v_jk; 03433 dCdV(j,k) += dcdr * dk * u_ki; 03434 } 03435 } 03436 } 03437 } 03438 03439 // ---------- these were previously methods of TMat --------------- 03440 03442 template<class T> 03443 T matRowDotVec(const TMat<T>& mat, int i, const TVec<T> v) 03444 { 03445 #ifdef BOUNDCHECK 03446 if (v.length()!=mat.width()) 03447 PLERROR("dotRow(%d,v), v.length_=%d != matrix width_=%d", 03448 i,v.length(),mat.width()); 03449 #endif 03450 T s = 0; 03451 T* rowi = mat.rowdata(i); 03452 T* v_=v.data(); 03453 int w=mat.width(); 03454 for (int j=0;j<w;j++) 03455 s += rowi[j] * v_[j]; 03456 return s; 03457 } 03458 03460 template<class T> 03461 T matColumnDotVec(const TMat<T>& mat, int j, const TVec<T> v) 03462 { 03463 #ifdef BOUNDCHECK 03464 if (v.length()!=mat.length()) 03465 PLERROR("dotColumn(%d,v), v.length_=%d != matrix length_=%d", 03466 j,v.length(),mat.length()); 03467 #endif 03468 T s = 0; 03469 T* colj = mat.data()+j; 03470 T* v_=v.data(); 03471 int l=mat.length(); 03472 for (int i=0;i<l;i++, colj+=mat.mod()) 03473 s += *colj * v_[i]; 03474 return s; 03475 } 03476 03478 template<class T> 03479 void matRowsDots(TVec<T> v, const TMat<T>& A, const TMat<T>& B) 03480 { 03481 #ifdef BOUNDCHECK 03482 if (A.length()!=v.length()) 03483 PLERROR("matRowDotsVec(v,A,B): v.length_=%d != A.length_=%d", 03484 v.length(),A.length()); 03485 if (A.length()!=B.length()) 03486 PLERROR("matRowDotsVec(v,A,B): A.length_=%d != B.length_=%d", 03487 A.length(),B.length()); 03488 if (A.width()!=B.width()) 03489 PLERROR("matRowDotsVec(v,A,B): A.width_=%d != B.width_=%d", 03490 A.width(),B.width()); 03491 #endif 03492 int l=A.length(), w=A.width(); 03493 T* vi = v.data(); 03494 for (int i=0;i<l;i++) 03495 { 03496 T s = 0; 03497 T* Aij = A[i]; 03498 T* Bij = B[i]; 03499 for (int j=0;j<w;j++) 03500 s += *Aij++ * *Bij++; 03501 *vi++ = s; 03502 } 03503 } 03504 03506 template<class T> 03507 void matRowsDotsAcc(TVec<T> v, const TMat<T>& A, const TMat<T>& B) 03508 { 03509 #ifdef BOUNDCHECK 03510 if (A.length()!=v.length()) 03511 PLERROR("matRowDotsVec(v,A,B): v.length_=%d != A.length_=%d", 03512 v.length(),A.length()); 03513 if (A.length()!=B.length()) 03514 PLERROR("matRowDotsVec(v,A,B): A.length_=%d != B.length_=%d", 03515 A.length(),B.length()); 03516 if (A.width()!=B.width()) 03517 PLERROR("matRowDotsVec(v,A,B): A.width_=%d != B.width_=%d", 03518 A.width(),B.width()); 03519 #endif 03520 int l=A.length(), w=A.width(); 03521 T* vi = v.data(); 03522 for (int i=0;i<l;i++) 03523 { 03524 T s = 0; 03525 T* Aij = A[i]; 03526 T* Bij = B[i]; 03527 for (int j=0;j<w;j++) 03528 s += *Aij++ * *Bij++; 03529 *vi++ += s; 03530 } 03531 } 03532 03535 template<class T> 03536 void fillItSymmetric(const TMat<T>& mat) { 03537 int m = mat.mod(); 03538 T* mat_data_to_fill; 03539 T* mat_data_to_copy; 03540 for (int i = 0; i < mat.length(); i++) { 03541 mat_data_to_fill = mat[i]; 03542 mat_data_to_copy = &mat[0][i]; 03543 for (int j = 0; j < i; j++) { 03544 *(mat_data_to_fill++) = *mat_data_to_copy; 03545 mat_data_to_copy += m; 03546 } 03547 } 03548 } 03549 03550 template<class T> 03551 void makeItSymmetric(const TMat<T>& mat, T max_dif) 03552 { 03553 if (!mat.isSquare()) 03554 PLERROR("at void makeItSymmetric, the matrix is not even square\n"); 03555 T dif; 03556 T value; 03557 bool warning_flag = false; 03558 int w=mat.width(); 03559 for (int i=0; i<mat.length()-1 ; i++) 03560 for (int j=i+1; j<w; j++) 03561 { 03562 dif = std::abs(mat[i][j] - mat[j][i]); 03563 if (dif > max_dif) 03564 { 03565 max_dif = dif; 03566 warning_flag = true; 03567 } 03568 value = (mat[i][j] + mat[j][i])/2; 03569 mat[i][j] = value; mat[j][i] = value; 03570 } 03571 if (warning_flag) 03572 PLWARNING("At void makeItSymmetric, the maximum difference %f is not affordable\n", max_dif); 03573 } 03574 03575 03576 /* DEPRECATED, use product(TVec, TMat, TVec) instead 03577 // y[i] = sum_j A[i,j] x[j] 03578 03579 template<class T> 03580 void product(const TMat<T>& mat, const TVec<T>& x, TVec<T>& y) 03581 { 03582 #ifdef BOUNDCHECK 03583 if (mat.length()!=y.length() || mat.width()!=x.length()) 03584 PLERROR("TMat(%d,%d)::product(TVec& x(%d),TVec& y(%d)), incompatible arguments", 03585 mat.length(),mat.width(),x.length(),y.length()); 03586 #endif 03587 T* x_=x.data(); 03588 T* y_=y.data(); 03589 for (int i=0;i<mat.length();i++) 03590 { 03591 T* Ai = mat[i]; 03592 T yi = 0; 03593 for (int j=0;j<mat.width();j++) 03594 yi += Ai[j] * x_[j]; 03595 y_[i]=yi; 03596 } 03597 } 03598 */ 03599 03601 template<class T> 03602 void product(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03603 { 03604 int n=m1.length(); 03605 int m=m1.width(); 03606 int l=m2.width(); 03607 #ifdef BOUNDCHECK 03608 if (n!=mat.length() || m!=m2.length() || l!=mat.width()) 03609 PLERROR("product(TMat, TMat, TMat), incompatible arguments:\n" 03610 "%dx%d <- %dx%d times %dx%d", 03611 mat.length(), mat.width(), n, m, m2.length(), l); 03612 #endif 03613 03614 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 03615 { 03616 // Size zero: no need to bother computing anything. 03617 // In such a case, the result of the matrix-matrix multiplication, if 03618 // not empty, is necessarily zero, since R^0 = {0}. 03619 if (!mat.isEmpty()) 03620 mat.clear(); 03621 return; 03622 } 03623 03624 for (int i=0;i<n;i++) 03625 { 03626 const T* m1i = m1[i]; 03627 T* mi = mat[i]; 03628 for (int j=0;j<l;j++) 03629 { 03630 T s=0; 03631 const T* m2kj = m2.data()+j; 03632 for (int k=0;k<m;k++,m2kj+=m2.mod()) 03633 s += m1i[k] * (*m2kj); 03634 mi[j] = s; 03635 } 03636 } 03637 } 03638 03640 template<class T> 03641 void productAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03642 { 03643 int n=m1.length(); 03644 int m=m1.width(); 03645 int l=m2.width(); 03646 #ifdef BOUNDCHECK 03647 if (n!=mat.length() || m!=m2.length() || l!=mat.width()) 03648 PLERROR("productAcc(TMat, TMat, TMat), incompatible arguments:\n" 03649 "%dx%d <- %dx%d times %dx%d", 03650 mat.length(), mat.width(), n, m, m2.length(), l); 03651 #endif 03652 03653 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 03654 { 03655 // Size zero: no need to bother computing anything. 03656 // In such a case, the result of the matrix-matrix multiplication, if 03657 // not empty, is necessarily zero, since R^0 = {0}. 03658 return; 03659 } 03660 03661 for (int i=0;i<n;i++) 03662 { 03663 const T* m1i = m1[i]; 03664 T* mi = mat[i]; 03665 for (int j=0;j<l;j++) 03666 { 03667 T s=0; 03668 T* m2kj = m2.data()+j; 03669 for (int k=0;k<m;k++,m2kj+=m2.mod()) 03670 s += m1i[k] * (*m2kj); 03671 mi[j] += s; 03672 } 03673 } 03674 } 03675 03677 // (Will use the transpose of m1 and/or m2 instead, 03678 // if you set the corresponding flags to true) 03679 template<class T> 03680 void productScaleAcc(const TMat<T>& mat, 03681 const TMat<T>& m1, bool transpose_m1, 03682 const TMat<T>& m2, bool transpose_m2, 03683 T alpha, T beta) 03684 { 03685 // Boundary checking is done in called functions 03686 if (transpose_m1) 03687 if (transpose_m2) // transpose_m1 && transpose_m2 03688 transposeTransposeProductScaleAcc(mat, m1, m2, alpha, beta); 03689 else // transpose_m1 && !transpose_m2 03690 transposeProductScaleAcc(mat, m1, m2, alpha, beta); 03691 else 03692 if (transpose_m2) // !transpose_m1 && transpose_m2 03693 productTransposeScaleAcc(mat, m1, m2, alpha, beta); 03694 else // !transpose_m1 && !transpose_m2 03695 productScaleAcc(mat, m1, m2, alpha, beta); 03696 } 03697 03699 template<class T> 03700 void productScaleAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2, 03701 T alpha, T beta) 03702 { 03703 int n=m1.length(); 03704 int m=m1.width(); 03705 int l=m2.width(); 03706 #ifdef BOUNDCHECK 03707 if (n!=mat.length() || m!=m2.length() || l!=mat.width()) 03708 PLERROR("productScaleAcc(TMat, TMat, TMat), incompatible arguments:\n" 03709 "%dx%d <- %dx%d times %dx%d", 03710 mat.length(), mat.width(), n, m, m2.length(), l); 03711 #endif 03712 03713 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 03714 { 03715 // Size zero: no need to bother computing anything. 03716 // In such a case, the result of the matrix-matrix multiplication, if 03717 // not empty, is necessarily zero, since R^0 = {0}. 03718 if (!mat.isEmpty()) 03719 mat *= beta; 03720 return; 03721 } 03722 03723 for (int i=0;i<n;i++) 03724 { 03725 const T* m1i = m1[i]; 03726 T* mi = mat[i]; 03727 for (int j=0;j<l;j++) 03728 { 03729 T s=0; 03730 T* m2kj = m2.data()+j; 03731 for (int k=0;k<m;k++,m2kj+=m2.mod()) 03732 s += m1i[k] * (*m2kj); 03733 mi[j] = alpha * s + beta * mi[j]; 03734 } 03735 } 03736 } 03737 03738 // result[i,j] += sum_k m1[i,k] * m2[k,j]^2 03739 template<class T> 03740 void product2Acc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03741 { 03742 #ifdef BOUNDCHECK 03743 if (m1.width()!=m2.length() || mat.length()!=m1.length() || mat.width()!=m2.width()) 03744 PLERROR("product2Acc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d", 03745 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03746 #endif 03747 int n=m1.length(); 03748 int m=m1.width(); 03749 int l=m2.width(); 03750 for (int i=0;i<n;i++) 03751 { 03752 const T* m1i = m1[i]; 03753 T* mi = mat[i]; 03754 for (int j=0;j<l;j++) 03755 { 03756 T s=0; 03757 T* m2kj = m2.data()+j; 03758 for (int k=0;k<m;k++,m2kj+=m2.mod()) 03759 s += m1i[k] * (*m2kj) * (*m2kj); 03760 mi[j] += s; 03761 } 03762 } 03763 } 03764 03765 // result[i,j] += sum_k m1[i,k]^2 * m2[k,j] 03766 template<class T> 03767 void squareProductAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03768 { 03769 #ifdef BOUNDCHECK 03770 if (m1.width()!=m2.length() || mat.length()!=m1.length() || mat.width()!=m2.width()) 03771 PLERROR("squareProductAcc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d", 03772 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03773 #endif 03774 int n=m1.length(); 03775 int m=m1.width(); 03776 int l=m2.width(); 03777 for (int i=0;i<n;i++) 03778 { 03779 const T* m1i = m1[i]; 03780 T* mi = mat[i]; 03781 for (int j=0;j<l;j++) 03782 { 03783 T s=0; 03784 T* m2kj = m2.data()+j; 03785 for (int k=0;k<m;k++,m2kj+=m2.mod()) 03786 { 03787 T m1ik=m1i[k]; 03788 s += m1ik*m1ik * (*m2kj); 03789 } 03790 mi[j] += s; 03791 } 03792 } 03793 } 03794 03795 // result[i][j] = v1[i] * v2[j] 03796 03797 template<class T> 03798 void externalProduct(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2) 03799 { 03800 #ifdef BOUNDCHECK 03801 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03802 PLERROR("externalProduct(Vec,Vec), incompatible arguments %dx%d= %d times %d", 03803 mat.length(),mat.width(),v1.length(), v2.length()); 03804 #endif 03805 const T* v_1=v1.data(); 03806 const T* v_2=v2.data(); 03807 int w=mat.width(); 03808 for (int i=0;i<mat.length();i++) 03809 { 03810 T* mi = mat[i]; 03811 T v1i = v_1[i]; 03812 for (int j=0;j<w;j++) 03813 mi[j] = v1i * v_2[j]; 03814 } 03815 } 03816 03817 // mat[i][j] += v1[i] * v2[j] 03818 template<class T> 03819 void externalProductAcc(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2) 03820 { 03821 #ifdef BOUNDCHECK 03822 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03823 PLERROR("externalProductAcc(Vec,Vec), incompatible arguments %dx%d= %d times %d", 03824 mat.length(),mat.width(),v1.length(), v2.length()); 03825 #endif 03826 03827 T* v_1=v1.data(); 03828 T* v_2=v2.data(); 03829 T* mp = mat.data(); 03830 int l = mat.length(); 03831 int w = mat.width(); 03832 03833 if(mat.isCompact()) 03834 { 03835 T* pv1 = v_1; 03836 for(int i=0; i<l; i++) 03837 { 03838 T* pv2 = v_2; 03839 T val = *pv1++; 03840 for(int j=0; j<w; j++) 03841 *mp++ += val * *pv2++; 03842 } 03843 } 03844 else 03845 { 03846 cerr << "!"; 03847 for (int i=0;i<l;i++) 03848 { 03849 T* mi = mat[i]; 03850 T v1i = v_1[i]; 03851 for (int j=0;j<w;j++) 03852 mi[j] += v1i * v_2[j]; 03853 } 03854 } 03855 } 03856 03857 // mat[i][j] += gamma * v1[i] * v2[j] 03858 template<class T> 03859 void externalProductScaleAcc(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2, T gamma) 03860 { 03861 Profiler::pl_profile_start("externalProductScaleAcc T"); 03862 03863 #ifdef BOUNDCHECK 03864 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03865 PLERROR("externalProductScaleAcc(Vec,Vec), incompatible arguments %dx%d= %d times %d", 03866 mat.length(),mat.width(),v1.length(), v2.length()); 03867 #endif 03868 const T* v_1=v1.data(); 03869 const T* v_2=v2.data(); 03870 int w=mat.width(); 03871 for (int i=0;i<mat.length();i++) 03872 { 03873 T* mi = mat[i]; 03874 T v1i = v_1[i]; 03875 for (int j=0;j<w;j++) 03876 mi[j] += gamma * v1i * v_2[j]; 03877 } 03878 Profiler::pl_profile_end("externalProductScaleAcc T"); 03879 } 03880 03881 // mat[i][j] = alpha * mat[i][j] + gamma * v1[i] * v2[j] 03882 template<class T> 03883 void externalProductScaleAcc(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2, T gamma, T alpha) 03884 { 03885 Profiler::pl_profile_start("externalProductScaleAcc T"); 03886 03887 #ifdef BOUNDCHECK 03888 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03889 PLERROR("externalProductScaleAcc(Vec,Vec), incompatible arguments %dx%d= %d times %d", 03890 mat.length(),mat.width(),v1.length(), v2.length()); 03891 #endif 03892 const T* v_1=v1.data(); 03893 const T* v_2=v2.data(); 03894 int w=mat.width(); 03895 for (int i=0;i<mat.length();i++) 03896 { 03897 T* mi = mat[i]; 03898 T v1i = v_1[i]; 03899 for (int j=0;j<w;j++) 03900 mi[j] = alpha*mi[j] + gamma * v1i * v_2[j]; 03901 } 03902 Profiler::pl_profile_end("externalProductScaleAcc T"); 03903 } 03904 03905 // mat[i][j] *= v1[i] * v2[j] 03906 template<class T> 03907 void externalProductMultUpdate(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2) 03908 { 03909 #ifdef BOUNDCHECK 03910 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03911 PLERROR("externalProductMultUpdate(mat,v1,v2), incompatible arguments %dx%d= %d times %d", 03912 mat.length(),mat.width(),v1.length(), v2.length()); 03913 #endif 03914 const T* v_1=v1.data(); 03915 const T* v_2=v2.data(); 03916 const int N = mat.length(); 03917 const int M = mat.width(); 03918 for (int i=0 ; i<N ; ++i) { 03919 T* mi = mat[i]; 03920 T v1i = v_1[i]; 03921 for (int j=0; j<M ; ++j) 03922 mi[j] *= v1i * v_2[j]; 03923 } 03924 } 03925 03926 03927 // mat[i][j] /= v1[i] * v2[j] 03928 template<class T> 03929 void externalProductDivUpdate(const TMat<T>& mat, const TVec<T>& v1, const TVec<T>& v2) 03930 { 03931 #ifdef BOUNDCHECK 03932 if (v1.length()!=mat.length() || mat.width()!=v2.length()) 03933 PLERROR("externalProductDivUpdate(mat,v1,v2), incompatible arguments %dx%d= %d times %d", 03934 mat.length(),mat.width(),v1.length(), v2.length()); 03935 #endif 03936 const T* v_1=v1.data(); 03937 const T* v_2=v2.data(); 03938 const int N = mat.length(); 03939 const int M = mat.width(); 03940 for (int i=0 ; i<N ; ++i) { 03941 T* mi = mat[i]; 03942 T v1i = v_1[i]; 03943 for (int j=0; j<M ; ++j) 03944 mi[j] /= v1i * v_2[j]; 03945 } 03946 } 03947 03948 03950 template<class T> 03951 void productTranspose(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03952 { 03953 int n=m1.length(); 03954 int m=m1.width(); 03955 int l=m2.length(); 03956 #ifdef BOUNDCHECK 03957 if (n!=mat.length() || m!=m2.width() || l!=mat.width()) 03958 PLERROR("productTranspose(TMat, TMat, TMat), incompatible arguments:\n" 03959 "%dx%d <- %dx%d times %dx%d'", 03960 mat.length(), mat.width(), n, m, l, m2.width()); 03961 #endif 03962 03963 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 03964 { 03965 // Size zero: no need to bother computing anything. 03966 // In such a case, the result of the matrix-matrix multiplication, if 03967 // not empty, is necessarily zero, since R^0 = {0}. 03968 if (!mat.isEmpty()) 03969 mat.clear(); 03970 return; 03971 } 03972 03973 for (int i=0;i<n;i++) 03974 { 03975 const T* m1i = m1[i]; 03976 T* mi = mat[i]; 03977 for (int j=0;j<l;j++) 03978 { 03979 T s=0; 03980 const T* m2j = m2[j]; 03981 for (int k=0;k<m;k++) 03982 s += m1i[k] * m2j[k]; 03983 mi[j] = s; 03984 } 03985 } 03986 } 03987 03988 // result[i,j] = sum_k m1[i,k]^2 * m2[j,k] 03989 template<class T> 03990 void squareProductTranspose(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 03991 { 03992 #ifdef BOUNDCHECK 03993 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 03994 PLERROR("squareProductTranspose(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 03995 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 03996 #endif 03997 int n=m1.length(); 03998 int m=m1.width(); 03999 int l=m2.length(); 04000 for (int i=0;i<n;i++) 04001 { 04002 const T* m1i = m1[i]; 04003 T* mi = mat[i]; 04004 for (int j=0;j<l;j++) 04005 { 04006 T s=0; 04007 const T* m2j = m2[j]; 04008 for (int k=0;k<m;k++) 04009 { 04010 T m1ik=m1i[k]; 04011 s += m1ik*m1ik * m2j[k]; 04012 } 04013 mi[j] = s; 04014 } 04015 } 04016 } 04017 04018 // result[i,j] = sum_k m1[i,k] * m2[j,k]^2 04019 template<class T> 04020 void product2Transpose(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 04021 { 04022 #ifdef BOUNDCHECK 04023 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 04024 PLERROR("product2Transpose(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 04025 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 04026 #endif 04027 int n=m1.length(); 04028 int m=m1.width(); 04029 int l=m2.length(); 04030 for (int i=0;i<n;i++) 04031 { 04032 const T* m1i = m1[i]; 04033 T* mi = mat[i]; 04034 for (int j=0;j<l;j++) 04035 { 04036 T s=0; 04037 const T* m2j = m2[j]; 04038 for (int k=0;k<m;k++) 04039 { 04040 T m2jk=m2j[k]; 04041 s += m1i[k] * m2jk*m2jk; 04042 } 04043 mi[j] = s; 04044 } 04045 } 04046 } 04047 04049 template<class T> 04050 void productTransposeAcc(const TMat<T>& mat, const TMat<T>& m1, 04051 const TMat<T>& m2) 04052 { 04053 int n=m1.length(); 04054 int m=m1.width(); 04055 int l=m2.length(); 04056 #ifdef BOUNDCHECK 04057 if (n!=mat.length() || m!=m2.width() || l!=mat.width()) 04058 PLERROR("productTransposeAcc(TMat, TMat, TMat), incompatible arguments" 04059 ":\n" 04060 "%dx%d <- %dx%d times %dx%d'", 04061 mat.length(), mat.width(), n, m, l, m2.width()); 04062 #endif 04063 04064 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 04065 { 04066 // Size zero: no need to bother computing anything. 04067 // In such a case, the result of the matrix-matrix multiplication, if 04068 // not empty, is necessarily zero, since R^0 = {0}. 04069 return; 04070 } 04071 04072 for (int i=0;i<n;i++) 04073 { 04074 const T* m1i = m1[i]; 04075 T* mi = mat[i]; 04076 for (int j=0;j<l;j++) 04077 { 04078 T s=0; 04079 const T* m2j = m2[j]; 04080 for (int k=0;k<m;k++) 04081 s += m1i[k] * m2j[k]; 04082 mi[j] += s; 04083 } 04084 } 04085 } 04086 04088 template<class T> 04089 void productTransposeScaleAcc(const TMat<T>& mat, const TMat<T>& m1, 04090 const TMat<T>& m2, T alpha, T beta) 04091 { 04092 int n=m1.length(); 04093 int m=m1.width(); 04094 int l=m2.length(); 04095 #ifdef BOUNDCHECK 04096 if (n!=mat.length() || m!=m2.width() || l!=mat.width()) 04097 PLERROR("productTransposeScaleAcc(TMat, TMat, TMat), incompatible" 04098 " arguments:\n" 04099 "%dx%d <- %dx%d times %dx%d'", 04100 mat.length(), mat.width(), n, m, l, m2.width()); 04101 #endif 04102 04103 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 04104 { 04105 // Size zero: no need to bother computing anything. 04106 // In such a case, the result of the matrix-matrix multiplication, if 04107 // not empty, is necessarily zero, since R^0 = {0}. 04108 if (!mat.isEmpty()) 04109 mat *= beta; 04110 return; 04111 } 04112 04113 for (int i=0;i<n;i++) 04114 { 04115 const T* m1i = m1[i]; 04116 T* mi = mat[i]; 04117 for (int j=0;j<l;j++) 04118 { 04119 T s=0; 04120 const T* m2j = m2[j]; 04121 for (int k=0;k<m;k++) 04122 s += m1i[k] * m2j[k]; 04123 mi[j] = alpha * s + beta * mi[j]; 04124 } 04125 } 04126 } 04127 04128 // result[i,j] += sum_k m1[i,k] * m2[j,k]^2 04129 template<class T> 04130 void product2TransposeAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 04131 { 04132 #ifdef BOUNDCHECK 04133 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 04134 PLERROR("product2TransposeAcc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 04135 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 04136 #endif 04137 int n=m1.length(); 04138 int m=m1.width(); 04139 int l=m2.length(); 04140 for (int i=0;i<n;i++) 04141 { 04142 const T* m1i = m1[i]; 04143 T* mi = mat[i]; 04144 for (int j=0;j<l;j++) 04145 { 04146 T s=0; 04147 const T* m2j = m2[j]; 04148 for (int k=0;k<m;k++) 04149 { 04150 T m2jk=m2j[k]; 04151 s += m1i[k] * m2jk*m2jk; 04152 } 04153 mi[j] += s; 04154 } 04155 } 04156 } 04157 04158 // result[i,j] += sum_k m1[i,k]^2 * m2[j,k] 04159 template<class T> 04160 void squareProductTransposeAcc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 04161 { 04162 #ifdef BOUNDCHECK 04163 if (m1.width()!=m2.width() || mat.length()!=m1.length() || mat.width()!=m2.length()) 04164 PLERROR("squareProductTransposeAcc(Mat,Mat), incompatible arguments %dx%d= %dx%d times %dx%d'", 04165 mat.length(),mat.width(),m1.length(),m1.width(), m2.length(),m2.width()); 04166 #endif 04167 int n=m1.length(); 04168 int m=m1.width(); 04169 int l=m2.length(); 04170 for (int i=0;i<n;i++) 04171 { 04172 const T* m1i = m1[i]; 04173 T* mi = mat[i]; 04174 for (int j=0;j<l;j++) 04175 { 04176 T s=0; 04177 const T* m2j = m2[j]; 04178 for (int k=0;k<m;k++) 04179 { 04180 T m1ik=m1i[k]; 04181 s += m1ik*m1ik * m2j[k]; 04182 } 04183 mi[j] += s; 04184 } 04185 } 04186 } 04187 04189 template<class T> 04190 void transposeProduct(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 04191 { 04192 int n=m1.width(); 04193 int m=m1.length(); 04194 int l=m2.width(); 04195 #ifdef BOUNDCHECK 04196 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 04197 PLERROR("transposeProduct(TMat, TMat, TMat), incompatible arguments:\n" 04198 "%dx%d <- %dx%d' times %dx%d", 04199 mat.length(), mat.width(), m, n, m2.length(), l); 04200 #endif 04201 04202 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 04203 { 04204 // Size zero: no need to bother computing anything. 04205 // In such a case, the result of the matrix-matrix multiplication, if 04206 // not empty, is necessarily zero, since R^0 = {0}. 04207 if (!mat.isEmpty()) 04208 mat.clear(); 04209 return; 04210 } 04211 04212 mat.clear(); 04213 for (int i=0;i<n;i++) 04214 { 04215 T* m1ki = m1.data()+i; 04216 T* mi = mat[i]; 04217 for (int k=0;k<m;k++,m1ki+=m1.mod()) 04218 { 04219 const T* m2k = m2[k]; 04220 T m1_ki = *m1ki; 04221 for (int j=0;j<l;j++) 04222 mi[j] += m1_ki * m2k[j]; 04223 } 04224 } 04225 } 04226 04227 // result[i,j] = sum_k m1[k,i] * m2[k,j]^2 04228 template<class T> 04229 void transposeProduct2(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 04230 { 04231 int n=m1.width(); 04232 int m=m1.length(); 04233 int l=m2.width(); 04234 #ifdef BOUNDCHECK 04235 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 04236 PLERROR("transposeProduct2(Mat,Mat), incompatible arguments " 04237 "%dx%d' times %dx%d into %dx%d", 04238 m1.length(),m1.width(), m2.length(),m2.width(), mat.length(), mat.width()); 04239 #endif 04240 mat.clear(); 04241 for (int i=0;i<n;i++) 04242 { 04243 T* m1ki = m1.data()+i; 04244 T* mi = mat[i]; 04245 for (int k=0;k<m;k++,m1ki+=m1.mod()) 04246 { 04247 const T* m2k = m2[k]; 04248 T m1_ki = *m1ki; 04249 for (int j=0;j<l;j++) 04250 { 04251 T m2kj=m2k[j]; 04252 mi[j] += m1_ki * m2kj*m2kj; 04253 } 04254 } 04255 } 04256 } 04257 04259 template<class T> 04260 void transposeProductAcc(const TMat<T>& mat, const TMat<T>& m1, 04261 const TMat<T>& m2) 04262 { 04263 int n=m1.width(); 04264 int m=m1.length(); 04265 int l=m2.width(); 04266 #ifdef BOUNDCHECK 04267 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 04268 PLERROR("transposeProductAcc(TMat, TMat, TMat), incompatible" 04269 " arguments:\n" 04270 "%dx%d <- %dx%d' times %dx%d", 04271 mat.length(), mat.width(), m, n, m2.length(), l); 04272 #endif 04273 04274 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 04275 { 04276 // Size zero: no need to bother computing anything. 04277 // In such a case, the result of the matrix-matrix multiplication, if 04278 // not empty, is necessarily zero, since R^0 = {0}. 04279 return; 04280 } 04281 04282 for (int i=0;i<n;i++) 04283 { 04284 T* m1ki = m1.data()+i; 04285 T* mi = mat[i]; 04286 for (int k=0;k<m;k++,m1ki+=m1.mod()) 04287 { 04288 const T* m2k = m2[k]; 04289 T m1_ki = *m1ki; 04290 for (int j=0;j<l;j++) 04291 mi[j] += m1_ki * m2k[j]; 04292 } 04293 } 04294 } 04295 04296 // mat[i,j] = alpha * sum_k m1[k,i] * m2[k,j] + beta * mat[i,j] 04297 template<class T> 04298 void transposeProductScaleAcc(const TMat<T>& mat, const TMat<T>& m1, 04299 const TMat<T>& m2, T alpha, T beta) 04300 { 04301 int n=m1.width(); 04302 int m=m1.length(); 04303 int l=m2.width(); 04304 #ifdef BOUNDCHECK 04305 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 04306 PLERROR("transposeProductScaleAcc(TMat, TMat, TMat), incompatible" 04307 " arguments:\n" 04308 "%dx%d <- %dx%d' times %dx%d", 04309 mat.length(), mat.width(), m, n, m2.length(), l); 04310 #endif 04311 04312 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 04313 { 04314 // Size zero: no need to bother computing anything. 04315 // In such a case, the result of the matrix-matrix multiplication, if 04316 // not empty, is necessarily zero, since R^0 = {0}. 04317 if (!mat.isEmpty()) 04318 mat *= beta; 04319 return; 04320 } 04321 04322 for (int i=0;i<n;i++) 04323 { 04324 T* m1ki = m1.data()+i; 04325 T* mi = mat[i]; 04326 04327 // initial scaling 04328 for (int j=0;j<l;j++) 04329 mi[j] *= beta; 04330 04331 for (int k=0;k<m;k++,m1ki+=m1.mod()) 04332 { 04333 const T* m2k = m2[k]; 04334 T m1_ki = *m1ki; 04335 for (int j=0;j<l;j++) 04336 mi[j] += alpha * m1_ki * m2k[j]; 04337 } 04338 } 04339 } 04340 04341 // result[i,j] += sum_k m1[k,i] * m2[k,j]^2 04342 template<class T> 04343 void transposeProduct2Acc(const TMat<T>& mat, const TMat<T>& m1, const TMat<T>& m2) 04344 { 04345 int n=m1.width(); 04346 int m=m1.length(); 04347 int l=m2.width(); 04348 #ifdef BOUNDCHECK 04349 if (m!=m2.length() || mat.length()!=n || mat.width()!=l) 04350 PLERROR("transposeProduct2Acc(Mat,Mat), incompatible arguments " 04351 "%dx%d' times %dx%d into %dx%d", 04352 m1.length(),m1.width(), m2.length(),m2.width(), mat.length(), mat.width()); 04353 #endif 04354 for (int i=0;i<n;i++) 04355 { 04356 T* m1ki = m1.data()+i; 04357 T* mi = mat[i]; 04358 for (int k=0;k<m;k++,m1ki+=m1.mod()) 04359 { 04360 const T* m2k = m2[k]; 04361 T m1_ki = *m1ki; 04362 for (int j=0;j<l;j++) 04363 { 04364 T m2kj = m2k[j]; 04365 mi[j] += m1_ki * m2kj * m2kj; 04366 } 04367 } 04368 } 04369 } 04370 04372 template<class T> 04373 void transposeTransposeProduct(const TMat<T>& mat, const TMat<T>& m1, 04374 const TMat<T>& m2) 04375 { 04376 int n=m1.width(); 04377 int m=m1.length(); 04378 int l=m2.length(); 04379 #ifdef BOUNDCHECK 04380 if (n!=mat.length() || m!=m2.width() || l!=mat.width()) 04381 PLERROR("transposeTransposeProduct(TMat, TMat, TMat), incompatible" 04382 " arguments:\n" 04383 "%dx%d <- %dx%d' times %dx%d'", 04384 mat.length(), mat.width(), m, n, l, m2.width()); 04385 #endif 04386 04387 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 04388 { 04389 // Size zero: no need to bother computing anything. 04390 // In such a case, the result of the matrix-matrix multiplication, if 04391 // not empty, is necessarily zero, since R^0 = {0}. 04392 if (!mat.isEmpty()) 04393 mat.clear(); 04394 return; 04395 } 04396 04397 for (int i=0;i<n;i++) 04398 { 04399 T* m1ki0 = m1.data()+i; 04400 T* mi = mat[i]; 04401 for (int j=0;j<l;j++) 04402 { 04403 T s=0; 04404 const T* m2j = m2[j]; 04405 T* m1ki = m1ki0; 04406 for (int k=0;k<m;k++,m1ki+=m1.mod()) 04407 s += (*m1ki) * m2j[k]; 04408 mi[j] = s; 04409 } 04410 } 04411 } 04412 04414 template<class T> 04415 void transposeTransposeProductAcc(const TMat<T>& mat, const TMat<T>& m1, 04416 const TMat<T>& m2) 04417 { 04418 int n=m1.width(); 04419 int m=m1.length(); 04420 int l=m2.length(); 04421 #ifdef BOUNDCHECK 04422 if (n!=mat.length() || m!=m2.width() || l!=mat.width()) 04423 PLERROR("transposeTransposeProductAcc(TMat, TMat, TMat), incompatible" 04424 " arguments:\n" 04425 "%dx%d <- %dx%d' times %dx%d'", 04426 mat.length(), mat.width(), m, n, l, m2.width()); 04427 #endif 04428 04429 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 04430 { 04431 // Size zero: no need to bother computing anything. 04432 // In such a case, the result of the matrix-matrix multiplication, if 04433 // not empty, is necessarily zero, since R^0 = {0}. 04434 return; 04435 } 04436 04437 for (int i=0;i<n;i++) 04438 { 04439 T* m1ki0 = m1.data()+i; 04440 T* mi = mat[i]; 04441 for (int j=0;j<l;j++) 04442 { 04443 T s=0; 04444 const T* m2j = m2[j]; 04445 T* m1ki = m1ki0; 04446 for (int k=0;k<m;k++,m1ki+=m1.mod()) 04447 s += (*m1ki) * m2j[k]; 04448 mi[j] += s; 04449 } 04450 } 04451 } 04452 04454 template<class T> 04455 void transposeTransposeProductScaleAcc(const TMat<T>& mat, const TMat<T>& m1, 04456 const TMat<T>& m2, T alpha, T beta) 04457 { 04458 int n=m1.width(); 04459 int m=m1.length(); 04460 int l=m2.length(); 04461 #ifdef BOUNDCHECK 04462 if (n!=mat.length() || m!=m2.width() || l!=mat.width()) 04463 PLERROR("transposeTransposeProductScaleAcc(TMat, TMat, TMat)," 04464 " incompatible arguments:\n" 04465 "%dx%d <- %dx%d' times %dx%d'", 04466 mat.length(), mat.width(), m, n, l, m2.width()); 04467 #endif 04468 04469 if (m1.isEmpty() || m2.isEmpty() || mat.isEmpty()) 04470 { 04471 // Size zero: no need to bother computing anything. 04472 // In such a case, the result of the matrix-matrix multiplication, if 04473 // not empty, is necessarily zero, since R^0 = {0}. 04474 if (!mat.isEmpty()) 04475 mat *= beta; 04476 return; 04477 } 04478 04479 for (int i=0;i<n;i++) 04480 { 04481 T* m1ki0 = m1.data()+i; 04482 T* mi = mat[i]; 04483 for (int j=0;j<l;j++) 04484 { 04485 T s=0; 04486 const T* m2j = m2[j]; 04487 T* m1ki = m1ki0; 04488 for (int k=0;k<m;k++,m1ki+=m1.mod()) 04489 s += (*m1ki) * m2j[k]; 04490 mi[j] = alpha * s + beta * mi[j]; 04491 } 04492 } 04493 } 04494 04495 template<class T> 04496 T trace(const TMat<T>& mat) 04497 { 04498 if (!mat.isSquare()) 04499 PLERROR( "In trace()\nThe matrix must be square." ); 04500 T tr = mat.firstElement(); 04501 for ( int i = 1; i < mat.length(); i++ ) 04502 tr += mat(i,i); 04503 return tr; 04504 } 04505 04507 template<class T> 04508 void regularizeMatrix(const TMat<T>& mat, T tolerance) 04509 { 04510 T reg; 04511 T* k; 04512 int shift; 04513 reg = tolerance * trace(mat); 04514 k = mat.data(); 04515 shift = mat.mod() + 1; 04516 for (int i = 0; i < mat.length(); i++) { 04517 *k += reg; 04518 k += shift; 04519 } 04520 } 04521 04522 04523 template<class T> 04524 void makeRowsSumTo1(const TMat<T>& mat) 04525 { 04526 for (int i = 0; i < mat.length(); i++) 04527 { 04528 TVec<T> row_i = mat(i); 04529 divide(row_i, sum(row_i), row_i); 04530 } 04531 } 04532 04533 // result[i,j] = x[i,j]*scale; 04534 template<class T> 04535 void multiply(const TMat<T>& result, const TMat<T>& x, T scale) 04536 { 04537 #ifdef BOUNDCHECK 04538 if (result.length()!=x.length() || result.width()!=x.width()) 04539 PLERROR("multiply incompatible dimensions: %dx%d <- %dx%d", 04540 result.length(),result.width(),x.length(),x.width()); 04541 #endif 04542 if(result.isCompact() && x.isCompact()) 04543 { 04544 typename TMat<T>::compact_iterator itm = result.compact_begin(); 04545 typename TMat<T>::compact_iterator itmend = result.compact_end(); 04546 typename TMat<T>::compact_iterator itx = x.compact_begin(); 04547 for(; itm!=itmend; ++itm, ++itx) 04548 *itm = *itx * scale; 04549 } 04550 else // use non-compact iterators 04551 { 04552 typename TMat<T>::iterator itm = result.begin(); 04553 typename TMat<T>::iterator itmend = result.end(); 04554 typename TMat<T>::iterator itx = x.begin(); 04555 for(; itm!=itmend; ++itm, ++itx) 04556 *itm = *itx * scale; 04557 } 04558 } 04559 04560 // result[i,j] = x[i,j]*y[i] or x[i,j]*y[j] (transpose case) 04561 template<class T> 04562 void multiply(TMat<T>& result, const TMat<T>& x, const TVec<T>& y, bool transpose=false) 04563 { 04564 PLASSERT_MSG(transpose && x.width()==y.length() || 04565 !transpose && x.length()==y.length(), 04566 "multiply matrix rows or columns by vector: incompatible dimensions"); 04567 result.resize(x.length(),x.width()); 04568 int w=x.width(); 04569 if(result.isCompact() && x.isCompact()) 04570 { 04571 typename TMat<T>::compact_iterator itm = result.compact_begin(); 04572 typename TMat<T>::compact_iterator itx = x.compact_begin(); 04573 typename TVec<T>::iterator ity = y.begin(); 04574 if (transpose) 04575 for (int i=0;i<x.length();i++) 04576 { 04577 ity = y.begin(); 04578 for (int j=0;j<w;j++,++itx,++itm,++ity) 04579 *itm = *itx * *ity; 04580 } 04581 else 04582 for (int i=0;i<x.length();i++,++ity) 04583 for (int j=0;j<w;j++,++itx,++itm) 04584 *itm = *itx * *ity; 04585 } 04586 else // use non-compact iterators 04587 { 04588 typename TMat<T>::iterator itm = result.begin(); 04589 typename TMat<T>::iterator itx = x.begin(); 04590 typename TVec<T>::iterator ity = y.begin(); 04591 if (transpose) 04592 for (int i=0;i<x.length();i++) 04593 { 04594 ity = y.begin(); 04595 for (int j=0;j<w;j++,++itx,++itm,++ity) 04596 *itm = *itx * *ity; 04597 } 04598 else 04599 for (int i=0;i<x.length();i++,++ity) 04600 for (int j=0;j<w;j++,++itx,++itm) 04601 *itm = *itx * *ity; 04602 } 04603 } 04604 04605 template<class T> 04606 inline TMat<T> operator*(const TMat<T>& m, const T& scalar) 04607 { 04608 TMat<T> result(m.length(),m.width()); 04609 multiply(result, m, scalar); 04610 return result; 04611 } 04612 04613 template<class T> 04614 inline TMat<T> operator*(const T& scalar, const TMat<T>& m) 04615 { return m * scalar;} 04616 04617 // Will not work properly for integers... 04618 template<class T> 04619 inline TMat<T> operator/(const TMat<T>& m, const T& scalar) 04620 { return m * (T(1)/scalar); } 04621 04622 // result[i,j] += x[i,j]*scale; 04623 template<class T> 04624 void multiplyAcc(const TMat<T>& mat, const TMat<T>& x, T scale) 04625 { 04626 #ifdef BOUNDCHECK 04627 if (mat.length()!=x.length() || mat.width()!=x.width()) 04628 PLERROR("multiplyAcc incompatible dimensions: %dx%d <- %dx%d", 04629 mat.length(),mat.width(),x.length(),x.width()); 04630 #endif 04631 if(mat.isCompact() && x.isCompact()) 04632 { 04633 typename TMat<T>::compact_iterator itm = mat.compact_begin(); 04634 typename TMat<T>::compact_iterator itmend = mat.compact_end(); 04635 typename TMat<T>::compact_iterator itx = x.compact_begin(); 04636 for(; itm!=itmend; ++itm, ++itx) 04637 *itm += *itx * scale; 04638 } 04639 else // use non-compact iterators 04640 { 04641 typename TMat<T>::iterator itm = mat.begin(); 04642 typename TMat<T>::iterator itmend = mat.end(); 04643 typename TMat<T>::iterator itx = x.begin(); 04644 for(; itm!=itmend; ++itm, ++itx) 04645 *itm += *itx * scale; 04646 } 04647 } 04648 04649 // result[i,j] += x[i,j]*y[i,j]; 04650 template<class T> 04651 void multiplyAcc(const TMat<T>& mat, const TMat<T>& x, const TMat<T>& y) 04652 { 04653 int n=mat.length()*mat.width(); 04654 if (mat.length()!=x.length() || mat.width()!=x.width() || y.length()!=mat.length() || y.width()!=mat.width()) 04655 PLERROR("multiplyAcc this has size=%dx%d, x is %dx%d, y is %dx%d", 04656 mat.length(),mat.width(),x.length(),x.width(),y.length(),y.width()); 04657 T* p=mat.data(); 04658 T* xp=x.data(); 04659 T* yp=y.data(); 04660 for (int i=0;i<n;i++) 04661 p[i] += xp[i] * yp[i]; 04662 } 04663 04664 // result[i,j] += x[i,j]*x[i,j]*scale; 04665 template<class T> 04666 void squareMultiplyAcc(const TMat<T>& mat, const TMat<T>& x, T scale) 04667 { 04668 int n=x.length()*x.width(); 04669 if (mat.length()*mat.width()!=n) 04670 PLERROR("squareMultiplyAcc this has size=%d and x has size=%d", 04671 mat.width()*mat.length(),n); 04672 T* p=mat.data(); 04673 T* xp=x.data(); 04674 for (int i=0;i<n;i++) 04675 { 04676 T xpi = xp[i]; 04677 p[i] += scale * xpi * xpi; 04678 } 04679 } 04680 04681 // result[i,j] += (x[i,j]-y[i,j])^2*scale; 04682 template<class T> 04683 void diffSquareMultiplyAcc(const TMat<T>& mat, const TMat<T>& x, const TMat<T>& y, T scale) 04684 { 04685 int n=x.length()*x.width(); 04686 if (mat.length()*mat.width()!=n) 04687 PLERROR("diffSquareMultiplyAcc this has size=%d and x has size=%d", 04688 mat.width()*mat.length(),n); 04689 T* p=mat.data(); 04690 T* xp=x.data(); 04691 T* yp=y.data(); 04692 for (int i=0;i<n;i++) 04693 { 04694 T diff = (xp[i]-yp[i]); 04695 p[i] += scale * diff * diff; 04696 } 04697 } 04698 04700 // swapRows // 04705 template<class T> 04706 void swapRows(const TMat<T>& mat, int i, int j) 04707 { 04708 if (i == j) 04709 return; 04710 mat.swapRows(i, j); 04711 } 04712 04714 // selectAndOrder // 04716 template<class T> 04717 TVec<T> selectAndOrder(const TMat<T>& mat, int pos, int colnum=0) 04718 { 04719 #ifdef BOUNDCHECK 04720 if (colnum<0 || colnum>=mat.width()) PLERROR("Bad column number (%d)", colnum); 04721 if (pos<0 || pos>=mat.length()) PLERROR("Bad position (%d)", pos); 04722 #endif 04723 04724 int l=0; 04725 int h=mat.length()-1; 04726 TMat<T> v = mat.column(colnum); 04727 04728 while (l<h) 04729 { 04730 T p = v((l+h)/2,0); 04731 int x = l; 04732 int y = h; 04733 04734 do 04735 { 04736 while (v(x,0)<p) x++; 04737 while (p<v(y,0)) y--; 04738 if (x<=y) 04739 { 04740 mat.swapRows(x,y); 04741 x++; 04742 y--; 04743 } 04744 } while (x<=y); 04745 04746 if (pos>=x) l=x; 04747 else h=x-1; 04748 } 04749 04750 return mat(l); 04751 } 04752 04753 04754 // result[i,i] += lambda 04755 template<class T> 04756 void addToDiagonal(const TMat<T>& mat, T lambda) 04757 { 04758 T *d = mat.data(); 04759 int l=mat.length(); 04760 for (int i=0;i<l;i++,d+=mat.mod()+1) *d+=lambda; 04761 } 04762 04763 04764 04765 // result[i,i] += lambda[i] 04766 04767 template<class T> 04768 void addToDiagonal(const TMat<T>& mat, const TVec<T>& lambda) 04769 { 04770 #ifdef BOUNDCHECK 04771 if (lambda.length()!=mat.length()) 04772 PLERROR("Mat(%d)::addToDiagonal(Vec(%d)) inconsistent lengths", 04773 mat.length(), lambda.length()); 04774 #endif 04775 T *l = lambda.data(); 04776 T *d = mat.data(); 04777 int le= mat.length(); 04778 for (int i=0;i<le;i++,d+=mat.mod()+1,l++) *d += *l; 04779 } 04780 04782 template<class T> 04783 void fillDiagonal(const TMat<T>& mat, T val) 04784 { 04785 int l=mat.length(); 04786 for (int i=0;i<l;i++) 04787 mat(i,i) = val; 04788 } 04789 04791 template<class T> 04792 void fillDiagonal(const TMat<T>& mat, const TVec<T>& v) 04793 { 04794 int l=mat.length(); 04795 for (int i=0;i<l;i++) 04796 mat(i,i) = v[i]; 04797 } 04798 04799 04801 template<class T> 04802 void diag(const TMat<T>& mat, const TVec<T>& d) 04803 { 04804 T* d_ = d.data(); 04805 int l=mat.length(); 04806 for (int i=0;i<l;i++) 04807 d_[i] = mat(i,i); 04808 } 04809 04810 template<class T> 04811 TVec<T> diag(const TMat<T>& mat) 04812 { 04813 TVec<T> d(mat.length()); 04814 diag(mat, d); 04815 return d; 04816 } 04817 04818 template<class T> 04819 void diagonalOfSquare(const TMat<T>& mat, const TVec<T>& d) 04820 { 04821 T* d_=d.data(); 04822 for (int i=0;i<mat.length();i++) 04823 d_[i]=pownorm(mat(i)); 04824 } 04825 04826 04827 template<class T> 04828 void projectOnOrthogonalSubspace(const TMat<T>& mat, TMat<T> orthonormal_subspace) 04829 { 04830 for (int i=0;i<mat.length();i++) 04831 { 04832 TVec<T> row_i = mat(i); 04833 projectOnOrthogonalSubspace(row_i, orthonormal_subspace); 04834 } 04835 } 04836 04837 04838 template<class T> 04839 void averageAcrossRowsAndColumns(const TMat<T>& mat, TVec<T>& avg_across_rows, TVec<T>& avg_across_columns, bool ignored) 04840 { 04841 avg_across_rows.resize(mat.width()); 04842 avg_across_columns.resize(mat.length()); 04843 avg_across_rows.clear(); 04844 avg_across_columns.clear(); 04845 T* row_i=mat.data(); 04846 int w=mat.width(); 04847 for (int i=0;i<mat.length();i++) 04848 { 04849 T& avg_cols_i=avg_across_columns[i]; 04850 T* avg_rows = avg_across_rows.data(); 04851 for (int j=0;j<w;j++) 04852 { 04853 T row_ij=row_i[j]; 04854 avg_cols_i += row_ij; 04855 avg_rows[j] += row_ij; 04856 } 04857 row_i+=mat.mod(); 04858 } 04859 avg_across_rows /= mat.length(); 04860 avg_across_columns /= mat.width(); 04861 } 04862 04863 04864 template<class T> 04865 void addToRows(const TMat<T>& mat, const TVec<T> row, bool ignored) 04866 { 04867 int l=mat.length(); 04868 for (int i=0;i<l;i++) 04869 { 04870 TVec<T> row_i = mat(i); 04871 row_i += row; 04872 } 04873 } 04874 04875 04876 template<class T> 04877 void addToColumns(const TMat<T>& mat, const TVec<T> col, bool ignored) 04878 { 04879 T* row_i=mat.data(); 04880 int w=mat.width(); 04881 for (int i=0;i<mat.length();i++) 04882 { 04883 T col_i=col[i]; 04884 for (int j=0;j<w;j++) 04885 row_i[j] += col_i; 04886 row_i+=mat.mod(); 04887 } 04888 } 04889 04890 template<class T> 04891 void substractFromRows(const TMat<T>& mat, const TVec<T> row, bool ignored) 04892 { 04893 for (int i=0;i<mat.length();i++) 04894 { 04895 TVec<T> row_i = mat(i); 04896 row_i -= row; 04897 } 04898 } 04899 04900 04901 04902 // Probably bugged!!! 04903 template<class T> 04904 void substractFromColumns(const TMat<T>& mat, const TVec<T> col, bool ignored) 04905 { 04906 T* row_i=mat.data(); 04907 int w=mat.width(); 04908 for (int i=0;i<mat.length();i++) 04909 { 04910 T col_i=col[i]; 04911 for (int j=0;j<w;j++) 04912 row_i[j] -= col_i; 04913 row_i+=mat.mod(); 04914 } 04915 } 04916 04917 04918 template<class T> 04919 void addToMat(const TMat<T>& mat, T scalar, bool ignored) 04920 { mat += scalar; } 04921 04922 04923 // -------------- taken and adapted from Mat_maths.cc ------------------ 04924 04927 template<class T> 04928 T sum(const TMat<T>& mat, bool ignore_missing) 04929 { 04930 double res = 0.0; 04931 T* m_i = mat.data(); 04932 int w=mat.width(); 04933 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 04934 { 04935 for(int j=0; j<w; j++) 04936 { 04937 if (!is_missing(m_i[j])) res += m_i[j]; 04938 else if (!ignore_missing) return MISSING_VALUE; 04939 } 04940 } 04941 return T(res); 04942 } 04943 04946 template<class T> 04947 T sum(const TMat<T>& mat) 04948 { 04949 T res = T(0); 04950 T* m_i = mat.data(); 04951 int w=mat.width(); 04952 04953 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 04954 for(int j=0; j<w; j++) 04955 res += m_i[j]; 04956 return res; 04957 } 04958 04959 template<class T> 04960 T product(const TMat<T>& mat) 04961 { 04962 double res = 1.0; 04963 T* m_i = mat.data(); 04964 int w=mat.width(); 04965 04966 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 04967 for(int j=0; j<w; j++) 04968 res *= m_i[j]; 04969 return T(res); 04970 } 04971 04972 template<class T> 04973 T sum_of_squares(const TMat<T>& mat) 04974 { 04975 double res = 0.0; 04976 T* m_i = mat.data(); 04977 int w=mat.width(); 04978 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 04979 for(int j=0; j<w; j++) 04980 { 04981 T v = m_i[j]; 04982 res += v*v; 04983 } 04984 return T(res); 04985 } 04986 04987 template<class T> 04988 T mean(const TMat<T>& mat) 04989 { 04990 #ifdef BOUNDCHECK 04991 if(mat.length()==0 || mat.width()==0) 04992 PLERROR("IN T mean(const TMat<T>& mat) mat has 0 size"); 04993 #endif 04994 double res = 0.0; 04995 T* m_i = mat.data(); 04996 int w=mat.width(); 04997 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 04998 for(int j=0; j<w; j++) 04999 res += m_i[j]; 05000 return T(res/(mat.length()*mat.width())); 05001 } 05002 05003 template<class T> 05004 T geometric_mean(const TMat<T>& mat) 05005 { 05006 #ifdef BOUNDCHECK 05007 if(mat.length()==0 || mat.width()==0) 05008 PLERROR("IN T geometric_mean(const TMat<T>& mat) mat has 0 size"); 05009 #endif 05010 double res = 0.0; 05011 T* m_i = mat.data(); 05012 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05013 for(int j=0; j<mat.width(); j++) 05014 { 05015 T mij = m_i[j]; 05016 if (mij<=0) 05017 PLERROR("geometric_mean(TMat<T>): argument %g <=0 at position (%d,%d)", 05018 mij,i,j); 05019 res += pl_log(m_i[j]); 05020 } 05021 return T(exp(res/(mat.length()*mat.width()))); 05022 } 05023 05024 template<class T> 05025 T variance(const TMat<T>& mat, T meanval) 05026 { 05027 #ifdef BOUNDCHECK 05028 if(mat.length()==0 || mat.width()==0) 05029 PLERROR("IN T variance(const TMat<T>& mat, T meanval) mat has 0 size"); 05030 #endif 05031 double res = 0.0; 05032 T* m_i = mat.data(); 05033 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05034 for(int j=0; j<mat.width(); j++) 05035 { 05036 double diff = m_i[j]-meanval; 05037 res += diff*diff; 05038 } 05039 return res/(mat.length()*mat.width()-1); 05040 } 05041 05042 template<class T> 05043 T correlation(const TMat<T>& mat) 05044 { 05045 int n = mat.length(); 05046 #ifdef BOUNDCHECK 05047 if(n==0 || mat.width()==0) 05048 PLERROR("In T correlation(const TMat<T>& mat) mat has 0 size"); 05049 #endif 05050 if (mat.width() != 2) 05051 PLERROR("In T correlation(const TMat<T>& mat), mat width (%d) must be 2", mat.width()); 05052 05053 double s_x=0, s_y=0, s_xy=0, s_x2=0, s_y2=0; 05054 for (int i=0; i<n; i++) 05055 { 05056 T x = mat(i,0); 05057 T y = mat(i,1); 05058 s_x += x; 05059 s_x2 += x*x; 05060 s_y += y; 05061 s_y2 += y*y; 05062 s_xy += x*y; 05063 } 05064 05065 return (n*s_xy - s_x*s_y)/sqrt((n*s_x2 - s_x*s_x)*(n*s_y2 - s_y*s_y)); 05066 } 05067 05068 template<class T> 05069 T correlation(const TVec<T>& x, const TVec<T>& y) 05070 { 05071 int n = x.length(); 05072 #ifdef BOUNDCHECK 05073 if(n==0 || y.length()==0) 05074 PLERROR("In T correlation(const TVec<T>& x, const TVec<T>& y), one Vec has 0 size"); 05075 #endif 05076 if (n != y.length()) 05077 PLERROR("In T correlation(const TVec<T>& x, const TVec<T>& y), both Vec must have same length (%d != %d)", n, y.length()); 05078 05079 double s_x=0, s_y=0, s_xy=0, s_x2=0, s_y2=0; 05080 for (int i=0; i<n; i++) 05081 { 05082 T x_val = x[i]; 05083 T y_val = y[i]; 05084 s_x += x_val; 05085 s_x2 += x_val*x_val; 05086 s_y += y_val; 05087 s_y2 += y_val*y_val; 05088 s_xy += x_val*y_val; 05089 } 05090 05091 return (n*s_xy - s_x*s_y)/sqrt((n*s_x2 - s_x*s_x)*(n*s_y2 - s_y*s_y)); 05092 } 05093 05095 template<class T> 05096 T min(const TMat<T>& mat) 05097 { 05098 #ifdef BOUNDCHECK 05099 if(mat.length()==0 || mat.width()==0) 05100 PLERROR("IN T min(const TMat<T>& mat) mat has 0 size"); 05101 #endif 05102 T* m_i = mat.data(); 05103 double minval = m_i[0]; 05104 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05105 for(int j=0; j<mat.width(); j++) 05106 if(m_i[j]<minval) 05107 minval = m_i[j]; 05108 return minval; 05109 } 05110 05112 template<class T> 05113 T min(const TMat<T>& mat, int& min_i, int& min_j) 05114 { 05115 PLASSERT(mat.size() != 0); 05116 05117 T* m_i = mat.data(); 05118 double minval = m_i[0]; 05119 min_i = 0; 05120 min_j = 0; 05121 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05122 for(int j=0; j<mat.width(); j++) 05123 if(m_i[j]<minval) 05124 { 05125 minval = m_i[j]; 05126 min_i = i; 05127 min_j = j; 05128 } 05129 return minval; 05130 } 05131 05133 template<class T> 05134 T max(const TMat<T>& mat) 05135 { 05136 #ifdef BOUNDCHECK 05137 if(mat.length()==0 || mat.width()==0) 05138 PLERROR("IN T max(const TMat<T>& mat) mat has 0 size"); 05139 #endif 05140 T* m_i = mat.data(); 05141 double maxval = m_i[0]; 05142 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05143 for(int j=0; j<mat.width(); j++) 05144 if(m_i[j]>maxval) 05145 maxval = m_i[j]; 05146 return maxval; 05147 } 05148 05150 template<class T> 05151 T max(const TMat<T>& mat, int& max_i, int& max_j) 05152 { 05153 PLASSERT(mat.size() != 0); 05154 05155 T* m_i = mat.data(); 05156 double maxval = m_i[0]; 05157 max_i = 0; 05158 max_j = 0; 05159 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05160 for(int j=0; j<mat.width(); j++) 05161 if(m_i[j]>maxval) 05162 { 05163 maxval = m_i[j]; 05164 max_i = i; 05165 max_j = j; 05166 } 05167 return maxval; 05168 } 05169 05171 template<class T> 05172 T minabs(const TMat<T>& mat) 05173 { 05174 #ifdef BOUNDCHECK 05175 if(mat.length()==0 || mat.width()==0) 05176 PLERROR("IN T minabs(const TMat<T>& mat) mat has 0 size"); 05177 #endif 05178 T* m_i = mat.data(); 05179 double minval = fabs(m_i[0]); 05180 int w=mat.width(); 05181 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05182 for(int j=0; j<w; j++) 05183 { 05184 T a=fabs(m_i[j]); 05185 if(a<minval) 05186 minval = a; 05187 } 05188 return minval; 05189 } 05190 05192 template<class T> 05193 T minabs(const TMat<T>& mat, int& min_i, int& min_j) 05194 { 05195 PLASSERT(mat.size() != 0); 05196 05197 T* m_i = mat.data(); 05198 double minval = fabs(m_i[0]); 05199 min_i = 0; 05200 min_j = 0; 05201 int w=mat.width(); 05202 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05203 for(int j=0; j<w; j++) 05204 { 05205 T a = fabs(m_i[j]); 05206 if(a<minval) 05207 { 05208 minval = a; 05209 min_i = i; 05210 min_j = j; 05211 } 05212 } 05213 return minval; 05214 } 05215 05217 template<class T> 05218 T maxabs(const TMat<T>& mat) 05219 { 05220 #ifdef BOUNDCHECK 05221 if(mat.length()==0 || mat.width()==0) 05222 PLERROR("IN T maxabs(const TMat<T>& mat) mat has 0 size"); 05223 #endif 05224 T* m_i = mat.data(); 05225 double maxval = fabs(m_i[0]); 05226 int w=mat.width(); 05227 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05228 for(int j=0; j<w; j++) 05229 { 05230 T a=fabs(m_i[j]); 05231 if(a>maxval) 05232 maxval = a; 05233 } 05234 return maxval; 05235 } 05236 05238 template<class T> 05239 T maxabs(const TMat<T>& mat, int& max_i, int& max_j) 05240 { 05241 PLASSERT(mat.size() != 0); 05242 05243 T* m_i = mat.data(); 05244 double maxval = fabs(m_i[0]); 05245 max_i = 0; 05246 max_j = 0; 05247 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05248 for(int j=0; j<mat.width(); j++) 05249 { 05250 T a = fabs(m_i[j]); 05251 if(a>maxval) 05252 { 05253 maxval = a; 05254 max_i = i; 05255 max_j = j; 05256 } 05257 } 05258 return maxval; 05259 } 05260 05262 template<class T> 05263 void argmin(const TMat<T>& mat, int& mini, int& minj) 05264 { 05265 #ifdef BOUNDCHECK 05266 if(mat.length()==0 || mat.width()==0) 05267 PLERROR("IN void argmin(const TMat<T>& mat, int& mini, iny& minj) mat has 0 size"); 05268 #endif 05269 T* m_i = mat.data(); 05270 mini=0; 05271 minj=0; 05272 double minval = m_i[0]; 05273 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05274 for(int j=0; j<mat.width(); j++) 05275 if(m_i[j]<minval) 05276 { 05277 minval = m_i[j]; 05278 mini = i; 05279 minj = j; 05280 } 05281 } 05282 05283 // Same as above with the max. 05284 template<class T> 05285 void argmax(const TMat<T>& mat, int& maxi, int& maxj) 05286 { 05287 #ifdef BOUNDCHECK 05288 if(mat.length()==0 || mat.width()==0) 05289 PLERROR("IN void argmax(const TMat<T>& mat, int& maxi, iny& maxj) mat has 0 size"); 05290 #endif 05291 T* m_i = mat.data(); 05292 maxi=0; 05293 maxj=0; 05294 double maxval = m_i[0]; 05295 for(int i=0; i<mat.length(); i++, m_i+=mat.mod()) 05296 for(int j=0; j<mat.width(); j++) 05297 if(m_i[j]>maxval) 05298 { 05299 maxval = m_i[j]; 05300 maxi = i; 05301 maxj = j; 05302 } 05303 } 05304 05306 template<class T> 05307 int argmin(const TMat<T>& m) 05308 { 05309 int imin, jmin; 05310 argmin(m,imin,jmin); 05311 return (imin*m.width()+jmin); 05312 } 05313 05315 template<class T> 05316 int argmax(const TMat<T>& m) 05317 { 05318 int imax, jmax; 05319 argmax(m,imax,jmax); 05320 return (imax*m.width()+jmax); 05321 } 05322 05328 // singlecolumn[i] = sum_j mat(j,i) 05329 template<class T> 05330 void rowSum(const TMat<T>& mat, const TMat<T>& singlecolumn) 05331 { 05332 #ifdef BOUNDCHECK 05333 if(singlecolumn.length()!=mat.length() || singlecolumn.width() != 1) 05334 PLERROR("IN void rowSum(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix"); 05335 #endif 05336 for(int i=0; i<mat.length(); i++) 05337 singlecolumn(i,0) = sum(mat(i)); 05338 } 05339 05340 // singlecolumn[i] += sum_j mat(j,i) 05341 template<class T> 05342 void rowSumAcc(const TMat<T>& mat, const TMat<T>& singlecolumn) 05343 { 05344 #ifdef BOUNDCHECK 05345 if(singlecolumn.length()!=mat.length() || singlecolumn.width() != 1) 05346 PLERROR("IN void rowSum(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix"); 05347 #endif 05348 for(int i=0; i<mat.length(); i++) 05349 singlecolumn(i,0) += sum(mat(i)); 05350 } 05351 05352 05353 template<class T> 05354 void rowSum(const TMat<T>& mat, const TVec<T>& colvec) 05355 { 05356 #ifdef BOUNDCHECK 05357 if(colvec.length()!=mat.length()) 05358 PLERROR("IN void rowSum(const TMat<T>& mat, const TVec<T>& colvec) colvec must have same length as mat"); 05359 #endif 05360 for(int i=0; i<mat.length(); i++) 05361 colvec[i] = sum(mat(i)); 05362 } 05363 05364 template<class T> 05365 void rowMean(const TMat<T>& mat, const TMat<T>& singlecolumn) 05366 { 05367 #ifdef BOUNDCHECK 05368 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 05369 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"); 05370 #endif 05371 for(int i=0; i<mat.length(); i++) 05372 singlecolumn(i,0) = mean(mat(i)); 05373 } 05374 05375 template<class T> 05376 void rowVariance(const TMat<T>& mat, const TMat<T>& singlecolumn, const TMat<T>& rowmean) 05377 { 05378 #ifdef BOUNDCHECK 05379 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || rowmean.length()!=mat.length() || rowmean.width()!=1 || mat.width()==0) 05380 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"); 05381 #endif 05382 for(int i=0; i<mat.length(); i++) 05383 singlecolumn(i,0) = variance(mat(i),rowmean(i,0)); 05384 } 05385 05386 template<class T> 05387 void rowSumOfSquares(const TMat<T>& mat, const TMat<T>& singlecolumn) 05388 { 05389 #ifdef BOUNDCHECK 05390 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1) 05391 PLERROR("IN void rowSumOfSquares(const TMat<T>& mat, TMat<T>& singlecolumn) singlecolumn must be a mat.length() x 1 matrix"); 05392 #endif 05393 int w=mat.width(); 05394 for (int i=0;i<mat.length();i++) 05395 { 05396 T ss=0; 05397 T* mi=mat[i]; 05398 for (int j=0;j<w;j++) { T mij=mi[j]; ss+=mij*mij; } 05399 singlecolumn(i,0)=ss; 05400 } 05401 } 05402 05403 template<class T> 05404 void rowMax(const TMat<T>& mat, const TMat<T>& singlecolumn) 05405 { 05406 #ifdef BOUNDCHECK 05407 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 05408 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"); 05409 #endif 05410 for(int i=0; i<mat.length(); i++) 05411 singlecolumn(i,0) = max(mat(i)); 05412 } 05413 05414 template<class T> 05415 void rowMax(const TMat<T>& mat, const TVec<T>& colvec) 05416 { 05417 #ifdef BOUNDCHECK 05418 if(colvec.length()!=mat.length()) 05419 PLERROR("IN void rowSum(const TMat<T>& mat, const TVec<T>& colvec) colvec must have same length as mat"); 05420 #endif 05421 for(int i=0; i<mat.length(); i++) 05422 colvec[i] = max(mat(i)); 05423 } 05424 05425 template<class T> 05426 void rowMin(const TMat<T>& mat, const TMat<T>& singlecolumn) 05427 { 05428 #ifdef BOUNDCHECK 05429 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 05430 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"); 05431 #endif 05432 for(int i=0; i<mat.length(); i++) 05433 singlecolumn(i,0) = min(mat(i)); 05434 } 05435 05436 05437 template<class T> 05438 void rowMin(const TMat<T>& mat, const TVec<T>& colvec) 05439 { 05440 #ifdef BOUNDCHECK 05441 if(colvec.length()!=mat.length()) 05442 PLERROR("IN void rowSum(const TMat<T>& mat, const TVec<T>& colvec) colvec must have same length as mat"); 05443 #endif 05444 for(int i=0; i<mat.length(); i++) 05445 colvec[i] = min(mat(i)); 05446 } 05447 05448 template<class T> 05449 void rowArgmax(const TMat<T>& mat, const TMat<T>& singlecolumn) 05450 { 05451 #ifdef BOUNDCHECK 05452 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 05453 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"); 05454 #endif 05455 for(int i=0; i<mat.length(); i++) 05456 singlecolumn(i,0) = argmax(mat(i)); 05457 } 05458 05459 template<class T> 05460 void rowArgmin(const TMat<T>& mat, const TMat<T>& singlecolumn) 05461 { 05462 #ifdef BOUNDCHECK 05463 if(singlecolumn.length()!=mat.length() || singlecolumn.width()!=1 || mat.width()==0) 05464 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"); 05465 #endif 05466 for(int i=0; i<mat.length(); i++) 05467 singlecolumn(i,0) = argmin(mat(i)); 05468 } 05469 05475 template<class T> 05476 void columnSum(const TMat<T>& mat, TVec<T>& result) 05477 { 05478 #ifdef BOUNDCHECK 05479 if(result.length()!=mat.width()) 05480 PLERROR("IN void columnSum(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat"); 05481 #endif 05482 int l = mat.length(); 05483 result << mat(0); 05484 for(int j=1; j<l; j++) 05485 result += mat(j); 05486 } 05487 05488 template<class T> 05489 void columnSumOfSquares(const TMat<T>& mat, TVec<T>& result) 05490 { 05491 #ifdef BOUNDCHECK 05492 if(result.length()!=mat.width()) 05493 PLERROR("IN void columnSumOfSquares(const TMat<T>& mat, TVec<T>& result) the length of result must equal the width of mat"); 05494 #endif 05495 for(int j=0; j<mat.width(); j++) 05496 result[j] = sum_of_squares(mat.column(j)); 05497 } 05498 05499 template<class T> 05500 void columnMean(const TMat<T>& mat, TVec<T>& result) 05501 { 05502 #ifdef BOUNDCHECK 05503 if(result.length()!=mat.width() || mat.length()==0) 05504 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"); 05505 #endif 05506 columnSum(mat,result); 05507 result /= real(mat.length()); 05508 } 05509 05510 template<class T> 05511 void columnWeightedMean(const TMat<T>& mat, TVec<T>& result) 05512 { 05513 #ifdef BOUNDCHECK 05514 if(result.length()!=mat.width()-1 || mat.length()<=1) 05515 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"); 05516 #endif 05517 TVec<T> column_j_vec(mat.length()), weights_vec(mat.length()); 05518 TMat<T> column_j_mat(mat.length(), 1), weights_mat(mat.length(), 1); 05519 for(int j=0; j<mat.width()-1; j++){ 05520 column_j_mat = mat.column(j); 05521 weights_mat = mat.column(mat.width()-1); 05522 column_j_vec = column_j_mat.toVecCopy(); 05523 weights_vec = weights_mat.toVecCopy(); 05524 result[j] = weighted_mean(column_j_vec, weights_vec); 05525 } 05526 } 05527 05528 template<class T> 05529 void columnVariance(const TMat<T>& mat, TVec<T>& result, const TVec<T>& columnmean) 05530 { 05531 #ifdef BOUNDCHECK 05532 if(result.length()!=mat.width() || columnmean.length()!=mat.width() || mat.length()==0) 05533 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"); 05534 #endif 05535 for(int j=0; j<mat.width(); j++) 05536 result[j] = variance(mat.column(j),columnmean[j]); 05537 } 05538 05539 template<class T> 05540 void columnWeightedVariance(const TMat<T>& mat, TVec<T>& result, const TVec<T>& column_weighted_mean) 05541 { 05542 #ifdef BOUNDCHECK 05543 if(result.length()!=mat.width()-1 || column_weighted_mean.length()!=mat.width()-1 || mat.length()<=1) 05544 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"); 05545 #endif 05546 T column_no_weighted_mean_j; 05547 TVec<T> column_j_vec(mat.length()), weights_vec(mat.length()); 05548 TMat<T> column_j_mat(mat.length(), 1), weights_mat(mat.length(), 1); 05549 for(int j=0; j<mat.width()-1; j++){ 05550 column_j_mat = mat.column(j); 05551 weights_mat = mat.column(mat.width()-1); 05552 column_j_vec = column_j_mat.toVecCopy(); 05553 weights_vec = weights_mat.toVecCopy(); 05554 column_no_weighted_mean_j = mean(mat.column(j)); 05555 result[j] = weighted_variance(column_j_vec, weights_vec, column_no_weighted_mean_j, column_weighted_mean[j]); 05556 } 05557 } 05558 05559 template<class T> 05560 void columnMax(const TMat<T>& mat, TVec<T>& result) 05561 { 05562 #ifdef BOUNDCHECK 05563 if(result.length()!=mat.width() || mat.length()==0) 05564 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"); 05565 #endif 05566 for(int j=0; j<mat.width(); j++) 05567 result[j] = max(mat.column(j)); 05568 } 05569 05570 template<class T> 05571 void columnMin(const TMat<T>& mat, TVec<T>& result) 05572 { 05573 #ifdef BOUNDCHECK 05574 if(result.length()!=mat.width() || mat.length()==0) 05575 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"); 05576 #endif 05577 for(int j=0; j<mat.width(); j++) 05578 result[j] = min(mat.column(j)); 05579 } 05580 05581 template<class T> 05582 void columnArgmax(const TMat<T>& mat, TVec<T>& result) 05583 { 05584 #ifdef BOUNDCHECK 05585 if(result.length()!=mat.width() || mat.length()==0) 05586 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"); 05587 #endif 05588 int imax, jmax; 05589 for(int j=0; j<mat.width(); j++) 05590 { 05591 argmax(mat.column(j), imax, jmax); 05592 result[j] = imax; 05593 } 05594 } 05595 05596 template<class T> 05597 void columnArgmin(const TMat<T>& mat, TVec<T>& result) 05598 { 05599 #ifdef BOUNDCHECK 05600 if(result.length()!=mat.width() || mat.length()==0) 05601 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"); 05602 #endif 05603 int imin, jmin; 05604 for(int j=0; j<mat.width(); j++) 05605 { 05606 argmin(mat.column(j), imin, jmin); 05607 result[j] = imin; 05608 } 05609 } 05610 05611 template<class T> 05612 T mahalanobis_distance(const TVec<T>& input, const TVec<T>& meanvec, const 05613 TMat<T>& inversecovmat) 05614 { 05615 TVec<T> diff = input-meanvec; 05616 return dot(diff,product(inversecovmat,diff)); 05617 } 05618 05620 template<class T> 05621 inline void computeMean(const TMat<T>& m, TVec<T>& meanvec) { columnMean(m,meanvec); } 05622 05624 template<class T> 05625 void computeMeanAndVariance(const TMat<T>& m, TVec<T>& meanvec, TVec<T>& variancevec) 05626 { 05627 columnMean(m,meanvec); 05628 columnVariance(m,variancevec,meanvec); 05629 } 05630 05637 template<class T> 05638 void computeInverseStandardDeviationFromMeanAndSquareMean(const TMat<T>& inverse_standard_deviation, 05639 const TMat<T>& means, 05640 const TMat<T>& mean_of_squares, 05641 real default_stddev = 1, 05642 real min_stddev = -1) 05643 { 05644 int n=inverse_standard_deviation.length(); 05645 int m=inverse_standard_deviation.width(); 05646 int invs_mod = inverse_standard_deviation.mod(); 05647 int mu_mod = means.mod(); 05648 int mu2_mod = mean_of_squares.mod(); 05649 #ifdef BOUNDCHECK 05650 if (means.length()!=n || means.width()!=m || mean_of_squares.length()!=n 05651 || mean_of_squares.width()!=m) 05652 PLERROR("In computeInverseStandardDeviationFromMeanAndSquareMean - Arguments have incompatible sizes"); 05653 #endif 05654 T* invs = inverse_standard_deviation.data(); 05655 T* mu = means.data(); 05656 T* mu2 = mean_of_squares.data(); 05657 for (int i=0;i<n;i++, invs += invs_mod, mu += mu_mod, mu2 += mu2_mod) { 05658 for (int j=0;j<m;j++) 05659 { 05660 real diff = mu2[j] - mu[j] * mu[j]; 05661 if (diff>0) { 05662 real sqrt_diff = sqrt(diff); 05663 if (sqrt_diff < min_stddev) // NB: Cannot happen if 'min_stddev' is -1. 05664 invs[j] = real(1.0 / default_stddev); 05665 else 05666 invs[j] = real(1.0 / sqrt_diff); 05667 } 05668 else { 05669 if (min_stddev < 0) 05670 // No minimum standard deviation provided, this is suspect. 05671 PLWARNING("In computeInverseStandardDeviationFromMeanAndSquareMean - Variance is not > 0"); 05672 invs[j] = real(1.0 / default_stddev); 05673 } 05674 } 05675 } 05676 } 05677 05678 05679 05680 template<class T> 05681 void computeCovar(const TMat<T>& m, const TVec<T>& meanvec, TMat<T>& covarmat) 05682 { 05683 int n = m.width(); 05684 covarmat.resize(n,n); 05685 transposeProduct(covarmat,m,m); 05686 covarmat /= T(m.length()); 05687 externalProductScaleAcc(covarmat,meanvec,meanvec,T(-1)); 05688 } 05689 05690 template<class T> 05691 void computeMeanAndCovar(const TMat<T>& m, TVec<T>& meanvec, TMat<T>& covarmat) 05692 { 05693 int n = m.width(); 05694 meanvec.resize(n); 05695 covarmat.resize(n,n); 05696 columnMean(m,meanvec); 05697 05698 transposeProduct(covarmat,m,m); 05699 covarmat /= T(m.length()); 05700 externalProductScaleAcc(covarmat,meanvec,meanvec,T(-1)); 05701 05702 /* 05703 Mat mm = m.copy(); 05704 mm -= meanvec; 05705 transposeProduct(covarmat,mm,mm); 05706 covarmat /= T(m.length()); 05707 */ 05708 } 05709 05711 template<class T> 05712 void computeMeanAndStddev(const TMat<T>& m, TVec<T>& meanvec, TVec<T>& stddevvec) 05713 { 05714 columnMean(m,meanvec); 05715 columnVariance(m,stddevvec,meanvec); 05716 int l=stddevvec.length(); 05717 for(int i=0; i<l; i++) 05718 stddevvec[i] = sqrt(stddevvec[i]); 05719 } 05720 05721 05724 template<class T> 05725 void computeColumnsMeanAndStddev(const TMat<T>& m, TMat<T>& meanvec, TMat<T>& stddevvec) 05726 { 05727 rowMean(m,meanvec); 05728 rowVariance(m,stddevvec,meanvec); 05729 int l=stddevvec.length(); 05730 for(int i=0; i<l; i++) 05731 stddevvec[i][0] = sqrt(stddevvec[i][0]); 05732 } 05733 05735 template<class T> 05736 void normalize(TMat<T>& m) 05737 { 05738 TVec<T> meanvec(m.width()); 05739 TVec<T> stddevvec(m.width()); 05740 computeMeanAndStddev(m,meanvec,stddevvec); 05741 m -= meanvec; 05742 m /= stddevvec; 05743 } 05744 05746 template<class T> 05747 void normalizeRows(const TMat<T>& m) 05748 { 05749 int l = m.length(); 05750 for(int i=0; i<l; i++) 05751 { 05752 TVec<T> v = m(i); 05753 v /= sum(v); 05754 } 05755 } 05756 05758 template<class T> 05759 void normalizeColumns(const TMat<T>& m) 05760 { 05761 int w = m.width(); 05762 for(int j=0; j<w; j++) 05763 { 05764 TMat<T> v = m.column(j); 05765 v /= sum(v); 05766 } 05767 } 05768 05770 template<class T> 05771 void normalize(TMat<T>& m, double n) 05772 { 05773 for(int i=0; i<m.length(); i++) 05774 { 05775 TVec<T> m_i = m(i); 05776 normalize(m_i,n); 05777 } 05778 } 05779 05780 template<class T> 05781 void operator+=(const TMat<T>& m, T scalar) 05782 { 05783 T* m_i = m.data(); 05784 int w = m.width(); 05785 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 05786 for(int j=0; j<w; j++) 05787 m_i[j] += scalar; 05788 } 05789 05790 template<class T> 05791 void operator*=(const TMat<T>& m, T scalar) 05792 { 05793 T* m_i = m.data(); 05794 int w = m.width(); 05795 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 05796 for(int j=0; j<w; j++) 05797 m_i[j] *= scalar; 05798 } 05799 05800 template<class T> 05801 inline void operator-=(const TMat<T>& m, T scalar) { m += (-scalar); } 05802 05803 template<class T> 05804 inline void operator/=(const TMat<T>& m, T scalar) { m *= (T(1)/scalar); } 05805 05806 template<class T> 05807 inline void operator/=(const TMat<T>& m, int scalar) { m *= (T(1)/scalar); } 05808 05809 05811 template<class T> 05812 void operator+=(const TMat<T>& m, const TVec<T>& v) 05813 { 05814 #ifdef BOUNDCHECK 05815 if(m.width()!=v.length()) 05816 PLERROR("IN operator+=(const TMat<T>& m, const TVec<T>& v) v must be as long as m is wide"); 05817 #endif 05818 T* m_i = m.data(); 05819 T* vv = v.data(); 05820 int w=m.width(); 05821 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 05822 for(int j=0; j<w; j++) 05823 m_i[j] += vv[j]; 05824 } 05825 05827 template<class T> 05828 void operator-=(const TMat<T>& m, const TVec<T>& v) 05829 { 05830 #ifdef BOUNDCHECK 05831 if(m.width()!=v.length()) 05832 PLERROR("IN operator-=(const TMat<T>& m, const TVec<T>& v) v must be as long as m is wide"); 05833 #endif 05834 T* m_i = m.data(); 05835 T* vv = v.data(); 05836 int w=m.width(); 05837 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 05838 for(int j=0; j<w; j++) 05839 m_i[j] -= vv[j]; 05840 } 05841 05843 template<class T> 05844 void operator*=(const TMat<T>& m, const TVec<T>& v) 05845 { 05846 #ifdef BOUNDCHECK 05847 if(m.width()!=v.length()) 05848 PLERROR("IN operator*=(const TMat<T>& m, const TVec<T>& v) v must be as long as m is wide"); 05849 #endif 05850 T* m_i = m.data(); 05851 T* vv = v.data(); 05852 int w=m.width(); 05853 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 05854 for(int j=0; j<w; j++) 05855 m_i[j] *= vv[j]; 05856 } 05857 05859 template<class T> 05860 void operator*=(const TMat<T>& m1, const TMat<T>& m2) 05861 { 05862 int n=m1.length(); 05863 int l=m1.width(); 05864 #ifdef BOUNDCHECK 05865 if(l!=m2.width() || n!=m2.length()) 05866 PLERROR("IN operator*=(const TMat<T>& m1(%d,%d), const TMat<T>& m2(%d,%d)) sizes differ", 05867 m1.length(),m1.width(),m2.length(),m2.width()); 05868 #endif 05869 T* m1_i = m1.data(); 05870 T* m2_i = m2.data(); 05871 for(int i=0; i<n; i++, m1_i+=m1.mod(),m2_i+=m2.mod()) 05872 for(int j=0; j<l; j++) 05873 m1_i[j] *= m2_i[j]; 05874 } 05875 05877 template<class T> 05878 void operator/=(const TMat<T>& m, const TVec<T>& v) 05879 { 05880 #ifdef BOUNDCHECK 05881 if(m.width()!=v.length()) 05882 PLERROR("IN operator/=(const TMat<T>& m, const TVec<T>& v) v must be as long as m is wide"); 05883 #endif 05884 T* m_i = m.data(); 05885 T* vv = v.data(); 05886 int w=m.width(); 05887 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 05888 for(int j=0; j<w; j++) 05889 m_i[j] /= vv[j]; 05890 } 05891 05893 template<class T> 05894 void operator/=(const TMat<T>& m1, const TMat<T>& m2) 05895 { 05896 int n=m1.length(); 05897 int l=m1.width(); 05898 #ifdef BOUNDCHECK 05899 if(l!=m2.width() || n!=m2.length()) 05900 PLERROR("IN operator/=(const TMat<T>& m1(%d,%d), const TMat<T>& m2(%d,%d)) sizes differ", 05901 m1.length(),m1.width(),m2.length(),m2.width()); 05902 #endif 05903 T* m1_i = m1.data(); 05904 T* m2_i = m2.data(); 05905 for(int i=0; i<n; i++, m1_i+=m1.mod(),m2_i+=m2.mod()) 05906 for(int j=0; j<l; j++) 05907 m1_i[j] /= m2_i[j]; 05908 } 05909 05910 template<class T> 05911 void operator+=(const TMat<T>& m1, const TMat<T>& m2) 05912 { 05913 int n=m1.length(); 05914 int l=m1.width(); 05915 #ifdef BOUNDCHECK 05916 if(m1.width()!=m2.width() || m1.length()!=m2.length()) 05917 PLERROR("IN operator+=(const TMat<T>& m1(%d,%d), const TMat<T>& m2(%d,%d)): m1 and m2 must have same dimensions", 05918 m1.length(),m1.width(),m2.length(),m2.width()); 05919 #endif 05920 T* m1_i = m1.data(); 05921 T* m2_i = m2.data(); 05922 for(int i=0; i<n; i++, m1_i+=m1.mod(),m2_i+=m2.mod()) 05923 for(int j=0; j<l; j++) 05924 m1_i[j] += m2_i[j]; 05925 } 05926 05927 template<class T> 05928 void operator-=(const TMat<T>& m1, const TMat<T>& m2) 05929 { 05930 int n=m1.length(); 05931 int l=m1.width(); 05932 #ifdef BOUNDCHECK 05933 if(m1.width()!=m2.width() || m1.length()!=m2.length()) 05934 PLERROR("IN operator+=(const TMat<T>& m1(%d,%d), const TMat<T>& m2(%d,%d)): m1 and m2 must have same dimensions", 05935 m1.length(),m1.width(),m2.length(),m2.width()); 05936 #endif 05937 if(m1.isNotEmpty()) // calc only if some data 05938 { 05939 T* m1_i = m1.data(); 05940 T* m2_i = m2.data(); 05941 for(int i=0; i<n; i++, m1_i+=m1.mod(),m2_i+=m2.mod()) 05942 for(int j=0; j<l; j++) 05943 m1_i[j] -= m2_i[j]; 05944 } 05945 } 05946 05947 template<class T> 05948 TMat<T> operator-(const TMat<T>& m1, const TMat<T>& m2) 05949 { 05950 TMat<T> result(m1.length(), m1.width()); 05951 substract(m1,m2,result); 05952 return result; 05953 } 05954 05955 template<class T> 05956 TMat<T> operator+(const TMat<T>& m1, const TMat<T>& m2) 05957 { 05958 TMat<T> result(m1.length(), m1.width()); 05959 add(m1,m2,result); 05960 return result; 05961 } 05962 05963 template<class T> 05964 void substract(const TMat<T>& m1, const TMat<T>& m2, TMat<T>& destination) 05965 { 05966 #ifdef BOUNDCHECK 05967 if(m1.width()!=m2.width() || m1.length()!=m2.length() 05968 || m1.width()!=destination.width() || m1.length()!=destination.length()) 05969 PLERROR("IN substract(m1(%d,%d), m2(%d,%d), dest(%d,%d)): args must have same dimensions", 05970 m1.length(),m1.width(),m2.length(),m2.width(),destination.length(), 05971 destination.width()); 05972 #endif 05973 T* m1_i = m1.data(); 05974 T* m2_i = m2.data(); 05975 T* d_i = destination.data(); 05976 int m1_mod = m1.mod(); 05977 int m2_mod = m2.mod(); 05978 int d_mod = destination.mod(); 05979 int w = m1.width(); 05980 for (int i=0;i<m1.length();i++,m1_i+=m1_mod,m2_i+=m2_mod,d_i+=d_mod) 05981 for (int j=0;j<w;j++) 05982 d_i[j] = m1_i[j] - m2_i[j]; 05983 } 05984 05985 template<class T> 05986 void add(const TMat<T>& m1, const TMat<T>& m2, TMat<T>& destination) 05987 { 05988 #ifdef BOUNDCHECK 05989 if(m1.width()!=m2.width() || m1.length()!=m2.length() 05990 || m1.width()!=destination.width() || m1.length()!=destination.length()) 05991 PLERROR("IN substract(m1(%d,%d), m2(%d,%d), dest(%d,%d)): args must have same dimensions", 05992 m1.length(),m1.width(),m2.length(),m2.width(),destination.length(), 05993 destination.width()); 05994 #endif 05995 T* m1_i = m1.data(); 05996 T* m2_i = m2.data(); 05997 T* d_i = destination.data(); 05998 int m1_mod = m1.mod(); 05999 int m2_mod = m2.mod(); 06000 int d_mod = destination.mod(); 06001 int w = m1.width(); 06002 for (int i=0;i<m1.length();i++,m1_i+=m1_mod,m2_i+=m2_mod,d_i+=d_mod) 06003 for (int j=0;j<w;j++) 06004 d_i[j] = m1_i[j] + m2_i[j]; 06005 } 06006 06008 template<class T> 06009 TMat<T> operator-(const TMat<T>& m) 06010 { 06011 TMat<T> opposite(m.length(),m.width()); 06012 T *m_i=m.data(); 06013 T *o_i=opposite.data(); 06014 int w=m.width(); 06015 for (int i=0;i<m.length();i++,m_i+=m.mod(),o_i+=opposite.mod()) 06016 for (int j=0;j<w;j++) 06017 o_i[j] = - m_i[j]; 06018 return opposite; 06019 } 06020 06022 template<class T> 06023 void negateElements(const TMat<T>& m) 06024 { 06025 T* m_i = m.data(); 06026 int w=m.width(); 06027 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 06028 for(int j=0; j<w; j++) 06029 m_i[j] = -m_i[j]; 06030 } 06031 06033 template<class T> 06034 void invertElements(const TMat<T>& m) 06035 { 06036 T* m_i = m.data(); 06037 int w=m.width(); 06038 for(int i=0; i<m.length(); i++, m_i+=m.mod()) 06039 for(int j=0; j<w; j++) 06040 m_i[j] = 1.0/m_i[j]; 06041 } 06042 06043 // result * m = identity 06044 // (works only if m.length()>=m.width()) 06045 template<class T> 06046 TMat<T> leftPseudoInverse(TMat<T>& m) 06047 { 06048 TMat<T> inv(m.width(), m.length()); 06049 leftPseudoInverse(m,inv); 06050 return inv; 06051 } 06052 06053 // result * m = identity 06054 // (works only if m.length()>=m.width()) 06055 template<class T> 06056 void leftPseudoInverse(const TMat<T>& m, TMat<T>& inv) 06057 { 06058 if (m.length()==m.width()) 06059 inverse(m,inv); 06060 if (m.length()<m.width()) 06061 PLERROR("leftPseudoInverse: matrix length(%d) must be >= width(%d)", 06062 m.length(), m.width()); 06063 PLERROR("SVD not implemented yet"); 06064 } 06065 06066 // m * result = identity 06067 // (works only if m.length()<=m.width()) 06068 template<class T> 06069 TMat<T> rightPseudoInverse(TMat<T>& m) 06070 { 06071 TMat<T> inv(m.width(), m.length()); 06072 rightPseudoInverse(m,inv); 06073 return inv; 06074 } 06075 06076 // m * result = identity 06077 // (works only if m.length()<=m.width()) 06078 template<class T> 06079 void rightPseudoInverse(const TMat<T>& m, TMat<T>& inv) 06080 { 06081 if (m.length()==m.width()) 06082 inverse(m,inv); 06083 if (m.length()>m.width()) 06084 PLERROR("rightPseudoInverse: matrix length(%d) must be <= width(%d)", 06085 m.length(), m.width()); 06086 PLERROR("SVD not implemented yet"); 06087 } 06088 06089 // find and return inv s.t. m * inv = inv * m = I = identity 06090 // (m must be square) 06091 template<class T> 06092 TMat<T> inverse(TMat<T>& m) 06093 { 06094 TMat<T> inv(m.length(),m.length()); 06095 inverse(m,inv); 06096 return inv; 06097 } 06098 06099 // find inv s.t. m * inv = inv * m = I = identity 06100 // (m must be square) 06101 template<class T> 06102 void inverse(const TMat<T>& m, TMat<T>& inv) 06103 { 06104 int n=m.length(); 06105 if (m.width()!=n) 06106 PLERROR("inverse(TMat<T>,TMat<T>): argument(%d,%d) must be square matrix", 06107 m.width(), n); 06108 inv.resize(n,n); 06109 if (n==1) 06110 inv.data()[0]=1.0/m.data()[0]; 06111 else 06112 PLERROR("matrix inverse not implemented yet"); 06113 } 06114 06115 // for square positive definite symmetric matrices A, 06116 // find X(n,m) s.t. A(n,n) X(n,m) = B(n,m). 06117 // This is obtained by doing a Cholesky decomposition 06118 // A = L L', with L lower diagonal, thus to solve 06119 // L L' X = B. 06120 // We use the CholeskySolve function which solves for x_i in L L' x_i = b_i 06121 // (on the columns x_i and b_i of X and B respectively). 06122 // Optionally provide pointers to the temporary matrix L(n,n) and vector y(n) 06123 // to avoid memory allocations. 06124 template<class T> 06125 void solveLinearSystemByCholesky(const TMat<T>& A, const TMat<T>& B, TMat<T>& X, TMat<T>* pL=0, TVec<T>* py=0) 06126 { 06127 int n=A.length(); 06128 int m=X.width(); 06129 if (X.length()!=n || A.width()!=n || B.length()!=n || B.width()!=m) 06130 PLERROR("solveLinearSystemByCholesky: A(%d,%d) * X(%d,%d) == B(%d,%d), incompatible", 06131 n,A.width(),X.length(),m,B.length(),B.width()); 06132 TMat<T>* L; 06133 TVec<T>* y; 06134 if (pL) L=pL; else L = new TMat<T>(n,n); 06135 if (py) y=py; else y = new TVec<T>(n); 06136 choleskyDecomposition(A,*L); 06137 choleskySolve(*L,B,X,*y); 06138 if (!pL) delete L; 06139 if (!py) delete y; 06140 } 06141 06142 // for square positive definite symmetric matrices A, 06143 // find X(n,m) s.t. X(n,m) A(m,m) = B(n,m). 06144 // This is obtained by doing a Cholesky decomposition 06145 // A = L L', with L lower diagonal, thus to solve 06146 // X L L' = B. 06147 // We use the CholeskySolve function which solves for x_i in L L' x_i = b_i: 06148 // L L' X' = B' 06149 // is solved on the rows of X (x_i) and the columns of B (b_i). 06150 // Optionally provide pointers to the temporary matrices L and y 06151 // to avoid memory allocations. 06152 template<class T> 06153 void solveTransposeLinearSystemByCholesky(const TMat<T>& A, const TMat<T>& B, TMat<T>& X,TMat<T>* pL=0, TVec<T>* py=0) 06154 { 06155 int n=X.length(); 06156 int m=X.width(); 06157 if (A.length()!=m || A.width()!=m || B.length()!=n || B.width()!=m) 06158 PLERROR("solveTransposeLinearSystemByCholesky: X(%d,%d) * A(%d,%d) == B(%d,%d), incompatible", 06159 n,m,A.length(),A.width(),B.length(),B.width()); 06160 TMat<T>* L; 06161 TVec<T>* y; 06162 if (pL) L=pL; else L = new TMat<T>(m,m); 06163 if (py) y=py; else y = new TVec<T>(m); 06164 choleskyDecomposition(A,*L); 06165 for (int i=0;i<n;i++) 06166 choleskySolve(*L,B(i),X(i),*y); 06167 if (!pL) delete L; 06168 if (!py) delete y; 06169 } 06170 06171 /* Perform a Cholesky decomposition of nxn symmetric positive definite 06172 matrix A, i.e., decompose it into 06173 A = L L' 06174 where L is a lower diagonal matrix (with zeros above the diagonal). 06175 L be used to solve a linear system A x = b, i.e., LL'x=b, with choleskySolve(L,b,x). 06176 See choleskySolve(TMat<T>*,TVec<T>*) for an example of use. 06177 06178 From the above equation, one obtains 06179 06180 for i=0..n-1 06181 L[i][i] = sqrt(A[i][i] - sum_{k=0}^{i-1} L[i][k]^2) 06182 for j=i+1... n-1 06183 L[j][i] = (1/L[i][i]) ( A[i][j] - sum_{k=0}^{i-1} L[i][k] L[j][k] ) 06184 06185 */ 06186 template<class T> 06187 void choleskyDecomposition(const TMat<T>& A, TMat<T>& L) 06188 { 06189 int n = A.length(); 06190 if (n!=A.width()) 06191 PLERROR("choleskyDecomposition: non-square matrix %dx%d\n",n,A.width()); 06192 L.resize(n,n); 06193 int i,j,k; 06194 T sum; 06195 bool restart=false; 06196 do 06197 { 06198 restart=false; 06199 for (i=0;i<n;i++) 06200 { 06201 const T* Ai = A[i]; 06202 T* Li = L[i]; 06203 T Lii=0; 06204 for (j=i;j<n;j++) 06205 { 06206 T* Lj = L[j]; 06207 for (sum=Ai[j],k=i-1;k>=0;k--) sum -= Li[k] * Lj[k]; 06208 if (i==j) 06209 { 06210 if (sum <= 0.0) 06211 { 06212 T eps = -1.1*sum; 06213 if (fast_exact_is_equal(sum,0.0)) eps=1e-8; 06214 PLWARNING("Cholesky decomposition would fail: add %g to diagonal",eps); 06215 // saveAscii("A.amat",A); 06216 T* Aii=A.data(); 06217 int addm=A.mod()+1; 06218 for (int ii=0;ii<n;ii++,Aii+=addm) *Aii += eps; 06219 restart=true; 06220 break; 06221 } 06222 Lii = sqrt(sum); 06223 } 06224 else Lj[i] = sum/Lii; 06225 } 06226 if (restart) break; 06227 Li[i] = Lii; 06228 } 06229 } 06230 while (restart); 06231 06232 } 06233 06234 /* Back-propagate through the call to choleskyDecomposition(A,L). 06235 The argument A holds the original symmetric positive definite 06236 matrix while is the lower diagonal matrix with L L' = A. 06237 Given the derivative of C wrt L, fill the derivative 06238 of C wrt A. dC_dA must have been cleared beforehand. 06239 We are given A, L, dC_dL, and write into dC_dA. 06240 Note that dC_dL is modified after the call 06241 because of the internal dependencies between the L's. 06242 06243 for i=n-1..0 06244 for j=n-1..i+1 06245 dC_dL[i][i] -= dC_dL[j][i] L[j][i] / L[i][i] 06246 dC_dA[i][j] += dC_dL[j][i] / L[i][i] 06247 for k=0..i-1 06248 dC_dL[i][k] -= dC_dL[j][i] L[j][k] / L[i][i] 06249 dC_dL[j][k] -= dC_dL[j][i] L[i][k] / L[i][i] 06250 dC_dA[i][i] += 0.5 * dC_dL[i][i] / L[i][i] 06251 for k=0..i-1 06252 dC_dL[i][k] -= dC_dL[i][i] L[i][k] / L[i][i] 06253 06254 */ 06255 template<class T> 06256 void bpropCholeskyDecomposition(const TMat<T>& A, const TMat<T>& L, 06257 TMat<T>& dC_dA, TMat<T>& dC_dL) 06258 { 06259 int n = A.length(); 06260 if (dC_dA) 06261 dC_dA.resize(n,n); 06262 int i,j,k; 06263 for (i=n-1;i>=0;i--) 06264 { 06265 const T* Li = L[i]; 06266 T* dC_dLi = dC_dL[i]; 06267 T* dC_dAi = dC_dA[i]; 06268 T invLii = 1.0/Li[i]; 06269 for (j=n-1;j>i;j--) 06270 { 06271 const T* Lj = L[j]; 06272 T* dC_dLj = dC_dL[j]; 06273 T dC_dLji = dC_dLj[i]; 06274 dC_dLi[i] -= dC_dLji * Lj[i] * invLii; 06275 dC_dAi[j] += dC_dLji * invLii; 06276 for (k=0;k<i;k++) 06277 { 06278 dC_dLi[k] -= dC_dLji * Lj[k] * invLii; 06279 dC_dLj[k] -= dC_dLji * Li[k] * invLii; 06280 } 06281 } 06282 T dC_dLii = dC_dLi[i]; 06283 dC_dAi[i] += 0.5 * dC_dLii * invLii; 06284 for (k=0;k<i;k++) 06285 dC_dLi[k] -= dC_dLii * Li[k] * invLii; 06286 } 06287 } 06288 06289 // Given L lower-diagonal, solve L y = b 06290 template<class T> 06291 void choleskyLeftSolve(const TMat<T>& L, const TVec<T>& b, const TVec<T>& y) 06292 { 06293 int i,k; 06294 T sum; 06295 int n = L.length(); 06296 #ifdef BOUNDCHECK 06297 if (L.width()!=n) 06298 PLERROR("choleskySolve: matrix L (%d x %d) is not square!", 06299 n, L.width()); 06300 if (b.length()!=n || y.length()!=n) 06301 PLERROR("choleskySolve: RHS vector b(%d) or unknown y(%d) incompatible with L(%d,%d)", 06302 b.length(),y.length(),n,n); 06303 #endif 06304 06305 if (n == 0) 06306 // Empty matrix, there is nothing that needs being solved. 06307 return; 06308 06309 T* bp = b.data(); 06310 T* yp = y.data(); 06311 06312 // solve L y = b (in variable x if y=0): 06313 // for i=0..n-1 06314 // y[i] = (b[i] - sum_{k<i} L[i][k] y[k])/L[i][i] 06315 for (i=0;i<n;i++) 06316 { 06317 const T* Li = L[i]; 06318 for (sum=bp[i],k=i-1;k>=0;k--) sum -= Li[k] * yp[k]; 06319 if (Li[i]==0) 06320 PLERROR("choleskyLeftSolve: found zero entry in diagonal of L (%d)",i); 06321 yp[i] = sum / Li[i]; 06322 } 06323 } 06324 06325 // Given L lower-diagonal, solve L' x = y 06326 template<class T> 06327 void choleskyRightSolve(const TMat<T>& L, TVec<T>& y, TVec<T>& x) 06328 { 06329 int i,k; 06330 T sum; 06331 int n = L.length(); 06332 #ifdef BOUNDCHECK 06333 if (L.width()!=n) 06334 PLERROR("choleskySolve: matrix L (%d x %d) is not square!", 06335 n, L.width()); 06336 if (x.length()!=n || y.length()!=n) 06337 PLERROR("choleskySolve: RHS vector y(%d) or unknown x(%d) incompatible with L(%d,%d)", 06338 y.length(),x.length(),n,n); 06339 #endif 06340 06341 if (n == 0) 06342 // Empty matrix, there is nothing that needs being solved. 06343 return; 06344 06345 T* xp = x.data(); 06346 T* yp = y.data(); 06347 06348 // for i=n-1..0 06349 // x[i] = (y[i] - sum_{k>i} L[k][i] x[k])/L[i][i] 06350 for (i=n-1;i>=0;i--) 06351 { 06352 for (sum=yp[i],k=i+1;k<n;k++) sum -= L[k][i] * xp[k]; 06353 xp[i] = sum / L[i][i]; 06354 } 06355 } 06356 06357 /* Solve the linear system A x = L L' x = b using a Cholesky decomposition 06358 of A into L L' performed with a prior call to choleskyDecomposition(A,L) 06359 (which on return has the matrix L, that is lower diagonal, and A = L L'). 06360 The solution of the linear system L L' x = b will be in x. 06361 See choleskySolve(TMat<T>*,TVec<T>*) for an example of use. 06362 The algorithm is first to solve L y = b, and then L' x = y. 06363 The argument y is optional and can be used to hold the intermediate 06364 solution to L y = b. 06365 06366 The solution to L L' x = b is obtained as follows: 06367 06368 * Solve L y = b by iterating once through the rows of L 06369 (store result in x): 06370 y[i] = (b[i] - sum_{k<i} L[i][k] y[k])/L[i][i] 06371 06372 * Solve L' x = y by iterating once (backwards) through the rows of L. 06373 x[i] = (y[i] - sum_{k>i} L[k][i] x[k])/L[i][i] 06374 06375 */ 06376 template<class T> 06377 void choleskySolve(const TMat<T>& L, TVec<T> b, TVec<T> x, TVec<T>& y) 06378 { 06379 // solve L y = b 06380 choleskyLeftSolve(L,b,y); 06381 // solve L' x = y 06382 choleskyRightSolve(L,y,x); 06383 } 06384 06385 // same as the previous choleskySolve but do it m times on the columns 06386 // of nxm matrices X and B. 06387 template<class T> 06388 void choleskySolve(const TMat<T>& L, const TMat<T>& B, TMat<T>& X, TVec<T>& y) 06389 { 06390 int i,k; 06391 T sum; 06392 int n = L.length(); 06393 int m = X.width(); 06394 if (L.width()!=n) 06395 PLERROR("choleskySolve: matrix L (%d x %d) is not square!", 06396 n, L.width()); 06397 if (B.length()!=n || B.width() !=m) 06398 PLERROR("choleskySolve: RHS matrix B(%d,%d) instead of (%d,%d) like X", 06399 B.length(),B.width(), n, m); 06400 if (X.length()!=n) 06401 PLERROR("choleskySolve: X(%d,%d) not compatible with L(%d,%d)", 06402 X.length(),m,n,n); 06403 if (y.length()!=n) 06404 PLERROR("choleskySolve: y(%d) not compatible with L(%d,%d)", 06405 y.length(),n,n); 06406 int bmod = B.mod(); 06407 int xmod = X.mod(); 06408 // loop over columns b and x of B and X 06409 for (int j=0;j<m;j++) 06410 { 06411 T* bp = B.data()+j; 06412 T* yp = y.data(); 06413 // solve L y = b (in variable x if y=0): 06414 // for i=0..n-1 06415 // y[i] = (b[i] - sum_{k<i} L[i][k] y[k])/L[i][i] 06416 for (i=0;i<n;i++,bp+=bmod) 06417 { 06418 const T* Li = L[i]; 06419 for (sum = *bp,k=i-1;k>=0;k--) sum -= Li[k] * yp[k]; 06420 yp[i] = sum / Li[i]; 06421 } 06422 // solve L' x = y 06423 // for i=n-1..0 06424 // x[i] = (y[i] - sum_{k>i} L[k][i] x[k])/L[i][i] 06425 for (i=n-1;i>=0;i--) 06426 { 06427 sum=yp[i]; 06428 if (i+1<n) 06429 { 06430 T* xp = &X(i+1,j); 06431 for (k=i+1;k<n;k++,xp+=xmod) sum -= L[k][i] * *xp; 06432 } 06433 X(i,j) = sum / L[i][i]; 06434 } 06435 } 06436 } 06437 06438 /* 06439 Back-propagate through the CholeskySolve(L,b,x,y) operation 06440 (the optional argument y of this call must have been provided). 06441 06442 dC_dL and dC_db must have been cleared beforehand. 06443 dC_dx will be modified (because of the dependencies between the x's. 06444 06445 (1) back-prop through step L' x = y: 06446 for i=0..n-1 06447 dC_dy[i] = dC_dx[i] / L[i][i] 06448 dC_dL[i][i] -= dC_dx[i] x[i] / L[i][i] 06449 for k=i+1..n 06450 dC_dx[k] -= dC_dx[i] L[k][i] / L[i][i] 06451 dC_dL[k][i] -= dC_dx[i] x[k] / L[i][i] 06452 06453 (2) back-prop through step L y = b: 06454 for i=n-1..0 06455 dC_db[i] = dC_dy[i] / L[i][i] 06456 dC_dL[i][i] -= dC_dy[i] y[i] / L[i][i] 06457 for k=0..i-1 06458 dC_dy[k] -= dC_dy[i] L[i][k] / L[i][i] 06459 dC_dL[i][k] -= dC_dy[i] * y[k] / L[i][i] 06460 */ 06461 template<class T> 06462 void bpropCholeskySolve(const TMat<T>& L, const TVec<T>& x, const TVec<T>& y, 06463 TMat<T>& dC_dL, TVec<T>& dC_db, TVec<T>& dC_dx) 06464 { 06465 int n = L.length(); 06466 int i,k; 06467 TVec<T> dC_dy(n); 06468 const T *xp = x.data(); 06469 const T *yp = y.data(); 06470 T *dC_dbp = dC_db.data(); 06471 T *dC_dxp = dC_dx.data(); 06472 T* dC_dyp = dC_dy.data(); 06473 06474 // (1) back-prop through step L' x = y: 06475 for (i=0;i<n;i++) 06476 { 06477 const T* Li = L[i]; 06478 T invLii = 1.0 / Li[i]; 06479 dC_dyp[i] = dC_dxp[i] * invLii; 06480 T dC_dxi = dC_dxp[i]; 06481 dC_dL[i][i] -= dC_dxp[i] * xp[i] * invLii; 06482 for (k=i+1;k<n;k++) 06483 { 06484 dC_dxp[k] -= dC_dxi * L[k][i] * invLii; 06485 dC_dL[k][i] -= dC_dxi * xp[k] * invLii; 06486 } 06487 } 06488 06489 // (2) back-prop through step L y = b: 06490 for (i=n-1;i>=0;i--) 06491 { 06492 const T* Li = L[i]; 06493 T* dC_dLi = dC_dL[i]; 06494 T invLii = 1.0 / Li[i]; 06495 T dC_dyi = dC_dyp[i]; 06496 T dC_dyi_over_Lii = dC_dyi * invLii; 06497 dC_dbp[i] += dC_dyi_over_Lii; 06498 dC_dLi[i] -= dC_dyi_over_Lii * yp[i]; 06499 for (k=0;k<i;k++) 06500 { 06501 dC_dyp[k] -= dC_dyi_over_Lii * Li[k]; 06502 dC_dLi[k] -= dC_dyi_over_Lii * yp[k]; 06503 } 06504 }; 06505 } 06506 06507 /* Use Cholesky decomposition to invert symmetric 06508 positive definite matrix A. 06509 Also returns the log of the determinant of A 06510 06511 We have L L' = A, and we want to solve L L' Ainv = I. 06512 06513 1) solve L Linv = I, i.e., invert L 06514 06515 for j=0..n-1 06516 Linv[j][j] = 1 / L[j][j] 06517 for i=j+1..n-1 06518 Linv[i][j] = - sum_{j<=k<i} L[i][k] Linv[k][j] / L[i][i] 06519 and 0 elsewhere (Linv is lower diagonal) 06520 06521 2) solve L' Ainv = Linv 06522 06523 for j=0..n-1 06524 for i=n-1..0 06525 Ainv[i][j] = (Linv[i][j] - sum_{k>i} L[k][i] Ainv[k][j])/L[i][i] 06526 06527 */ 06528 template<class T> 06529 real choleskyInvert(const TMat<T>& A, TMat<T>& Ainv) 06530 { 06531 int n= A.length(); 06532 TMat<T> L(n,n); 06533 Ainv.resize(n,n); 06534 06535 choleskyDecomposition(A,L); 06536 // now L L' = A 06537 06538 real logdet = pl_log(fabs(L(0,0))); 06539 for(int i=1; i<n; i++) 06540 logdet += pl_log(fabs(L(i,i))); 06541 logdet *= 2; 06542 06543 // Compute Linv and put its transpose above L's diagonal. 06544 // and put Linv[i][i] = 1 / L[i][i] in L's diagonal. 06545 int i,j; 06546 T *Lii = L.data(); 06547 for (i=0;i<n;i++,Lii+=1+n) 06548 *Lii = 1.0 / *Lii; 06549 06550 for (j=0;j<n;j++) 06551 { 06552 T *Linv_xj = L[j]; // Linv' in L's upper triangle 06553 for (i=j+1;i<n;i++) 06554 { 06555 T sum=0.0; 06556 T* Li = L[i]; 06557 int k; 06558 for (k=j;k<i;k++) sum -= Li[k] * Linv_xj[k]; 06559 Linv_xj[i] = sum * Li[i]; // * not / because inverse already done above 06560 } 06561 } 06562 // recall: now Linv above and on diagonal of L, L below it, 06563 06564 // compute A's inverse 06565 for (j=0;j<n;j++) 06566 { 06567 T* Linv_xj = L[j]; 06568 for (i=n-1;i>=j;i--) 06569 { 06570 T sum = Linv_xj[i]; // this is Linv[i][j] 06571 int k; 06572 for (k=i+1;k<n;k++) 06573 sum -= L[k][i] * Ainv[k][j]; 06574 Ainv[i][j] = sum * L[i][i]; 06575 } 06576 for (i=j-1;i>=0;i--) // symmetric part 06577 Ainv[i][j] = Ainv[j][i]; 06578 }; 06579 06580 return logdet; 06581 } 06582 06583 /* Solve a linear system of equations A x = b, when A is 06584 symmetric positive definite. Return x. */ 06585 template<class T> 06586 TVec<T> choleskySolve(const TMat<T>& A, const TVec<T>& b) 06587 { 06588 int n = A.length(); 06589 TMat<T> L(n,n); 06590 TVec<T> x(n); 06591 choleskyDecomposition(A,L); 06592 choleskySolve(L,b,x); 06593 return x; 06594 } 06595 06596 /* return inverse of positive definite matrix A 06597 using Cholesky decomposition. No side-effect on A. */ 06598 template<class T> 06599 TMat<T> choleskyInvert(const TMat<T>& A) 06600 { 06601 int n=A.length(); 06602 TMat<T> Ainv(n,n); 06603 choleskyInvert(A,Ainv); 06604 return Ainv; 06605 } 06606 06607 template<class T> 06608 void LU_decomposition(TMat<T>& A, TVec<T>& Trow, int& detsign, TVec<T>* p=0) 06609 { 06610 int n=A.length(); 06611 if (n!=A.width()) 06612 PLERROR("LU_decomposition: matrix A(%d,%d) should be square", n,A.width()); 06613 TVec<T>* pivot = (p==0)?new TVec<T>(n):p; 06614 T* pv = pivot->data(); 06615 detsign = 1; 06616 for (int i=0;i<n;i++) 06617 { 06618 T max_abs = maxabs(A(i)); 06619 if (max_abs==0) 06620 PLERROR("LU_decomposition: %d-th row has only zeros",i); 06621 pv[i] = 1.0 / max_abs; 06622 } 06623 int mod = A.mod(); 06624 for (int j=0;j<n;j++) 06625 { 06626 for (int i=0;i<j;i++) 06627 { 06628 T* Ai = A[i]; 06629 T* Akj = A.data()+j; 06630 T Uij = Ai[j]; 06631 for (int k=0;k<i;k++,Akj+=mod) 06632 Uij -= Ai[k] * *Akj; 06633 Ai[j] = Uij; 06634 } 06635 T max_abs = 0; 06636 int maxi = 0; 06637 for (int i=j;i<n;i++) 06638 { 06639 T* Ai = A[i]; 06640 T* Akj = A.data()+j; 06641 T Lij = Ai[j]; 06642 for (int k=0;k<j;k++,Akj+=mod) 06643 Lij -= Ai[k] * *Akj; 06644 Ai[j] = Lij; 06645 T piv = fabs(Lij) * pv[i]; 06646 if (piv >= max_abs) 06647 { 06648 maxi = i; 06649 max_abs = piv; 06650 } 06651 } 06652 if (j!=maxi) 06653 // swap row j and row maxi 06654 { 06655 A.swapRows(j,maxi); 06656 pv[maxi]=pv[j]; 06657 detsign = -detsign; 06658 } 06659 Trow[j] = maxi; 06660 T& Ajj = A(j,j); 06661 if (Ajj==0) Ajj=1e-20; // some regularization of singular matrices 06662 if (j<n-1) 06663 { 06664 T denom = 1.0/Ajj; 06665 T* Aij = &A(j+1,j); 06666 for (int i=j+1;i<n;i++, Aij+=mod) 06667 *Aij *= denom; 06668 } 06669 } 06670 if (p == 0) delete pivot; 06671 } 06672 06675 template<class T> 06676 T det(const TMat<T>& A, bool log_det = false) 06677 { 06678 // Work storage. 06679 static TMat<T> LU; 06680 static TVec<T> Trow, p; 06681 06682 int n = A.length(); 06683 if (n!=A.width()) 06684 PLERROR("det(const TMat<T>& A): A(%d,%d) is not square!",n,A.width()); 06685 for (int i=0;i<n;i++) 06686 { 06687 T max_abs = maxabs(A(i)); 06688 if (max_abs==0) 06689 return 0.0; 06690 } 06691 LU.resize(A.length(), A.width()); 06692 LU << A; 06693 Trow.resize(n); 06694 p.resize(n); 06695 int detsign; 06696 LU_decomposition(LU, Trow, detsign, &p); 06697 return det(LU, detsign, log_det); 06698 } 06699 06703 template<class T> 06704 T det(const TMat<T>& LU, int detsign, bool log_det = false) 06705 { 06706 T determinant = detsign; 06707 bool minus = false; 06708 if (log_det) { 06709 if (detsign < 0) { 06710 minus = !minus; 06711 detsign = - detsign; 06712 } 06713 determinant = pl_log(double(detsign)); 06714 } 06715 int mod = LU.mod(); 06716 int n = LU.width(); 06717 if (n!=LU.width()) 06718 PLERROR("det(const TMat<T>& LU, int detsign): LU(%d,%d) is not square!",n,LU.width()); 06719 T* LUii = LU.data(); 06720 if (log_det) { 06721 for (int i=0;i<n;i++, LUii+=1+mod) { 06722 real LUii_ = *LUii; 06723 if (LUii_ < 0) { 06724 minus = !minus; 06725 LUii_ = - LUii_; 06726 } 06727 determinant += pl_log(LUii_); 06728 } 06729 } else { 06730 for (int i=0;i<n;i++, LUii+=1+mod) 06731 determinant *= *LUii; 06732 } 06733 if (log_det && minus) 06734 // The determinant is negative: its log should be NaN. 06735 determinant = MISSING_VALUE; 06736 return determinant; 06737 } 06738 06739 // dest[i,j] = 1 if src[i,j]==v, 0 otherwise 06740 template<class T> 06741 void equals(const TMat<T>& src, T v, TMat<T>& dest) 06742 { 06743 int l=src.length(); 06744 int w=src.width(); 06745 #ifdef BOUNDCHECK 06746 if (l!=dest.length() || w!=dest.width()) 06747 PLERROR("equals(TMat<T>(%d,%d),T,TMat<T>(%d,%d)) args of unequal dimensions", 06748 src.length(),src.width(),dest.length(),dest.width()); 06749 #endif 06750 for (int i=0;i<l;i++) 06751 { 06752 const T* s=src[i]; 06753 T* d=dest[i]; 06754 for (int j=0;j<w;j++) 06755 if (s[j]==v) d[j]=1.0; else d[j]=0.0; 06756 } 06757 } 06758 06759 // dest[i,j] = src[i,j] 06760 template<class T> 06761 void transpose(const TMat<T> src, TMat<T> dest) 06762 { 06763 int l=src.length(); 06764 int w=src.width(); 06765 #ifdef BOUNDCHECK 06766 if (w!=dest.length() || l!=dest.width()) 06767 PLERROR("transpose(TMat<T>(%d,%d),T,TMat<T>(%d,%d)) args of unequal dimensions", 06768 src.length(),src.width(),dest.length(),dest.width()); 06769 #endif 06770 int dmod=dest.mod(); 06771 for (int i=0;i<l;i++) 06772 { 06773 const T* si=src[i]; 06774 T* dji= &dest(0,i); 06775 for (int j=0;j<w;j++,dji+=dmod) 06776 *dji = si[j]; 06777 } 06778 } 06779 06780 // res[i,j] = src[i,j] 06781 template<class T> 06782 TMat<T> transpose(const TMat<T>& src) 06783 { 06784 TMat<T> res(src.width(),src.length()); 06785 transpose(src,res); 06786 return res; 06787 } 06788 06790 template<class T, class U> 06791 void apply(U (*func)(T), const TMat<T>& source, TMat<U>& destination) 06792 { 06793 int l=source.length(); 06794 int w=source.width(); 06795 if (l!=destination.length() || w!=destination.width()) 06796 PLERROR("apply: source(%d,%d) TMat<T> and destination(%d,%d) TMat<U> must have same length and width", 06797 l,w,destination.length(),destination.width()); 06798 for(int i=0; i<l; i++) { 06799 for(int j=0; j<w; j++) 06800 destination(i,j)=func(source(i,j)); 06801 } 06802 } 06803 06806 template<class T, class U> 06807 void apply(const TMat<T>& source, TMat<U>& destination, U (*func)(T)) 06808 { 06809 apply(func, source, destination); 06810 } 06811 06812 // Apply a vector operation to each row of matrices, result in rows of a matrix 06813 template<class T> 06814 void apply(T (*func)(const TVec<T>&), const TMat<T>& m, TMat<T>& dest) 06815 { 06816 if (dest.length()!=m.length()) 06817 PLERROR("apply: m.length_=%d, dest.length_=%d", 06818 m.length(),dest.length()); 06819 int l=m.length(); 06820 for (int i=0;i<l;i++) 06821 dest(i,0)=func(m(i)); 06822 } 06823 06824 template<class T> 06825 void apply(T (*func)(const TVec<T>&,const TVec<T>&), const TMat<T>& m1, const TMat<T>& m2, 06826 TMat<T>& dest) 06827 { 06828 if (dest.length()!=m1.length() || m1.length()!=m2.length()) 06829 PLERROR("apply: m1.length_=%d, m2.length_=%d, dest.length_=%d", 06830 m1.length(),m2.length(),dest.length()); 06831 for (int i=0;i<m1.length();i++) 06832 dest(i,0)=func(m1(i),m2(i)); 06833 } 06834 06835 // Perform a traditional linear regression (but with weight decay), 06836 // without bias term. i.e. find weights such that: 06837 // 06838 // norm(weights*inputs - outputs) + weight_decay*norm(weights) 06839 // 06840 // is minimized, 06841 // 06842 // This is achieved by solving the following linear system: 06843 // 06844 // (X' X + weight_decay I) * weights = X' outputs 06845 06846 template<class T> 06847 void linearRegressionNoBias(TMat<T> inputs, TMat<T> outputs, T weight_decay, 06848 TMat<T> weights) 06849 { 06850 int inputsize = inputs.width(); 06851 int outputsize = outputs.width(); 06852 int l = inputs.length(); 06853 if(outputs.length()!=l) 06854 PLERROR("In linearRegressionNoBias: inputs and outputs should have the same length"); 06855 if(weights.length()!=inputsize || weights.width()!=outputsize) 06856 PLERROR("In linearRegressionNoBias: weights should be a (inputsize x outputsize) matrix (%d x %d)",inputsize,outputsize); 06857 static TMat<T> XtX; 06858 XtX.resize(inputsize,inputsize); 06859 transposeProduct(XtX, inputs,inputs); 06860 static TMat<T> XtY; 06861 XtY.resize(inputsize,outputsize); 06862 transposeProduct(XtY, inputs,outputs); 06863 for(int i=0; i<inputsize; i++) 06864 XtX(i,i) += weight_decay; 06865 solveLinearSystemByCholesky(XtX,XtY,weights); 06866 } 06867 06868 06869 // Perform a traditional linear regression (but with weight decay), 06870 // i.e. find bias and weights such that 06871 // 06872 // norm(bias + weights*inputs - outputs) + weight_decay*norm(weights) 06873 // 06874 // is minimized, where theta'=(biases;weights') {biases in first row} 06875 // 06876 // This is achieved by solving the following linear system: 06877 // 06878 // (X' X + weight_decay I) * theta = X' outputs 06879 // 06880 // where X = augmented inputs, i.e. X(t) = (1,inputs(t)) 06881 // 06882 template<class T> 06883 void linearRegression(TMat<T> inputs, TMat<T> outputs, T weight_decay, 06884 TMat<T> theta_t) 06885 { 06886 int l = inputs.length(); 06887 int n_inputs = inputs.width(); 06888 int n_outputs = outputs.width(); 06889 if (outputs.length()!=l) 06890 PLERROR("linearRegression: inputs.length_=%d while outputs.length_=%d", 06891 l,outputs.length()); 06892 if (theta_t.length()!=n_inputs+1 || theta_t.width()!=n_outputs) 06893 PLERROR("linearRegression: theta_t(%d,%d) should be (n_inputs(%d)+1)xn_outputs(%d)", 06894 theta_t.length(),theta_t.width(),n_inputs,n_outputs); 06895 06896 int n=n_inputs+1; 06897 06898 static TMat<T> XtX; 06899 XtX.resize(n,n); 06900 XtX.clear(); 06901 static TMat<T> XtY; 06902 XtY.resize(n,n_outputs); 06903 XtY.clear(); 06904 // compute X' X and X'Y: 06905 // XtX(i,j) = sum_t X[t,i]*X[t,j] (with X[t,0]=1, X[t,i+1]=inputs[t,i]) 06906 // YtY(i,j) = sum_t X[t,i]*Y[t,j] 06907 // 06908 int xmod=inputs.mod(); 06909 int ymod=outputs.mod(); 06910 T *xt = inputs.data(); 06911 T *yt = outputs.data(); 06912 XtX(0,0) = l; // we know the answer ahead of time for element (0,0) 06913 for (int t=0;t<l;t++,xt+=xmod,yt+=ymod) 06914 { 06915 T* xx0 = XtX.data(); 06916 T* xy0 = XtY.data(); 06917 for (int j=0;j<n_outputs;j++) 06918 xy0[j] += yt[j]; 06919 T *xxi = xx0+n; // start the inner matrix at (1,0) 06920 T *xyi = xy0+n_outputs; // start xy at (1,0) 06921 for (int i=0;i<n_inputs;i++,xxi+=n,xyi+=n_outputs) 06922 { 06923 T xti = xt[i]; 06924 xxi[0]+=xti; 06925 T *xxip=xxi+1; 06926 for (int j=0;j<i;j++) 06927 xxip[j] += xti*xt[j]; 06928 xxip[i]+=xti*xti; 06929 for (int j=0;j<n_outputs;j++) 06930 xyi[j] += xti * yt[j]; 06931 } 06932 } 06933 // now do the symmetric part of XtX 06934 T* xx = XtX.data(); 06935 T* xxi = xx+n; 06936 for (int i=1;i<n;i++,xxi+=n) 06937 { 06938 T *xx_i=xx+i; 06939 for (int j=0;j<i;j++,xx_i+=n) 06940 *xx_i = xxi[j]; 06941 } 06942 06943 // add weight_decay on the diagonal of XX' (except for the bias) 06944 T* xxii = &XtX(1,1); 06945 for (int i=0;i<n_inputs;i++,xxii+=1+n) 06946 *xxii += weight_decay; 06947 06948 // now solve by Cholesky decomposition 06949 solveLinearSystemByCholesky(XtX,XtY,theta_t); 06950 } 06951 06952 // Compute a linear fitting of 2 dimensional data resulting 06953 // in parameters m et b for y = mx + b 06954 // 1 1 06955 // Cost function used: C = - Sum[t] { (m * x_t + b - y_t)^2 } + - weight_decay * m^2 06956 // 2 2 06957 06958 template<class T> 06959 void linearRegression(TVec<T> inputs, TVec<T> outputs, T weight_decay, TVec<T> theta_t) 06960 { 06961 int npts = inputs.length(); 06962 06963 if (outputs.length()!=npts) 06964 PLERROR("linearRegression: inputs.length_=%d while outputs.length_=%d", 06965 inputs.length(),outputs.length()); 06966 if (theta_t.length()!=2) 06967 PLERROR("linearRegression: theta_t(%d) should be 2", theta_t.length()); 06968 06969 T sum_x = 0, sum_y = 0, sum_xy = 0, sum_x2 = 0, sum2_x = 0, sum2_y = 0; 06970 06971 for (int i = 0; i < npts; ++i) { 06972 sum_x += inputs[i]; 06973 sum_y += outputs[i]; 06974 sum_xy += inputs[i] * outputs[i]; 06975 sum_x2 += inputs[i] * inputs[i]; 06976 } 06977 sum2_x = sum_x * sum_x; 06978 sum2_y = sum_y * sum_y; 06979 06980 // m 06981 theta_t[1] = (sum_xy - (sum_x * sum_y) / npts) / (sum_x2 + weight_decay - sum2_x / npts); 06982 // b 06983 theta_t[0] = (sum_y - theta_t[1] * sum_x) / npts; 06984 } 06985 06986 06987 template<class T> 06988 TMat<T> smooth(TMat<T> data, int windowsize) 06989 { 06990 TVec<T> sumvec(data.width()); 06991 TMat<T> result(data.length(), data.width()); 06992 int currentwindowsize = windowsize/2; 06993 for(int k=0; k<currentwindowsize; k++) 06994 sumvec += data(k); 06995 result(0) << sumvec; 06996 //result(0) /= (T)currentwindowsize; 06997 TVec<T> res0 = result(0); 06998 res0 /= (T)currentwindowsize; 06999 result(0) << res0; 07000 07001 for(int i=0; i<data.length(); i++) 07002 { 07003 int lowi = i-(windowsize-1)/2; // lowest index of window rows (inclusive) 07004 int highi = i+windowsize/2; // highest index of window rows (inclusive) 07005 if(lowi-1>=0) // remove row lowi-1 if it exists 07006 { 07007 sumvec -= data(lowi-1); 07008 currentwindowsize--; 07009 } 07010 if(highi<data.length()) // add row highi if it exists 07011 { 07012 sumvec += data(highi); 07013 currentwindowsize++; 07014 } 07015 result(i) << sumvec; 07016 //result(i) /= (T)currentwindowsize; 07017 TVec<T> resi = result(i); 07018 resi /= (T)currentwindowsize; 07019 result(i) << resi; 07020 } 07021 07022 07023 return result; 07024 } 07025 07026 07027 template<class T> 07028 TMat<T> square(const TMat<T>& m) 07029 { 07030 TMat<T> res(m.length(), m.width()); 07031 int w=m.width(); 07032 for(int i=0; i<m.length(); i++) 07033 for(int j=0; j<w; j++) 07034 res(i,j) = square(m(i,j)); 07035 return res; 07036 } 07037 07038 template<class T> 07039 TMat<T> sqrt(const TMat<T>& m) 07040 { 07041 TMat<T> res(m.length(), m.width()); 07042 int w=m.width(); 07043 for(int i=0; i<m.length(); i++) 07044 for(int j=0; j<w; j++) 07045 res(i,j) = sqrt(m(i,j)); 07046 return res; 07047 } 07048 07049 template<class T> 07050 inline void affineMatrixInitialize(TMat<T> W, bool output_on_columns=true, real scale=1.0) 07051 { 07052 int n_inputs = output_on_columns?W.width():W.length(); 07053 real delta = scale/n_inputs; 07054 fill_random_uniform(W,-delta,delta); 07055 W(0).clear(); 07056 } 07057 07058 template<class T> 07059 TMat<T> grep(TMat<T> data, int col, TVec<T> values, bool exclude=false) 07060 { 07061 TMat<T> result(data.length(),data.width()); 07062 int length=0; 07063 07064 for(int i=0; i<data.length(); i++) 07065 { 07066 bool contains = values.contains(data(i,col)); 07067 if( (!exclude && contains) || (exclude && !contains) ) 07068 result(length++) << data(i); 07069 } 07070 result.resize(length,result.width()); 07071 result.compact(); // use less memory 07072 return result; 07073 } 07074 07075 07076 template<class T> 07077 void convolve(TMat<T> m, TMat<T> mask, TMat<T> result) 07078 { 07079 if(result.length() != m.length()-mask.length()+1 || result.width() != m.width()-mask.width()+1) 07080 PLERROR("In convolve(TMat<T> m, TMat<T> mask, TMat<T> result), result does not have the appropriate dimensions"); 07081 T sum; 07082 for(int i=0; i<result.length(); i++) 07083 for(int j=0; j<result.width(); j++) 07084 { 07085 T* maskptr = mask.data(); 07086 T* mptr = m[i]+j; 07087 sum = 0.0; 07088 int w=mask.width(); 07089 07090 for(int l=0; l<mask.length(); l++, maskptr += mask.mod(), mptr += m.mod()) 07091 for(int c=0; c<w; c++) 07092 sum += maskptr[c] * mptr[c]; 07093 result(i,j) = sum; 07094 } 07095 } 07096 07097 template<class T> 07098 void subsample(TMat<T> m, int thesubsamplefactor, TMat<T> result) 07099 { 07100 T sum; 07101 int norm = thesubsamplefactor * thesubsamplefactor; 07102 for(int i=0; i<result.length(); i++) 07103 for(int j=0; j<result.width(); j++) 07104 { 07105 T* mptr = m[thesubsamplefactor*i]+thesubsamplefactor*j; 07106 sum = 0.0; 07107 for(int l=0; l<thesubsamplefactor; l++, mptr += m.mod()) 07108 for(int c=0; c<thesubsamplefactor; c++) 07109 sum += mptr[c]; 07110 result(i,j) = sum/norm; 07111 } 07112 } 07113 07114 07115 template<class T> 07116 void classification_confusion_matrix(TMat<T> outputs, TMat<T> target_classes, TMat<T> confusion_matrix) 07117 { 07118 int argmax, target; 07119 T v_max, tmp; 07120 07121 for (int i=0; i<outputs.length(); i++) { 07122 // Find argmax(outputs) 07123 v_max = outputs(i,0); 07124 argmax = 0; 07125 for (int j=1; j<outputs.width(); ++j) { 07126 tmp = outputs(i,j); 07127 if (tmp > v_max) { 07128 argmax = j; 07129 v_max = tmp; 07130 } 07131 } 07132 // Update confusion matrix 07133 target = (int) target_classes(i,0); 07134 confusion_matrix(argmax, target) ++; 07135 } 07136 } 07137 07151 template<class T> 07152 int GramSchmidtOrthogonalization(TMat<T> A, T tolerance=1e-6) 07153 { 07154 int n_basis = 0; 07155 for (int i=0;i<A.length();i++) 07156 { 07157 TVec<T> Ai=A(i); 07158 if (n_basis!=i) 07159 { 07160 TVec<T> Ab = A(n_basis); 07161 Ab << Ai; 07162 Ai=Ab; 07163 } 07164 if (i>0) 07165 projectOnOrthogonalSubspace(Ai, A.subMatRows(0,n_basis)); 07166 T normAi = norm(Ai); 07167 if (normAi>1e-6) 07168 { 07169 if (normAi!=1) 07170 Ai/=normAi; 07171 n_basis++; 07172 } 07173 // else ignore row i 07174 } 07175 return n_basis; 07176 } 07177 07179 07181 template<class T> 07182 inline TVec<T> product(const TMat<T>& m, const TVec<T>& v) 07183 { TVec<T> res(m.length()); product(res, m,v); return res; } 07184 07186 template<class T> 07187 inline TVec<T> transposeProduct(const TMat<T>& m, const TVec<T>& v) 07188 { TVec<T> res(m.width()); transposeProduct(res, m,v); return res; } 07189 07191 template<class T> 07192 inline TMat<T> product(const TMat<T>& m1, const TMat<T>& m2) 07193 { TMat<T> res(m1.length(),m2.width()); product(res, m1,m2); return res; } 07194 07196 template<class T> 07197 inline TMat<T> transposeProduct(const TMat<T>& m1, const TMat<T>& m2) 07198 { TMat<T> res(m1.width(),m2.width()); transposeProduct(res, m1,m2); return res; } 07199 07201 template<class T> 07202 inline TMat<T> productTranspose(const TMat<T>& m1, const TMat<T>& m2) 07203 { TMat<T> res(m1.length(),m2.length()); productTranspose(res, m1,m2); return res; } 07204 07206 template<class T> 07207 inline TMat<T> operator+(const TMat<T>& m, const TVec<T>& v) 07208 { TMat<T> res = m.copy(); res+=v; return res; } 07209 07211 template<class T> 07212 inline TMat<T> operator-(const TMat<T>& m, const TVec<T>& v) 07213 { TMat<T> res = m.copy(); res-=v; return res; } 07214 07216 template<class T> 07217 inline TMat<T> operator*(const TMat<T>& m, const TVec<T>& v) 07218 { TMat<T> res = m.copy(); res*=v; return res; } 07219 07221 template<class T> 07222 inline TMat<T> operator/(const TMat<T>& m, const TVec<T>& v) 07223 { TMat<T> res = m.copy(); res/=v; return res; } 07224 07226 template<class T> 07227 inline TMat<T> operator/(const TMat<T>& m1, const TMat<T>& m2) 07228 { TMat<T> res = m1.copy(); res/=m2; return res; } 07229 07230 template<class T> 07231 inline void choleskySolve(const TMat<T>& L, TVec<T> b, TVec<T> x) 07232 { TVec<T> y(b.size()); choleskySolve(L,b,x,y); } 07233 07235 template<class T> 07236 inline TMat<T> grep(TMat<T> data, int col, T value, bool exclude=false) 07237 { return grep(data,col,TVec<T>(1,value),exclude); } 07238 07239 template<class T> 07240 void addIfNonMissing(const TVec<T>& source, const TVec<int>& nnonmissing, TVec<T> destination) 07241 { 07242 #ifdef BOUNDCHECK 07243 if (source.length()!=nnonmissing.length() || source.length()!=destination.length()) 07244 PLERROR("addIfNonMissing: all arguments should have the same length, got %d,%d,%d\n", 07245 source.length(),nnonmissing.length(),destination.length()); 07246 #endif 07247 T* s=source.data(); 07248 T* d=destination.data(); 07249 int* n=nnonmissing.data(); 07250 int size=source.length(); 07251 for (int i=0;i<size;i++) 07252 if (finite(s[i])) 07253 { 07254 d[i] += s[i]; 07255 n[i]++; 07256 } 07257 } 07258 07259 template<class T> 07260 void addXandX2IfNonMissing(const TVec<T>& source, const TVec<int>& nnonmissing, TVec<T> somme, TVec<T> somme2) 07261 { 07262 #ifdef BOUNDCHECK 07263 if (source.length()!=nnonmissing.length() || source.length()!=somme.length() || source.length()!=somme2.length()) 07264 PLERROR("addIfNonMissing: all arguments should have the same length, got %d,%d,%d,%d\n", 07265 source.length(),nnonmissing.length(),somme.length(),somme2.length()); 07266 #endif 07267 T* s=source.data(); 07268 T* s1=somme.data(); 07269 T* s2=somme.data(); 07270 int* n=nnonmissing.data(); 07271 int size=source.length(); 07272 for (int i=0;i<size;i++) 07273 if (finite(s[i])) 07274 { 07275 s1[i] += s[i]; 07276 s2[i] += s[i]*s[i]; 07277 n[i]++; 07278 } 07279 } 07280 07281 // input_gradient[j] = sum_i weights[i,j]*output_gradient[i] 07282 // weights[i,j] -= learning_rate * output_gradient[i] * input[j] 07283 template<class T> 07284 void layerBpropUpdate(TVec<T> input_gradient, TMat<T> weights, const TVec<T>& input, 07285 const TVec<T>& output_gradient, real learning_rate) 07286 { 07287 int n_inputs = input_gradient.length(); 07288 int n_outputs = output_gradient.length(); 07289 #ifdef BOUNDCHECK 07290 if (weights.length() != n_outputs || weights.width() != n_inputs 07291 || input.length() != n_inputs) 07292 PLERROR("layerBpropUpdate: arguments have incompatible sizes"); 07293 #endif 07294 input_gradient.clear(); 07295 T* in_g = input_gradient.data(); 07296 T* out_g = output_gradient.data(); 07297 T* inp = input.data(); 07298 for (int i=0;i<n_outputs;i++) 07299 { 07300 T* Wi = weights[i]; 07301 T out_gi = out_g[i]; 07302 for (int j=0;j<n_inputs;j++) 07303 { 07304 in_g[j] += Wi[j] * out_gi; 07305 Wi[j] -= learning_rate * out_gi * inp[j]; 07306 } 07307 } 07308 } 07309 07310 07311 // input_gradient[j] = sum_i weights[i,j]*output_gradient[i] 07312 // weights[i,j] -= learning_rate * (output_gradient[i] * input[j] + weight_decay * weights[i,j]) 07313 template<class T> 07314 void layerL2BpropUpdate(TVec<T> input_gradient, TMat<T> weights, const TVec<T>& input, 07315 const TVec<T>& output_gradient, real learning_rate, T weight_decay) 07316 { 07317 int n_inputs = input_gradient.length(); 07318 int n_outputs = output_gradient.length(); 07319 #ifdef BOUNDCHECK 07320 if (weights.length() != n_outputs || weights.width() != n_inputs 07321 || input.length() != n_inputs) 07322 PLERROR("layerL2BpropUpdate: arguments have incompatible sizes"); 07323 #endif 07324 input_gradient.clear(); 07325 T* in_g = input_gradient.data(); 07326 T* out_g = output_gradient.data(); 07327 T* inp = input.data(); 07328 for (int i=0;i<n_outputs;i++) 07329 { 07330 T* Wi = weights[i]; 07331 T out_gi = out_g[i]; 07332 for (int j=0;j<n_inputs;j++) 07333 { 07334 T Wij = Wi[j]; 07335 in_g[j] += Wij * out_gi; 07336 Wi[j] -= learning_rate * (out_gi * inp[j] + weight_decay * Wij); 07337 } 07338 } 07339 } 07340 07341 // like layerL2BpropUpdate but weights is given transposed (not reflected in the formula below). 07342 // input_gradient[j] = sum_i weights[j,i]*output_gradient[i] 07343 // weights[i,j] -= learning_rate * (output_gradient[i] * input[j] + weight_decay * weights[i,j]) 07344 template<class T> 07345 void transposedLayerL2BpropUpdate(TVec<T> input_gradient, TMat<T> weights, const TVec<T>& input, 07346 const TVec<T>& output_gradient, real learning_rate, T weight_decay) 07347 { 07348 int n_inputs = input_gradient.length(); 07349 int n_outputs = output_gradient.length(); 07350 #ifdef BOUNDCHECK 07351 if (weights.width() != n_outputs || weights.length() != n_inputs 07352 || input.length() != n_inputs) 07353 PLERROR("layerL2BpropUpdate: arguments have incompatible sizes"); 07354 #endif 07355 input_gradient.clear(); 07356 T* in_g = input_gradient.data(); 07357 T* out_g = output_gradient.data(); 07358 T* inp = input.data(); 07359 for (int j=0;j<n_inputs;j++) 07360 { 07361 T* Wj = weights[j]; 07362 T inp_j = inp[j]; 07363 for (int i=0;i<n_outputs;i++) 07364 { 07365 T out_gi = out_g[i]; 07366 T Wji = Wj[i]; 07367 in_g[j] += Wji * out_gi; 07368 Wj[i] -= learning_rate * (out_gi * inp_j + weight_decay * Wji); 07369 } 07370 } 07371 } 07372 07373 // input_gradient[j] = sum_i weights[i,j]*output_gradient[i] 07374 // weights[i,j] -= learning_rate * (output_gradient[i] * input[j] + weight_decay * sign(weights[i,j])) 07375 template<class T> 07376 void layerL1BpropUpdate(TVec<T> input_gradient, TMat<T> weights, const TVec<T>& input, 07377 const TVec<T>& output_gradient, real learning_rate, T weight_decay) 07378 { 07379 int n_inputs = input_gradient.length(); 07380 int n_outputs = output_gradient.length(); 07381 #ifdef BOUNDCHECK 07382 if (weights.length() != n_outputs || weights.width() != n_inputs 07383 || input.length() != n_inputs) 07384 PLERROR("layerL1BpropUpdate: arguments have incompatible sizes"); 07385 #endif 07386 input_gradient.clear(); 07387 T* in_g = input_gradient.data(); 07388 T* out_g = output_gradient.data(); 07389 T* inp = input.data(); 07390 for (int i=0;i<n_outputs;i++) 07391 { 07392 T* Wi = weights[i]; 07393 T out_gi = out_g[i]; 07394 for (int j=0;j<n_inputs;j++) 07395 { 07396 T Wij = Wi[j]; 07397 in_g[j] += Wij * out_gi; 07398 Wi[j] -= learning_rate * (out_gi * inp[j] + weight_decay * sign(Wij)); 07399 } 07400 } 07401 } 07402 07403 // like layerL1BpropUpdate but weights is given transposed. 07404 // input_gradient[j] = sum_i weights[j,i]*output_gradient[i] 07405 // weights[i,j] -= learning_rate * (output_gradient[i] * input[j] + weight_decay * sign(weights[i,j])) 07406 template<class T> 07407 void transposedLayerL1BpropUpdate(TVec<T> input_gradient, TMat<T> weights, const TVec<T>& input, 07408 const TVec<T>& output_gradient, real learning_rate, T weight_decay) 07409 { 07410 int n_inputs = input_gradient.length(); 07411 int n_outputs = output_gradient.length(); 07412 #ifdef BOUNDCHECK 07413 if (weights.width() != n_outputs || weights.length() != n_inputs 07414 || input.length() != n_inputs) 07415 PLERROR("layerL1BpropUpdate: arguments have incompatible sizes"); 07416 #endif 07417 input_gradient.clear(); 07418 T* in_g = input_gradient.data(); 07419 T* out_g = output_gradient.data(); 07420 T* inp = input.data(); 07421 for (int j=0;j<n_inputs;j++) 07422 { 07423 T* Wj = weights[j]; 07424 T inp_j = inp[j]; 07425 for (int i=0;i<n_outputs;i++) 07426 { 07427 T out_gi = out_g[i]; 07428 T Wji = Wj[i]; 07429 in_g[j] += Wji * out_gi; 07430 Wj[i] -= learning_rate * (out_gi * inp_j + weight_decay * sign(Wji)); 07431 } 07432 } 07433 } 07434 07437 template<class T> 07438 void identityMatrix(TMat<T> m) 07439 { 07440 int l=m.length(); 07441 int w=m.width(); 07442 for (int i=0;i<l;i++) 07443 { 07444 T* mi = m[i]; 07445 for (int j=0;j<w;j++) 07446 if (j==i) 07447 mi[j]=1; 07448 else 07449 mi[j]=0; 07450 } 07451 } 07452 07454 template<class T> 07455 TMat<T> identityMatrix(int n, int m=-1) 07456 { 07457 if (m<0) m=n; 07458 TMat<T> result(n,m); 07459 identityMatrix(result); 07460 return result; 07461 } 07462 07463 07464 } // end of namespace PLearn 07465 07466 07467 // Norman: replaced the code below with this wrapper 07468 SET_HASH_FUNCTION(PLearn::TVec<T>, T, v, sumsquare(v)) 07469 SET_HASH_WITH_FUNCTION(PLearn::Vec, v, sumsquare(v)) 07470 07471 //#if __GNUC__==3 && __GNUC_MINOR__>0 07472 //namespace __gnu_cxx { 07473 //#else 07474 //namespace std { 07475 //#endif 07476 // 07477 //template<class T> 07478 //struct hash<PLearn::TVec<T> > 07479 //{ 07480 // size_t operator()(PLearn::TVec<T> v) const { return hash<T>()(sumsquare(v));} 07481 //}; 07482 07483 //} // end of namespace std 07484 07485 07486 #endif // TMat_maths_impl_H 07487 07488 07489 /* 07490 Local Variables: 07491 mode:c++ 07492 c-basic-offset:4 07493 c-file-style:"stroustrup" 07494 c-file-offsets:((innamespace . 0)(inline-open . 0)) 07495 indent-tabs-mode:nil 07496 fill-column:79 07497 End: 07498 */ 07499 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :