op_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
00024 template<typename eT>
00025 inline
00026 void
00027 op_cor::direct_cor(Mat<eT>& out, const Mat<eT>& A, const u32 norm_type)
00028 {
00029 arma_extra_debug_sigprint();
00030
00031 if(A.is_vec())
00032 {
00033 out.set_size(1,1);
00034 out[0] = eT(1);
00035 }
00036 else
00037 {
00038 const u32 N = A.n_rows;
00039 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
00040
00041 const Row<eT> acc = sum(A);
00042 const Row<eT> sd = stddev(A);
00043
00044 out = (trans(A) * A);
00045 out -= (trans(acc) * acc)/eT(N);
00046 out /= norm_val;
00047 out /= trans(sd) * sd;
00048 }
00049 }
00050
00051
00052
00053 template<typename T>
00054 inline
00055 void
00056 op_cor::direct_cor(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const u32 norm_type)
00057 {
00058 arma_extra_debug_sigprint();
00059
00060 typedef typename std::complex<T> eT;
00061
00062 if(A.is_vec())
00063 {
00064 out.set_size(1,1);
00065 out[0] = eT(1);
00066 }
00067 else
00068 {
00069 const u32 N = A.n_rows;
00070 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
00071
00072 const Row<eT> acc = sum(A);
00073 const Row<T> sd = stddev(A);
00074
00075 out = trans(conj(A)) * A;
00076 out -= (trans(conj(acc)) * acc)/eT(N);
00077 out /= norm_val;
00078
00079
00080 out /= conv_to< Mat<eT> >::from(trans(sd) * sd);
00081 }
00082 }
00083
00084
00085
00086 template<typename T1>
00087 inline
00088 void
00089 op_cor::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cor>& in)
00090 {
00091 arma_extra_debug_sigprint();
00092
00093 typedef typename T1::elem_type eT;
00094
00095 const unwrap_check<T1> tmp(in.m, out);
00096 const Mat<eT>& A = tmp.M;
00097
00098 const u32 norm_type = in.aux_u32_a;
00099
00100 op_cor::direct_cor(out, A, norm_type);
00101 }
00102
00103
00104
00105