Chris@49
|
1 // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
|
Chris@49
|
2 // Copyright (C) 2009-2011 Conrad Sanderson
|
Chris@49
|
3 // Copyright (C) 2009-2010 Dimitrios Bouzas
|
Chris@49
|
4 //
|
Chris@49
|
5 // This Source Code Form is subject to the terms of the Mozilla Public
|
Chris@49
|
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
|
Chris@49
|
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
Chris@49
|
8
|
Chris@49
|
9
|
Chris@49
|
10
|
Chris@49
|
11 //! \addtogroup op_cor
|
Chris@49
|
12 //! @{
|
Chris@49
|
13
|
Chris@49
|
14
|
Chris@49
|
15
|
Chris@49
|
16 template<typename eT>
|
Chris@49
|
17 inline
|
Chris@49
|
18 void
|
Chris@49
|
19 op_cor::direct_cor(Mat<eT>& out, const Mat<eT>& A, const uword norm_type)
|
Chris@49
|
20 {
|
Chris@49
|
21 arma_extra_debug_sigprint();
|
Chris@49
|
22
|
Chris@49
|
23 if(A.is_empty())
|
Chris@49
|
24 {
|
Chris@49
|
25 out.reset();
|
Chris@49
|
26 return;
|
Chris@49
|
27 }
|
Chris@49
|
28
|
Chris@49
|
29 if(A.is_vec())
|
Chris@49
|
30 {
|
Chris@49
|
31 out.set_size(1,1);
|
Chris@49
|
32 out[0] = eT(1);
|
Chris@49
|
33 }
|
Chris@49
|
34 else
|
Chris@49
|
35 {
|
Chris@49
|
36 const uword N = A.n_rows;
|
Chris@49
|
37 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
|
Chris@49
|
38
|
Chris@49
|
39 const Row<eT> acc = sum(A);
|
Chris@49
|
40 const Row<eT> sd = stddev(A);
|
Chris@49
|
41
|
Chris@49
|
42 out = (trans(A) * A);
|
Chris@49
|
43 out -= (trans(acc) * acc)/eT(N);
|
Chris@49
|
44 out /= norm_val;
|
Chris@49
|
45 out /= trans(sd) * sd;
|
Chris@49
|
46 }
|
Chris@49
|
47 }
|
Chris@49
|
48
|
Chris@49
|
49
|
Chris@49
|
50
|
Chris@49
|
51 template<typename T>
|
Chris@49
|
52 inline
|
Chris@49
|
53 void
|
Chris@49
|
54 op_cor::direct_cor(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const uword norm_type)
|
Chris@49
|
55 {
|
Chris@49
|
56 arma_extra_debug_sigprint();
|
Chris@49
|
57
|
Chris@49
|
58 typedef typename std::complex<T> eT;
|
Chris@49
|
59
|
Chris@49
|
60 if(A.is_empty())
|
Chris@49
|
61 {
|
Chris@49
|
62 out.reset();
|
Chris@49
|
63 return;
|
Chris@49
|
64 }
|
Chris@49
|
65
|
Chris@49
|
66 if(A.is_vec())
|
Chris@49
|
67 {
|
Chris@49
|
68 out.set_size(1,1);
|
Chris@49
|
69 out[0] = eT(1);
|
Chris@49
|
70 }
|
Chris@49
|
71 else
|
Chris@49
|
72 {
|
Chris@49
|
73 const uword N = A.n_rows;
|
Chris@49
|
74 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
|
Chris@49
|
75
|
Chris@49
|
76 const Row<eT> acc = sum(A);
|
Chris@49
|
77 const Row<T> sd = stddev(A);
|
Chris@49
|
78
|
Chris@49
|
79 out = trans(A) * A; // out = strans(conj(A)) * A;
|
Chris@49
|
80 out -= (trans(acc) * acc)/eT(N); // out -= (strans(conj(acc)) * acc)/eT(N);
|
Chris@49
|
81 out /= norm_val;
|
Chris@49
|
82
|
Chris@49
|
83 //out = out / (trans(sd) * sd);
|
Chris@49
|
84 out /= conv_to< Mat<eT> >::from(trans(sd) * sd);
|
Chris@49
|
85 }
|
Chris@49
|
86 }
|
Chris@49
|
87
|
Chris@49
|
88
|
Chris@49
|
89
|
Chris@49
|
90 template<typename T1>
|
Chris@49
|
91 inline
|
Chris@49
|
92 void
|
Chris@49
|
93 op_cor::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cor>& in)
|
Chris@49
|
94 {
|
Chris@49
|
95 arma_extra_debug_sigprint();
|
Chris@49
|
96
|
Chris@49
|
97 typedef typename T1::elem_type eT;
|
Chris@49
|
98
|
Chris@49
|
99 const unwrap_check<T1> tmp(in.m, out);
|
Chris@49
|
100 const Mat<eT>& A = tmp.M;
|
Chris@49
|
101
|
Chris@49
|
102 const uword norm_type = in.aux_uword_a;
|
Chris@49
|
103
|
Chris@49
|
104 op_cor::direct_cor(out, A, norm_type);
|
Chris@49
|
105 }
|
Chris@49
|
106
|
Chris@49
|
107
|
Chris@49
|
108
|
Chris@49
|
109 //! @}
|