PLearn 0.1
TMat_maths_impl.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // PLearn (A C++ Machine Learning Library)
00004 // Copyright (C) 1998 Pascal Vincent
00005 // Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal
00006 // 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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines