fml  0.1-0
Fused Matrix Library
cor.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_GPU_STATS_COR_H
6 #define FML_GPU_STATS_COR_H
7 #pragma once
8 
9 
10 #include <stdexcept>
11 
12 #include "../../_internals/omp.hh"
13 
14 #include "../gpumat.hh"
15 #include "../gpuvec.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, true, 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, true, x);
61  dimops::scale(true, true, y);
62 
63  const REAL alpha = 1. / ((REAL) (x.nrows()-1));
64  linalg::matmult(true, false, alpha, x, y, cov);
65  }
66 
67 
68 
69  namespace internals
70  {
71  template <typename REAL>
72  __global__ void kernel_cor_vecvec(const len_t n, const REAL *x,
73  const REAL *y, const REAL *meanx, const REAL *meany, REAL *cp, REAL *normx,
74  REAL *normy)
75  {
76  int i = blockDim.x*blockIdx.x + threadIdx.x;
77  int j = blockDim.y*blockIdx.y + threadIdx.y;
78 
79  if (i < n && j < n)
80  {
81  const REAL xi_mm = x[i] - (*meanx);
82  const REAL yi_mm = y[i] - (*meany);
83 
84  atomicAdd(cp, xi_mm * yi_mm);
85  atomicAdd(normx, xi_mm * xi_mm);
86  atomicAdd(normy, yi_mm * yi_mm);
87  }
88  }
89  }
90 
92  template <typename REAL>
93  REAL cor(const gpuvec<REAL> &x, const gpuvec<REAL> &y)
94  {
95  const len_t n = x.size();
96  if (n != y.size())
97  throw std::runtime_error("non-conformal arguments");
98 
99  const REAL meanx = x.sum() / n;
100  const REAL meany = y.sum() / n;
101  gpuscalar<REAL> meanx_gpu(x.get_card(), meanx);
102  gpuscalar<REAL> meany_gpu(x.get_card(), meany);
103 
104  REAL cp = 0, normx = 0, normy = 0;
105  gpuscalar<REAL> cp_gpu(x.get_card(), cp);
106  gpuscalar<REAL> normx_gpu(x.get_card(), normx);
107  gpuscalar<REAL> normy_gpu(x.get_card(), normy);
108 
109  internals::kernel_cor_vecvec<<<x.get_griddim(), x.get_blockdim()>>>(n,
110  x.data_ptr(), y.data_ptr(), meanx_gpu.data_ptr(), meany_gpu.data_ptr(),
111  cp_gpu.data_ptr(), normx_gpu.data_ptr(), normy_gpu.data_ptr());
112 
113  cp_gpu.get_val(&cp);
114  normx_gpu.get_val(&normx);
115  normy_gpu.get_val(&normy);
116 
117  return cp / sqrt(normx * normy);
118  }
119 }
120 }
121 
122 
123 #endif
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::gpuvec
Vector class for data held on a single GPU.
Definition: gpuvec.hh:32
fml::unimat::nrows
len_t nrows() const
Number of rows.
Definition: unimat.hh:36
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::stats::cor
void cor(cpumat< REAL > &x, cpumat< REAL > &cov)
Covariance.
Definition: cor.hh:43
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::gpumat
Matrix class for data held on a single GPU.
Definition: gpumat.hh:35
fml::gpuvec::sum
T sum() const
Sum the vector.
Definition: gpuvec.hh:472
fml::gpuscalar
Definition: gpuscalar.hh:16
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