fml  0.1-0
Fused Matrix Library
cov.hh
1 // This file is part of fml which is released under the Boost Software
2 // License, Version 1.0. See accompanying file LICENSE or copy at
3 // https://www.boost.org/LICENSE_1_0.txt
4 
5 #ifndef FML_CPU_STATS_COV_H
6 #define FML_CPU_STATS_COV_H
7 #pragma once
8 
9 
10 #include <stdexcept>
11 
12 #include "../../_internals/omp.hh"
13 
14 #include "../cpumat.hh"
15 #include "../cpuvec.hh"
16 
17 #include "../dimops.hh"
18 
19 #include "../linalg/crossprod.hh"
20 #include "../linalg/matmult.hh"
21 
22 
23 namespace fml
24 {
25 namespace stats
26 {
42  template <typename REAL>
44  {
45  dimops::scale(true, false, x);
46 
47  const REAL alpha = 1. / ((REAL) (x.nrows()-1));
48  linalg::crossprod(alpha, x, cov);
49  }
50 
51 
52 
54  template <typename REAL>
56  {
57  if (x.nrows() != y.nrows())
58  throw std::runtime_error("non-conformable arguments");
59 
60  dimops::scale(true, false, x);
61  dimops::scale(true, false, y);
62 
63  const REAL alpha = 1. / ((REAL) (x.nrows()-1));
64  linalg::matmult(true, false, alpha, x, y, cov);
65  }
66 
67 
68 
70  template <typename REAL>
71  REAL cov(const cpuvec<REAL> &x, const cpuvec<REAL> &y)
72  {
73  const len_t n = x.size();
74  if (n != y.size())
75  throw std::runtime_error("non-conformal arguments");
76 
77  const REAL *x_d = x.data_ptr();
78  const REAL *y_d = y.data_ptr();
79 
80  REAL sum_xy = 0, sum_x = 0, sum_y = 0;
81 
82  #pragma omp parallel for reduction(+: sum_xy, sum_x, sum_y) if(n>fml::omp::OMP_MIN_SIZE)
83  for (len_t i=0; i<n; i++)
84  {
85  const REAL xi = x_d[i];
86  const REAL yi = y_d[i];
87 
88  sum_xy += xi*yi;
89  sum_x += xi;
90  sum_y += yi;
91  }
92 
93  return (sum_xy - (sum_x*sum_y*((REAL) 1./n))) * ((REAL) 1./(n-1));
94  }
95 }
96 }
97 
98 
99 #endif
fml::cpumat
Matrix class for data held on a single CPU.
Definition: cpumat.hh:36
fml::linalg::crossprod
void crossprod(const REAL alpha, const cpumat< REAL > &x, cpumat< REAL > &ret)
Computes lower triangle of alpha*x^T*x.
Definition: crossprod.hh:37
fml::univec::data_ptr
T * data_ptr()
Pointer to the internal array.
Definition: univec.hh:28
fml::unimat::nrows
len_t nrows() const
Number of rows.
Definition: unimat.hh:36
fml::cpuvec
Vector class for data held on a single CPU.
Definition: cpuvec.hh:31
fml::dimops::scale
void scale(const bool rm_mean, const bool rm_sd, cpumat< REAL > &x)
Remove the mean and/or the sd from a matrix.
Definition: dimops.hh:376
fml
Core namespace.
Definition: dimops.hh:10
fml::univec::size
len_t size() const
Number of elements in the vector.
Definition: univec.hh:26
fml::stats::cov
void cov(cpumat< REAL > &x, cpumat< REAL > &cov)
Covariance.
Definition: cov.hh:43
fml::linalg::matmult
void matmult(const bool transx, const bool transy, const REAL alpha, const cpumat< REAL > &x, const cpumat< REAL > &y, cpumat< REAL > &ret)
Computes ret = alpha*op(x)*op(y) where op(A) is A or A^T.
Definition: matmult.hh:43