fn_misc.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2009 NICTA
00002 // 
00003 // Authors:
00004 // - Conrad Sanderson (conradsand at ieee dot org)
00005 // 
00006 // This file is part of the Armadillo C++ library.
00007 // It is provided without any warranty of fitness
00008 // for any purpose. You can redistribute this file
00009 // and/or modify it under the terms of the GNU
00010 // Lesser General Public License (LGPL) as published
00011 // by the Free Software Foundation, either version 3
00012 // of the License or (at your option) any later version.
00013 // (see http://www.opensource.org/licenses for more info)
00014 
00015 
00016 //! \addtogroup fn_misc
00017 //! @{
00018 
00019 //! \brief
00020 //! Generate a vector with 'num' elements.
00021 //! The values of the elements linearly increase from 'start' upto (and including) 'end'.
00022 
00023 template<typename eT>
00024 inline
00025 Mat<eT>
00026 linspace(const eT start, const eT end, const u32 num, const u32 dim = 0)
00027   {
00028   arma_extra_debug_sigprint();
00029   
00030   arma_debug_check( (num < 2), "linspace(): num must be >= 2");
00031   
00032   Mat<eT> x;
00033   
00034   if(dim == 0)
00035     {
00036     x.set_size(num,1);  // column vector
00037     }
00038   else
00039     {
00040     x.set_size(1,num);  // row vector
00041     }
00042   
00043   
00044   const eT delta = (end-start)/(num-1);
00045   
00046   x[0] = start;
00047   
00048   for(u32 i=1; i<num; ++i)
00049     {
00050     x[i] = x[i-1] + delta;
00051     }
00052   
00053   return x; 
00054   }
00055 
00056 
00057 
00058 inline
00059 mat
00060 linspace(const double start, const double end, const u32 num, const u32 dim = 0)
00061   {
00062   arma_extra_debug_sigprint();
00063   return linspace<double>(start, end, num, dim);
00064   }
00065 
00066 
00067 
00068 
00069 
00070 
00071 //
00072 // reshape
00073 
00074 template<typename T1>
00075 inline
00076 Mat<typename T1::elem_type>
00077 reshape(const Base<typename T1::elem_type,T1>& X, const u32 in_n_rows, const u32 in_n_cols, const u32 dim = 0)
00078   {
00079   arma_extra_debug_sigprint();
00080 
00081   typedef typename T1::elem_type eT;
00082   
00083   Mat<eT> out;
00084 
00085   const unwrap<T1> A_tmp(X.get_ref());
00086   const Mat<eT>& A = A_tmp.M;
00087 
00088   const u32 in_n_elem = in_n_rows * in_n_cols;
00089 
00090   arma_debug_check( (A.n_elem != in_n_elem), "reshape(): incompatible dimensions");
00091   arma_debug_check( (dim > 1), "reshape(): dim must be 0 or 1");
00092 
00093   if(dim == 0)
00094     {
00095     out = A;
00096 
00097     access::rw(out.n_rows) = in_n_rows;
00098     access::rw(out.n_cols) = in_n_cols;
00099 
00100     return out;
00101     }
00102   else
00103     {
00104     out.set_size(in_n_rows, in_n_cols);
00105     
00106     eT* out_mem = out.memptr();
00107     u32 i = 0;
00108     
00109     for(u32 row=0; row<A.n_rows; ++row)
00110       {
00111       for(u32 col=0; col<A.n_cols; ++col)
00112         {
00113         out_mem[i] = A.at(row,col);
00114         ++i;
00115         }
00116       }
00117     
00118     return out;
00119     }
00120   }
00121 
00122 
00123 
00124 //
00125 // real
00126 
00127 template<typename T, typename T1>
00128 inline
00129 Mat<T>
00130 real(const Base<std::complex<T>, T1>& X)
00131   {
00132   arma_extra_debug_sigprint();
00133   
00134   typedef std::complex<T> eT;
00135 
00136   const unwrap<T1> A_tmp(X.get_ref());
00137   const Mat<eT>& A = A_tmp.M;
00138   
00139   Mat<T> out(A.n_rows, A.n_cols);
00140   
00141   const eT* A_mem = A.mem;
00142   T* out_mem = out.memptr();
00143   
00144   for(u32 i=0; i<out.n_elem; ++i)
00145     {
00146     out_mem[i] = std::real(A_mem[i]);
00147     }
00148   
00149   return out;
00150   }
00151 
00152 
00153 
00154 //
00155 // imag
00156 
00157 template<typename T, typename T1>
00158 inline
00159 Mat<T>
00160 imag(const Base<std::complex<T>,T1>& X)
00161   {
00162   arma_extra_debug_sigprint();
00163   
00164   typedef std::complex<T> eT;
00165 
00166   const unwrap<T1> A_tmp(X.get_ref());
00167   const Mat<eT>& A = A_tmp.M;
00168   
00169   Mat<T> out(A.n_rows, A.n_cols);
00170   
00171   const eT* A_mem = A.mem;
00172   T* out_mem = out.memptr();
00173   
00174   for(u32 i=0; i<out.n_elem; ++i)
00175     {
00176     out_mem[i] = std::imag(A_mem[i]);
00177     }
00178   
00179   return out;
00180   }
00181 
00182 
00183 
00184 //
00185 // log_add
00186 
00187 template<typename eT>
00188 inline
00189 eT
00190 log_add(eT log_a, eT log_b)
00191   {
00192   if(log_a < log_b)
00193     {
00194     std::swap(log_a, log_b);
00195     }
00196   
00197   const eT negdelta = log_b - log_a;
00198   
00199   if( (negdelta < Math<eT>::log_min()) || std::isnan(negdelta) )
00200     {
00201     return log_a;
00202     }
00203   else
00204     {
00205     return (log_a + log1p(std::exp(negdelta)));
00206     }
00207   }
00208 
00209 
00210 
00211 template<typename eT>
00212 inline
00213 eT
00214 trunc_log(const eT x)
00215   {
00216   if(std::numeric_limits<eT>::is_iec559)
00217     {
00218     if(x == std::numeric_limits<eT>::infinity())
00219       {
00220       return Math<eT>::log_max();
00221       }
00222     if(x <= 0)
00223       {
00224       return Math<eT>::log_min();
00225       }
00226     }
00227   else
00228     {
00229     return std::log(x);
00230     }
00231   }
00232 
00233 
00234 
00235 template<typename eT>
00236 inline
00237 eT
00238 trunc_exp(const eT x)
00239   {
00240   if(std::numeric_limits<eT>::is_iec559 && (x >= Math<eT>::log_max() ))
00241     {
00242     return std::numeric_limits<eT>::max();
00243     }
00244   else
00245     {
00246     return std::exp(x);
00247     }
00248   }
00249 
00250 
00251 
00252 //
00253 // log
00254 
00255 template<typename T1>
00256 inline
00257 const Op<T1, op_log>
00258 log(const Base<typename T1::elem_type,T1>& A)
00259   {
00260   arma_extra_debug_sigprint();
00261   
00262   return Op<T1, op_log>(A.get_ref());
00263   }
00264 
00265 
00266 
00267 //
00268 // trunc_log
00269 
00270 template<typename T1>
00271 inline
00272 const Op<T1, op_trunc_log>
00273 trunc_log(const Base<typename T1::elem_type,T1>& A)
00274   {
00275   arma_extra_debug_sigprint();
00276   
00277   return Op<T1, op_trunc_log>(A.get_ref());
00278   }
00279 
00280 
00281 
00282 //
00283 // log10
00284 
00285 template<typename T1>
00286 inline
00287 const Op<T1, op_log10>
00288 log10(const Base<typename T1::elem_type,T1>& A)
00289   {
00290   arma_extra_debug_sigprint();
00291   
00292   return Op<T1, op_log10>(A.get_ref());
00293   }
00294 
00295 
00296 
00297 //
00298 // exp
00299 
00300 template<typename T1>
00301 inline
00302 const Op<T1, op_exp>
00303 exp(const Base<typename T1::elem_type,T1>& A)
00304   {
00305   arma_extra_debug_sigprint();
00306   
00307   return Op<T1, op_exp>(A.get_ref());
00308   }
00309 
00310 
00311 
00312 //
00313 // trunc_exp
00314 
00315 template<typename T1>
00316 inline
00317 const Op<T1, op_trunc_exp>
00318 trunc_exp(const Base<typename T1::elem_type,T1>& A)
00319   {
00320   arma_extra_debug_sigprint();
00321   
00322   return Op<T1, op_trunc_exp>(A.get_ref());
00323   }
00324 
00325 
00326 
00327 //
00328 // abs
00329 
00330 template<typename T1>
00331 inline
00332 Mat<typename T1::pod_type>
00333 abs(const Base<typename T1::elem_type,T1>& X)
00334   {
00335   arma_extra_debug_sigprint();
00336   
00337   const unwrap<T1> A_tmp(X.get_ref());
00338 
00339   // if T1 is a complex matrix,
00340   // pod_type is the underlying type used by std::complex;
00341   // otherwise pod_type is the same as elem_type
00342   
00343   typedef typename T1::elem_type  in_eT;
00344   typedef typename T1::pod_type  out_eT;
00345 
00346   const Mat<in_eT>& A = A_tmp.M;
00347   
00348   Mat<out_eT> out(A.n_rows, A.n_cols);
00349   
00350   const in_eT* A_mem   = A.mem;
00351   out_eT*      out_mem = out.memptr();
00352   
00353   for(u32 i=0; i<out.n_elem; ++i)
00354     {
00355     out_mem[i] = std::abs(A_mem[i]);
00356     }
00357   
00358   return out;
00359   }
00360 
00361 
00362 
00363 //
00364 // fabs
00365 
00366 template<typename T1>
00367 inline
00368 Mat<typename T1::pod_type>
00369 fabs(const Base<typename T1::elem_type,T1>& A)
00370   {
00371   arma_extra_debug_sigprint();
00372   
00373   return abs(A);
00374   }
00375 
00376 
00377 
00378 //
00379 // square
00380 
00381 template<typename T1>
00382 inline
00383 const Op<T1, op_square>
00384 square(const Base<typename T1::elem_type,T1>& A)
00385   {
00386   arma_extra_debug_sigprint();
00387   
00388   return Op<T1, op_square>(A.get_ref());
00389   }
00390 
00391 
00392 
00393 //
00394 // sqrt
00395 
00396 template<typename T1>
00397 inline
00398 const Op<T1, op_sqrt>
00399 sqrt(const Base<typename T1::elem_type,T1>& A)
00400   {
00401   arma_extra_debug_sigprint();
00402   
00403   return Op<T1, op_sqrt>(A.get_ref());
00404   }
00405 
00406 
00407 
00408 // pow
00409 
00410 template<typename T1>
00411 inline
00412 const Op<T1, op_pow>
00413 pow(const Base<typename T1::elem_type,T1>& A, const typename T1::elem_type exponent)
00414   {
00415   arma_extra_debug_sigprint();
00416   
00417   return Op<T1, op_pow>(A.get_ref(), exponent);
00418   }
00419 
00420 
00421 
00422 // pow, specialised handling (non-complex exponent for complex matrices)
00423 
00424 template<typename T1>
00425 inline
00426 const Op<T1, op_pow>
00427 pow(const Base<typename T1::elem_type,T1>& A, const typename T1::elem_type::value_type exponent)
00428   {
00429   arma_extra_debug_sigprint();
00430   
00431   return Op<T1, op_pow>(A.get_ref(), eT(exponent));
00432   }
00433 
00434 
00435 
00436 // pow_s32  (integer exponent)
00437 
00438 template<typename T1>
00439 inline
00440 const Op<T1, op_pow_s32>
00441 pow(const Base<typename T1::elem_type,T1>& A, const s32 exponent)
00442   {
00443   arma_extra_debug_sigprint();
00444   
00445   if(exponent >= 0)
00446     {
00447     return Op<T1, op_pow_s32>(A.get_ref(), exponent, 0);
00448     }
00449   else
00450     {
00451     return Op<T1, op_pow_s32>(A.get_ref(), -exponent, 1);
00452     }
00453   }
00454 
00455 
00456 
00457 // conj
00458 
00459 template<typename T, typename T1>
00460 inline
00461 const Op<T1, op_conj>
00462 conj(const Base<std::complex<T>,T1>& A)
00463   {
00464   arma_extra_debug_sigprint();
00465 
00466   return Op<T1, op_conj>(A.get_ref());
00467   }
00468 
00469 
00470 
00471 template<typename T1>
00472 inline
00473 const T1&
00474 conj(const Op<T1, op_conj>& A)
00475   {
00476   arma_extra_debug_sigprint();
00477   
00478   return A.m;
00479   }
00480 
00481 
00482 
00483 //! the conjugate of the transpose of a complex matrix is the same as the hermitian transpose
00484 template<typename T1>
00485 inline
00486 const Op<T1, op_htrans>
00487 conj(const Op<T1, op_trans>& A)
00488   {
00489   arma_extra_debug_sigprint();
00490   
00491   arma_type_check< is_complex<typename T1::elem_type>::value == false >::apply();
00492 
00493   return Op<T1, op_htrans>(A.m);
00494   }
00495 
00496 
00497 //! @}