glue_cor_meat.hpp
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 template<typename eT>
00024 inline
00025 void
00026 glue_cor::direct_cor(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const u32 norm_type)
00027 {
00028 arma_extra_debug_sigprint();
00029
00030 if(A.is_vec() && B.is_vec())
00031 {
00032 arma_debug_check( (A.n_elem != B.n_elem), "cor(): the number of elements in A and B must match" );
00033
00034 const eT* A_ptr = A.memptr();
00035 const eT* B_ptr = B.memptr();
00036
00037 eT A_acc = eT(0);
00038 eT B_acc = eT(0);
00039 eT out_acc = eT(0);
00040
00041 const u32 N = A.n_elem;
00042
00043 for(u32 i=0; i<N; ++i)
00044 {
00045 const eT A_tmp = A_ptr[i];
00046 const eT B_tmp = B_ptr[i];
00047
00048 A_acc += A_tmp;
00049 B_acc += B_tmp;
00050
00051 out_acc += A_tmp * B_tmp;
00052 }
00053
00054 out_acc -= (A_acc * B_acc)/eT(N);
00055
00056 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
00057
00058 out.set_size(1,1);
00059 out[0] = out_acc/norm_val;
00060
00061 const Mat<eT> stddev_A = (A.n_rows == 1) ? stddev(trans(A)) : stddev(A);
00062 const Mat<eT> stddev_B = (B.n_rows == 1) ? stddev(trans(B)) : stddev(B);
00063
00064 out /= stddev_A * stddev_B;
00065 }
00066 else
00067 {
00068 arma_debug_assert_same_size(A, B, "cor()");
00069
00070 const u32 N = A.n_rows;
00071 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
00072
00073 out = trans(A) * B;
00074 out -= (trans(sum(A)) * sum(B))/eT(N);
00075 out /= norm_val;
00076 out /= trans(stddev(A)) * stddev(B);
00077 }
00078 }
00079
00080
00081
00082 template<typename T>
00083 inline
00084 void
00085 glue_cor::direct_cor(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const Mat< std::complex<T> >& B, const u32 norm_type)
00086 {
00087 arma_extra_debug_sigprint();
00088
00089 typedef typename std::complex<T> eT;
00090
00091 if(A.is_vec() && B.is_vec())
00092 {
00093 arma_debug_check( (A.n_elem != B.n_elem), "cor(): the number of elements in A and B must match" );
00094
00095 const eT* A_ptr = A.memptr();
00096 const eT* B_ptr = B.memptr();
00097
00098 eT A_acc = eT(0);
00099 eT B_acc = eT(0);
00100 eT out_acc = eT(0);
00101
00102 const u32 N = A.n_elem;
00103
00104 for(u32 i=0; i<N; ++i)
00105 {
00106 const eT A_tmp = A_ptr[i];
00107 const eT B_tmp = B_ptr[i];
00108
00109 A_acc += A_tmp;
00110 B_acc += B_tmp;
00111
00112 out_acc += std::conj(A_tmp) * B_tmp;
00113 }
00114
00115 out_acc -= (std::conj(A_acc) * B_acc)/eT(N);
00116
00117 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
00118
00119 out.set_size(1,1);
00120 out[0] = out_acc/norm_val;
00121
00122 const Mat<T> stddev_A = (A.n_rows == 1) ? stddev(trans(A)) : stddev(A);
00123 const Mat<T> stddev_B = (B.n_rows == 1) ? stddev(trans(B)) : stddev(B);
00124
00125 out /= conv_to< Mat<eT> >::from( stddev_A * stddev_B );
00126 }
00127 else
00128 {
00129 arma_debug_assert_same_size(A, B, "cor()");
00130
00131 const u32 N = A.n_rows;
00132 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
00133
00134 out = trans(conj(A)) * B;
00135 out -= (trans(conj(sum(A))) * sum(B))/eT(N);
00136 out /= norm_val;
00137 out /= conv_to< Mat<eT> >::from( trans(stddev(A)) * stddev(B) );
00138 }
00139 }
00140
00141
00142
00143 template<typename T1, typename T2>
00144 inline
00145 void
00146 glue_cor::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_cor>& X)
00147 {
00148 arma_extra_debug_sigprint();
00149
00150 typedef typename T1::elem_type eT;
00151
00152 const unwrap_check<T1> A_tmp(X.A, out);
00153 const unwrap_check<T2> B_tmp(X.B, out);
00154
00155 const Mat<eT>& A = A_tmp.M;
00156 const Mat<eT>& B = B_tmp.M;
00157
00158 const u32 norm_type = X.aux_u32;
00159
00160 if(&A != &B)
00161 {
00162 glue_cor::direct_cor(out, A, B, norm_type);
00163 }
00164 else
00165 {
00166 op_cor::direct_cor(out, A, norm_type);
00167 }
00168
00169 }
00170
00171
00172
00173