$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 00017 00018 00019 00020 template<typename eT> 00021 inline 00022 void 00023 glue_cov::direct_cov(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const uword norm_type) 00024 { 00025 arma_extra_debug_sigprint(); 00026 00027 if(A.is_vec() && B.is_vec()) 00028 { 00029 arma_debug_check( (A.n_elem != B.n_elem), "cov(): the number of elements in A and B must match" ); 00030 00031 const eT* A_ptr = A.memptr(); 00032 const eT* B_ptr = B.memptr(); 00033 00034 eT A_acc = eT(0); 00035 eT B_acc = eT(0); 00036 eT out_acc = eT(0); 00037 00038 const uword N = A.n_elem; 00039 00040 for(uword i=0; i<N; ++i) 00041 { 00042 const eT A_tmp = A_ptr[i]; 00043 const eT B_tmp = B_ptr[i]; 00044 00045 A_acc += A_tmp; 00046 B_acc += B_tmp; 00047 00048 out_acc += A_tmp * B_tmp; 00049 } 00050 00051 out_acc -= (A_acc * B_acc)/eT(N); 00052 00053 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); 00054 00055 out.set_size(1,1); 00056 out[0] = out_acc/norm_val; 00057 } 00058 else 00059 { 00060 arma_debug_assert_same_size(A, B, "cov()"); 00061 00062 const uword N = A.n_rows; 00063 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); 00064 00065 out = trans(A) * B; 00066 out -= (trans(sum(A)) * sum(B))/eT(N); 00067 out /= norm_val; 00068 } 00069 } 00070 00071 00072 00073 template<typename T> 00074 inline 00075 void 00076 glue_cov::direct_cov(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const Mat< std::complex<T> >& B, const uword norm_type) 00077 { 00078 arma_extra_debug_sigprint(); 00079 00080 typedef typename std::complex<T> eT; 00081 00082 if(A.is_vec() && B.is_vec()) 00083 { 00084 arma_debug_check( (A.n_elem != B.n_elem), "cov(): the number of elements in A and B must match" ); 00085 00086 const eT* A_ptr = A.memptr(); 00087 const eT* B_ptr = B.memptr(); 00088 00089 eT A_acc = eT(0); 00090 eT B_acc = eT(0); 00091 eT out_acc = eT(0); 00092 00093 const uword N = A.n_elem; 00094 00095 for(uword i=0; i<N; ++i) 00096 { 00097 const eT A_tmp = A_ptr[i]; 00098 const eT B_tmp = B_ptr[i]; 00099 00100 A_acc += A_tmp; 00101 B_acc += B_tmp; 00102 00103 out_acc += std::conj(A_tmp) * B_tmp; 00104 } 00105 00106 out_acc -= (std::conj(A_acc) * B_acc)/eT(N); 00107 00108 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); 00109 00110 out.set_size(1,1); 00111 out[0] = out_acc/norm_val; 00112 } 00113 else 00114 { 00115 arma_debug_assert_same_size(A, B, "cov()"); 00116 00117 const uword N = A.n_rows; 00118 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); 00119 00120 out = trans(A) * B; // out = strans(conj(A)) * B; 00121 out -= (trans(sum(A)) * sum(B))/eT(N); // out -= (strans(conj(sum(A))) * sum(B))/eT(N); 00122 out /= norm_val; 00123 } 00124 } 00125 00126 00127 00128 template<typename T1, typename T2> 00129 inline 00130 void 00131 glue_cov::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_cov>& X) 00132 { 00133 arma_extra_debug_sigprint(); 00134 00135 typedef typename T1::elem_type eT; 00136 00137 const unwrap_check<T1> A_tmp(X.A, out); 00138 const unwrap_check<T2> B_tmp(X.B, out); 00139 00140 const Mat<eT>& A = A_tmp.M; 00141 const Mat<eT>& B = B_tmp.M; 00142 00143 const uword norm_type = X.aux_uword; 00144 00145 if(&A != &B) 00146 { 00147 glue_cov::direct_cov(out, A, B, norm_type); 00148 } 00149 else 00150 { 00151 op_cov::direct_cov(out, A, norm_type); 00152 } 00153 00154 } 00155 00156 00157