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()) || (arma_isfinite(negdelta) == false) )
00200     {
00201     return log_a;
00202     }
00203   else
00204     {
00205     #if defined(ARMA_HAVE_LOG1P)
00206       return (log_a + log1p(std::exp(negdelta)));
00207     #else
00208       return (log_a + std::log(1.0 + std::exp(negdelta)));
00209     #endif
00210     }
00211   }
00212 
00213 
00214 
00215 template<typename eT>
00216 inline
00217 eT
00218 trunc_log(const eT x)
00219   {
00220   if(std::numeric_limits<eT>::is_iec559)
00221     {
00222     if(x == std::numeric_limits<eT>::infinity())
00223       {
00224       return Math<eT>::log_max();
00225       }
00226     if(x <= 0)
00227       {
00228       return Math<eT>::log_min();
00229       }
00230     }
00231   else
00232     {
00233     return std::log(x);
00234     }
00235   }
00236 
00237 
00238 
00239 template<typename eT>
00240 inline
00241 eT
00242 trunc_exp(const eT x)
00243   {
00244   if(std::numeric_limits<eT>::is_iec559 && (x >= Math<eT>::log_max() ))
00245     {
00246     return std::numeric_limits<eT>::max();
00247     }
00248   else
00249     {
00250     return std::exp(x);
00251     }
00252   }
00253 
00254 
00255 
00256 //
00257 // log
00258 
00259 template<typename T1>
00260 inline
00261 const Op<T1, op_log>
00262 log(const Base<typename T1::elem_type,T1>& A)
00263   {
00264   arma_extra_debug_sigprint();
00265   
00266   return Op<T1, op_log>(A.get_ref());
00267   }
00268 
00269 
00270 
00271 //
00272 // trunc_log
00273 
00274 template<typename T1>
00275 inline
00276 const Op<T1, op_trunc_log>
00277 trunc_log(const Base<typename T1::elem_type,T1>& A)
00278   {
00279   arma_extra_debug_sigprint();
00280   
00281   return Op<T1, op_trunc_log>(A.get_ref());
00282   }
00283 
00284 
00285 
00286 //
00287 // log10
00288 
00289 template<typename T1>
00290 inline
00291 const Op<T1, op_log10>
00292 log10(const Base<typename T1::elem_type,T1>& A)
00293   {
00294   arma_extra_debug_sigprint();
00295   
00296   return Op<T1, op_log10>(A.get_ref());
00297   }
00298 
00299 
00300 
00301 //
00302 // exp
00303 
00304 template<typename T1>
00305 inline
00306 const Op<T1, op_exp>
00307 exp(const Base<typename T1::elem_type,T1>& A)
00308   {
00309   arma_extra_debug_sigprint();
00310   
00311   return Op<T1, op_exp>(A.get_ref());
00312   }
00313 
00314 
00315 
00316 //
00317 // trunc_exp
00318 
00319 template<typename T1>
00320 inline
00321 const Op<T1, op_trunc_exp>
00322 trunc_exp(const Base<typename T1::elem_type,T1>& A)
00323   {
00324   arma_extra_debug_sigprint();
00325   
00326   return Op<T1, op_trunc_exp>(A.get_ref());
00327   }
00328 
00329 
00330 
00331 //
00332 // abs
00333 
00334 template<typename T1>
00335 inline
00336 Mat<typename T1::pod_type>
00337 abs(const Base<typename T1::elem_type,T1>& X)
00338   {
00339   arma_extra_debug_sigprint();
00340   
00341   const unwrap<T1> A_tmp(X.get_ref());
00342 
00343   // if T1 is a complex matrix,
00344   // pod_type is the underlying type used by std::complex;
00345   // otherwise pod_type is the same as elem_type
00346   
00347   typedef typename T1::elem_type  in_eT;
00348   typedef typename T1::pod_type  out_eT;
00349 
00350   const Mat<in_eT>& A = A_tmp.M;
00351   
00352   Mat<out_eT> out(A.n_rows, A.n_cols);
00353   
00354   const in_eT* A_mem   = A.mem;
00355   out_eT*      out_mem = out.memptr();
00356   
00357   for(u32 i=0; i<out.n_elem; ++i)
00358     {
00359     out_mem[i] = std::abs(A_mem[i]);
00360     }
00361   
00362   return out;
00363   }
00364 
00365 
00366 
00367 //
00368 // fabs
00369 
00370 template<typename T1>
00371 inline
00372 Mat<typename T1::pod_type>
00373 fabs(const Base<typename T1::elem_type,T1>& A)
00374   {
00375   arma_extra_debug_sigprint();
00376   
00377   return abs(A);
00378   }
00379 
00380 
00381 
00382 //
00383 // square
00384 
00385 template<typename T1>
00386 inline
00387 const Op<T1, op_square>
00388 square(const Base<typename T1::elem_type,T1>& A)
00389   {
00390   arma_extra_debug_sigprint();
00391   
00392   return Op<T1, op_square>(A.get_ref());
00393   }
00394 
00395 
00396 
00397 //
00398 // sqrt
00399 
00400 template<typename T1>
00401 inline
00402 const Op<T1, op_sqrt>
00403 sqrt(const Base<typename T1::elem_type,T1>& A)
00404   {
00405   arma_extra_debug_sigprint();
00406   
00407   return Op<T1, op_sqrt>(A.get_ref());
00408   }
00409 
00410 
00411 
00412 // pow
00413 
00414 template<typename T1>
00415 inline
00416 const Op<T1, op_pow>
00417 pow(const Base<typename T1::elem_type,T1>& A, const typename T1::elem_type exponent)
00418   {
00419   arma_extra_debug_sigprint();
00420   
00421   return Op<T1, op_pow>(A.get_ref(), exponent);
00422   }
00423 
00424 
00425 
00426 // pow, specialised handling (non-complex exponent for complex matrices)
00427 
00428 template<typename T1>
00429 inline
00430 const Op<T1, op_pow>
00431 pow(const Base<typename T1::elem_type,T1>& A, const typename T1::elem_type::value_type exponent)
00432   {
00433   arma_extra_debug_sigprint();
00434   
00435   return Op<T1, op_pow>(A.get_ref(), eT(exponent));
00436   }
00437 
00438 
00439 
00440 // pow_s32  (integer exponent)
00441 
00442 template<typename T1>
00443 inline
00444 const Op<T1, op_pow_s32>
00445 pow(const Base<typename T1::elem_type,T1>& A, const s32 exponent)
00446   {
00447   arma_extra_debug_sigprint();
00448   
00449   if(exponent >= 0)
00450     {
00451     return Op<T1, op_pow_s32>(A.get_ref(), exponent, 0);
00452     }
00453   else
00454     {
00455     return Op<T1, op_pow_s32>(A.get_ref(), -exponent, 1);
00456     }
00457   }
00458 
00459 
00460 
00461 // conj
00462 
00463 template<typename T, typename T1>
00464 inline
00465 const Op<T1, op_conj>
00466 conj(const Base<std::complex<T>,T1>& A)
00467   {
00468   arma_extra_debug_sigprint();
00469 
00470   return Op<T1, op_conj>(A.get_ref());
00471   }
00472 
00473 
00474 
00475 template<typename T1>
00476 inline
00477 const T1&
00478 conj(const Op<T1, op_conj>& A)
00479   {
00480   arma_extra_debug_sigprint();
00481   
00482   return A.m;
00483   }
00484 
00485 
00486 
00487 //! the conjugate of the transpose of a complex matrix is the same as the hermitian transpose
00488 template<typename T1>
00489 inline
00490 const Op<T1, op_htrans>
00491 conj(const Op<T1, op_trans>& A)
00492   {
00493   arma_extra_debug_sigprint();
00494   
00495   arma_type_check< is_complex<typename T1::elem_type>::value == false >::apply();
00496 
00497   return Op<T1, op_htrans>(A.m);
00498   }
00499 
00500 
00501 //! @}