Scopira 20080306

vectormath.h

00001 
00002 /*
00003  *  Copyright (c) 2002-2011    National Research Council
00004  *
00005  *  All rights reserved.
00006  *
00007  *  This material is confidential and proprietary information of
00008  *  National Research Council Canada ("Confidential Information").
00009  *  This Confidential Information may only be used and reproduced
00010  *  in accordance with the terms of the license agreement.
00011  *
00012  */
00013 
00014 #ifndef __INCLUDED_SCOPIRA_BASEKIT_VECTORMATH_H__
00015 #define __INCLUDED_SCOPIRA_BASEKIT_VECTORMATH_H__
00016 
00017 #include <math.h>
00018 
00019 #include <scopira/basekit/narray.h>
00020 
00021 // THIS FILE HAS BEEN FULLY DOCUMENTED
00022 
00023 namespace scopira
00024 {
00025   namespace basekit
00026   {
00035     template <class T>
00036       void add_scalar(T &v, typename T::data_type val);
00045     template <class T>
00046       void sub_scalar(T &v, typename T::data_type val);
00055     template <class T>
00056       void mul_scalar(T &v, typename T::data_type val);
00065     template <class T>
00066       void div_scalar(T &v, typename T::data_type val);
00076     template <class T>
00077       void invdiv_scalar(T &v, typename T::data_type val);
00078  
00079 
00080     // not counst parms for muffinization
00081     // probably wont need a switcher function because of micro-functions
00082     // median (TODO)
00083 
00092     template <class V, class T>
00093       size_t min(const V &v, T& out);
00102     template <class V, class T>
00103       size_t max(const V &v, T& out);
00111     template <class V, class T>
00112       void sum(const V &v, T& out);
00120     template <class V, class T>
00121       void mean(const V &v, T& out);
00129     template <class V, class T>
00130       void variance(const V &v, T& out);
00138     template <class V, class T>
00139       void stddev(const V &v, T& out)
00140         { variance<V,T>(v, out); out = static_cast<T>(::sqrt(out)); }
00148     template <class LH, class RH>
00149       void add_vector(LH &target, const RH& delta);
00157     template <class LH, class RH>
00158       void sub_vector(LH &target, const RH& delta);
00166     template <class LH, class RH>
00167       void mul_vector(LH &target, const RH& delta);
00175     template <class LH, class RH>
00176       void div_vector(LH &target, const RH& delta);
00187     template <class LH, class RH>
00188       void abs_vector(LH &target, const RH& source);
00196     template <class V>
00197       void print_vector(scopira::tool::oflow_i &o, const V &vv);
00205     template <class T>
00206       void pearson_scalar(T &x, T &y, double& rho);
00207 
00216     template <class T>
00217       void construct_index(narray<T> &array, narray<int> &index);  
00218 
00219     template <class DATAT, class REFT>
00220       class vecindex_sortable;
00221   }
00222 }
00223 
00224 template <class T>
00225   void scopira::basekit::add_scalar(T &v, typename T::data_type val)
00226 {
00227   typedef typename T::iterator iter;
00228   iter ii=v.begin(), endii = v.end();
00229 
00230   for (; ii != endii; ++ii)
00231     (*ii) += val;
00232 }
00233 
00234 template <class T>
00235   void scopira::basekit::sub_scalar(T &v, typename T::data_type val)
00236 {
00237   typedef typename T::iterator iter;
00238   iter ii=v.begin(), endii = v.end();
00239 
00240   for (; ii != endii; ++ii)
00241     (*ii) -= val;
00242 }
00243 
00244 template <class T>
00245   void scopira::basekit::mul_scalar(T &v, typename T::data_type val)
00246 {
00247   typedef typename T::iterator iter;
00248   iter ii=v.begin(), endii = v.end();
00249 
00250   for (; ii != endii; ++ii)
00251     (*ii) *= val;
00252 }
00253 
00254 template <class T>
00255   void scopira::basekit::div_scalar(T &v, typename T::data_type val)
00256 {
00257   typedef typename T::iterator iter;
00258   iter ii=v.begin(), endii = v.end();
00259 
00260   for (; ii != endii; ++ii)
00261     (*ii) /= val;
00262 }
00263 
00264 template <class T>
00265   void scopira::basekit::invdiv_scalar(T &v, typename T::data_type val)
00266 {
00267   typedef typename T::iterator iter;
00268   iter ii=v.begin(), endii = v.end();
00269 
00270   for (; ii != endii; ++ii)
00271     (*ii) = val / (*ii);
00272 }
00273 
00274 template <class V, class T>
00275   size_t scopira::basekit::min(const V &v, T& out)
00276 {
00277   typedef typename V::const_iterator iter;
00278   iter ii, endii;
00279   size_t idx, cur;
00280   T m;
00281 
00282   ii = v.begin();
00283   endii = v.end();
00284   assert(ii != endii);
00285   m = *ii;
00286   cur = idx = 0;
00287   for (++ii, ++cur; ii!=endii; ++ii, ++cur)
00288     if (*ii < m) {
00289       m = *ii;
00290       idx = cur;
00291     }
00292   out = m;
00293   return idx;
00294 }
00295 
00296 template <class V, class T>
00297   size_t scopira::basekit::max(const V &v, T& out)
00298 {
00299   typedef typename V::const_iterator iter;
00300   iter ii, endii;
00301   size_t idx, cur;
00302   T m;
00303 
00304   ii = v.begin();
00305   endii = v.end();
00306   assert(ii != endii);
00307   m = *ii;
00308   cur = idx = 0;
00309   for (++ii, ++cur; ii!=endii; ++ii, ++cur)
00310     if (*ii > m) {
00311       m = *ii;
00312       idx = cur;
00313     }
00314   out = m;
00315   return idx;
00316 }
00317 
00318 template <class V, class T>
00319   void scopira::basekit::sum(const V &v, T& out)
00320 {
00321   typedef typename V::const_iterator iter;
00322   iter ii, endii;
00323   T m;
00324 
00325   m = 0;
00326   ii = v.begin();
00327   endii = v.end();
00328   for (; ii!=endii; ++ii)
00329     m += *ii;
00330   out = m;
00331 }
00332 
00333 template <class V, class T>
00334   void scopira::basekit::mean(const V &v, T& out)
00335 {
00336   if (v.empty()) {
00337     out = 0;
00338     return;
00339   }
00340   sum(v, out);
00341   out = out / v.size();
00342 }
00343 
00344 template <class V, class T>
00345   void scopira::basekit::variance(const V &v, T& out)
00346 {
00347   typedef typename V::const_iterator iter;
00348   iter ii, endii;
00349   size_t mx;
00350   T avg, sum, tt;
00351 
00352   mx = v.size();
00353   if (mx <= 1) {
00354     out = 0; //special, stupid case
00355     return;
00356   }
00357 
00358   mean(v, avg);
00359 
00360   sum = 0;
00361   endii = v.end();
00362   for (ii=v.begin(); ii != endii; ++ii) {
00363     tt = *ii - avg;
00364     sum += tt * tt;
00365   }
00366 
00367   out = sum / (mx - 1);
00368 }
00369 
00370 template <class LH, class RH>
00371   void scopira::basekit::add_vector(LH &target, const RH& delta)
00372 {
00373   typename LH::iterator tar, endtar;
00374   typename RH::const_iterator src, endsrc;
00375 
00376   tar = target.begin();
00377   endtar = target.end();
00378   src = delta.begin();
00379   endsrc = delta.end();
00380   while (tar != endtar) {
00381     assert(src != endsrc);
00382     *tar += *src;
00383     ++tar;
00384     ++src;
00385   }
00386 }
00387 
00388 template <class LH, class RH>
00389   void scopira::basekit::sub_vector(LH &target, const RH& delta)
00390 {
00391   typename LH::iterator tar, endtar;
00392   typename RH::const_iterator src, endsrc;
00393 
00394   tar = target.begin();
00395   endtar = target.end();
00396   src = delta.begin();
00397   endsrc = delta.end();
00398   while (tar != endtar) {
00399     assert(src != endsrc);
00400     *tar -= *src;
00401     ++tar;
00402     ++src;
00403   }
00404 }
00405 
00406 template <class LH, class RH>
00407   void scopira::basekit::mul_vector(LH &target, const RH& delta)
00408 {
00409   typename LH::iterator tar, endtar;
00410   typename RH::const_iterator src, endsrc;
00411 
00412   tar = target.begin();
00413   endtar = target.end();
00414   src = delta.begin();
00415   endsrc = delta.end();
00416   while (tar != endtar) {
00417     assert(src != endsrc);
00418     *tar *= *src;
00419     ++tar;
00420     ++src;
00421   }
00422 }
00423 
00424 template <class LH, class RH>
00425   void scopira::basekit::div_vector(LH &target, const RH& delta)
00426 {
00427   typename LH::iterator tar, endtar;
00428   typename RH::const_iterator src, endsrc;
00429 
00430   tar = target.begin();
00431   endtar = target.end();
00432   src = delta.begin();
00433   endsrc = delta.end();
00434   while (tar != endtar) {
00435     assert(src != endsrc);
00436     *tar /= *src;
00437     ++tar;
00438     ++src;
00439   }
00440 }
00441 
00442 template <class LH, class RH>
00443   void scopira::basekit::abs_vector(LH &target, const RH& source)
00444 {
00445   typename LH::iterator tar, endtar;
00446   typename RH::const_iterator src, endsrc;
00447 
00448   tar = target.begin();
00449   endtar = target.end();
00450   src = source.begin();
00451   endsrc = source.end();
00452   while (tar != endtar) {
00453     assert(src != endsrc);
00454     *tar = (*src < 0.0) ? (*src * -1) : *src;
00455     ++tar;
00456     ++src;
00457   }
00458 }
00459 
00460 template <class V>
00461   void scopira::basekit::print_vector(scopira::tool::oflow_i &o, const V &vv)
00462 {
00463   typedef typename V::const_iterator iter;
00464   iter ii = vv.begin(), endii = vv.end();
00465   o << '(';
00466   for (; ii != endii; ++ii)
00467     o << *ii << ',';
00468   o << ")\n";
00469 }
00470 
00482 // see Evident Scorpio implementation
00483 template <class T>
00484   void scopira::basekit::pearson_scalar(T &x, T &y, double& rho)
00485 {
00486   double sX;    // sum of X pts
00487   double sXX;   // sum of squares
00488   double sXY;
00489   double sY;
00490   double sYY;
00491   double top, btm;
00492 
00493   sX = sXX = sY = sYY = sXY = 0;
00494   int n = x.size();
00495 
00496   assert(x.size() == y.size());
00497 
00498   for (long i=0; i<n; i++) 
00499   {
00500     sX += x[i];
00501     sXX += x[i] * x[i];
00502     sXY += x[i] * y[i];
00503     sY += y[i];
00504     sYY += y[i] * y[i];
00505   }
00506 
00507   top = sXY - (sX*sY/n);
00508   btm = (sXX-sX*sX/n) * (sYY-sY*sY/n);
00509 
00510   if (top==0 || btm==0)
00511     rho = 0;
00512   else
00513     rho = top / sqrt(btm);
00514 }
00515 
00516 // constructs an index for a vector
00517 template <class T>
00518   void scopira::basekit::construct_index(narray<T> &array, narray<int> &index)
00519 {
00520   int i, j, k, n, l=0, ir, temp_index, temp;
00521   int stack_index = -1;
00522   T array_val;
00523   narray<int> stack;
00524   
00525   stack.resize(50);
00526   n = array.size();
00527   ir = n-1;
00528   
00529   for (j=0; j<n; j++)
00530     index[j] = j;
00531   while (true) {
00532     if (ir-l < 7) {
00533       for (j=l+1; j<=ir; j++) {
00534         temp_index = index[j];
00535         array_val = array[temp_index];
00536         for (i=j-1; i>=1; i--) {
00537           if (array[index[i]] <= array_val)
00538             break;
00539           index[i+1] = index[i]; 
00540         }
00541         index[i+1] = temp_index;
00542       }
00543       if (stack_index < 0)
00544         break;
00545       ir = stack[stack_index--];
00546       l = stack[stack_index--];
00547     } else {
00548       k = (l+ir) >> 1;
00549       temp = index[k];  index[k] = index[l+1];  index[l+1] = temp;
00550       if (array[index[l+1]] > array[index[ir]]) {
00551         temp = index[l+1];  index[l+1] = index[ir];  index[ir] = temp;
00552       }
00553       if (array[index[l]] > array[index[ir]]) {
00554         temp = index[l];  index[l] = index[ir];  index[ir] = temp;
00555       }
00556       if (array[index[l+1]] > array[index[l]]) {
00557         temp = index[l+1];  index[l+1] = index[l];  index[l] = temp;
00558       }
00559       i = l+1;
00560       j = ir;
00561       temp_index = index[l+1];
00562       array_val = array[temp_index];
00563       while (true) {
00564         i++;
00565         while (array[index[i]] < array_val)
00566           i++;
00567         j--;
00568         while (array[index[j]] > array_val)
00569           j--;
00570         if (j < i)
00571           break;
00572         temp = index[i];  index[i] = index[j];  index[j] = temp;
00573       }
00574       index[l+1] = index[j];
00575       index[j] = temp_index;
00576       stack_index += 2;
00577       if (stack_index >= 50) {
00578         //KERROR << "Stack size too small in constructing index\n";
00579         return;
00580       }
00581       if (ir-i+1 >- j-l) {
00582         stack[stack_index] = ir;
00583         stack[stack_index-1] = i;
00584         ir = j-1;
00585       } else {
00586         stack[stack_index] = j-1;
00587         stack[stack_index-1] = l;
00588         l = i;
00589       }
00590     }
00591   }
00592 }
00593 
00604 template <class DATAT, class REFT>
00605   class scopira::basekit::vecindex_sortable
00606 {
00607   public:
00608     typedef DATAT data_type;
00609     typedef REFT ref_type;
00610   private:
00611     data_type &dm_data;
00612     const ref_type &dm_ref;
00613   public:
00614     vecindex_sortable(DATAT &data, REFT &ref) : dm_data(data), dm_ref(ref) { }
00615     int compare_element(int lidx, int ridx) const {
00616       typedef typename ref_type::data_type t;
00617       t l, r;
00618 
00619       l = dm_ref[dm_data[lidx]];
00620       r = dm_ref[dm_data[ridx]];
00621       if (l<r)
00622         return -1;
00623       else if (l>r)
00624         return 1;
00625       else
00626         return 0;
00627     }
00628     void swap_element(int lidx, int ridx) {
00629       typedef typename data_type::data_type t;
00630       t tmp;
00631 
00632       tmp = dm_data[lidx];
00633       dm_data[lidx] = dm_data[ridx];
00634       dm_data[ridx] = tmp;
00635     }
00636 };
00637 
00638 #endif
00639 
00640