annotate armadillo-3.900.4/include/armadillo_bits/op_cor_meat.hpp @ 84:55a047986812 tip

Update library URI so as not to be document-local
author Chris Cannam
date Wed, 22 Apr 2020 14:21:57 +0100
parents 1ec0e2823891
children
rev   line source
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 //! @}