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_cov
|
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_cov::direct_cov(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_vec())
|
Chris@49
|
24 {
|
Chris@49
|
25 if(A.n_rows == 1)
|
Chris@49
|
26 {
|
Chris@49
|
27 out = var(trans(A), norm_type);
|
Chris@49
|
28 }
|
Chris@49
|
29 else
|
Chris@49
|
30 {
|
Chris@49
|
31 out = var(A, norm_type);
|
Chris@49
|
32 }
|
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
|
Chris@49
|
41 out = trans(A) * A;
|
Chris@49
|
42 out -= (trans(acc) * acc)/eT(N);
|
Chris@49
|
43 out /= norm_val;
|
Chris@49
|
44 }
|
Chris@49
|
45 }
|
Chris@49
|
46
|
Chris@49
|
47
|
Chris@49
|
48
|
Chris@49
|
49 template<typename T>
|
Chris@49
|
50 inline
|
Chris@49
|
51 void
|
Chris@49
|
52 op_cov::direct_cov(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const uword norm_type)
|
Chris@49
|
53 {
|
Chris@49
|
54 arma_extra_debug_sigprint();
|
Chris@49
|
55
|
Chris@49
|
56 typedef typename std::complex<T> eT;
|
Chris@49
|
57
|
Chris@49
|
58 if(A.is_vec())
|
Chris@49
|
59 {
|
Chris@49
|
60 if(A.n_rows == 1)
|
Chris@49
|
61 {
|
Chris@49
|
62 const Mat<T> tmp_mat = var(trans(A), norm_type);
|
Chris@49
|
63 out.set_size(1,1);
|
Chris@49
|
64 out[0] = tmp_mat[0];
|
Chris@49
|
65 }
|
Chris@49
|
66 else
|
Chris@49
|
67 {
|
Chris@49
|
68 const Mat<T> tmp_mat = var(A, norm_type);
|
Chris@49
|
69 out.set_size(1,1);
|
Chris@49
|
70 out[0] = tmp_mat[0];
|
Chris@49
|
71 }
|
Chris@49
|
72 }
|
Chris@49
|
73 else
|
Chris@49
|
74 {
|
Chris@49
|
75 const uword N = A.n_rows;
|
Chris@49
|
76 const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
|
Chris@49
|
77
|
Chris@49
|
78 const Row<eT> acc = sum(A);
|
Chris@49
|
79
|
Chris@49
|
80 out = trans(A) * A; // out = strans(conj(A)) * A;
|
Chris@49
|
81 out -= (trans(acc) * acc)/eT(N); // out -= (strans(conj(acc)) * acc)/eT(N);
|
Chris@49
|
82 out /= norm_val;
|
Chris@49
|
83 }
|
Chris@49
|
84 }
|
Chris@49
|
85
|
Chris@49
|
86
|
Chris@49
|
87
|
Chris@49
|
88 template<typename T1>
|
Chris@49
|
89 inline
|
Chris@49
|
90 void
|
Chris@49
|
91 op_cov::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_cov>& in)
|
Chris@49
|
92 {
|
Chris@49
|
93 arma_extra_debug_sigprint();
|
Chris@49
|
94
|
Chris@49
|
95 typedef typename T1::elem_type eT;
|
Chris@49
|
96
|
Chris@49
|
97 const unwrap_check<T1> tmp(in.m, out);
|
Chris@49
|
98 const Mat<eT>& A = tmp.M;
|
Chris@49
|
99
|
Chris@49
|
100 const uword norm_type = in.aux_uword_a;
|
Chris@49
|
101
|
Chris@49
|
102 op_cov::direct_cov(out, A, norm_type);
|
Chris@49
|
103 }
|
Chris@49
|
104
|
Chris@49
|
105
|
Chris@49
|
106
|
Chris@49
|
107 //! @}
|