$search
00001 // Copyright (C) 2009-2011 NICTA (www.nicta.com.au) 00002 // Copyright (C) 2009-2011 Conrad Sanderson 00003 // Copyright (C) 2009-2010 Dimitrios Bouzas 00004 // 00005 // This file is part of the Armadillo C++ library. 00006 // It is provided without any warranty of fitness 00007 // for any purpose. You can redistribute this file 00008 // and/or modify it under the terms of the GNU 00009 // Lesser General Public License (LGPL) as published 00010 // by the Free Software Foundation, either version 3 00011 // of the License or (at your option) any later version. 00012 // (see http://www.opensource.org/licenses for more info) 00013 00014 00015 00018 00019 00020 00021 template<typename eT> 00022 inline 00023 void 00024 op_cor::direct_cor(Mat<eT>& out, const Mat<eT>& A, const uword norm_type) 00025 { 00026 arma_extra_debug_sigprint(); 00027 00028 if(A.is_empty()) 00029 { 00030 out.reset(); 00031 return; 00032 } 00033 00034 if(A.is_vec()) 00035 { 00036 out.set_size(1,1); 00037 out[0] = eT(1); 00038 } 00039 else 00040 { 00041 const uword N = A.n_rows; 00042 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); 00043 00044 const Row<eT> acc = sum(A); 00045 const Row<eT> sd = stddev(A); 00046 00047 out = (trans(A) * A); 00048 out -= (trans(acc) * acc)/eT(N); 00049 out /= norm_val; 00050 out /= trans(sd) * sd; 00051 } 00052 } 00053 00054 00055 00056 template<typename T> 00057 inline 00058 void 00059 op_cor::direct_cor(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const uword norm_type) 00060 { 00061 arma_extra_debug_sigprint(); 00062 00063 typedef typename std::complex<T> eT; 00064 00065 if(A.is_empty()) 00066 { 00067 out.reset(); 00068 return; 00069 } 00070 00071 if(A.is_vec()) 00072 { 00073 out.set_size(1,1); 00074 out[0] = eT(1); 00075 } 00076 else 00077 { 00078 const uword N = A.n_rows; 00079 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); 00080 00081 const Row<eT> acc = sum(A); 00082 const Row<T> sd = stddev(A); 00083 00084 out = trans(A) * A; // out = strans(conj(A)) * A; 00085 out -= (trans(acc) * acc)/eT(N); // out -= (strans(conj(acc)) * acc)/eT(N); 00086 out /= norm_val; 00087 00088 //out = out / (trans(sd) * sd); 00089 out /= conv_to< Mat<eT> >::from(trans(sd) * sd); 00090 } 00091 } 00092 00093 00094 00095 template<typename T1> 00096 inline 00097 void 00098 op_cor::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cor>& in) 00099 { 00100 arma_extra_debug_sigprint(); 00101 00102 typedef typename T1::elem_type eT; 00103 00104 const unwrap_check<T1> tmp(in.m, out); 00105 const Mat<eT>& A = tmp.M; 00106 00107 const uword norm_type = in.aux_uword_a; 00108 00109 op_cor::direct_cor(out, A, norm_type); 00110 } 00111 00112 00113