fml  0.1-0
Fused Matrix Library
dimops.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_DIMOPS_H
6 #define FML_CPU_DIMOPS_H
7 #pragma once
8 
9 
10 #include <cmath>
11 
12 #include "cpumat.hh"
13 #include "cpuvec.hh"
14 
15 
17 namespace dimops
18 {
19  namespace internals
20  {
21  template <typename REAL>
22  static inline void col_sum(const len_t j, const len_t m, const REAL *x, REAL &sum)
23  {
24  sum = 0;
25 
26  #pragma omp simd reduction(+:sum)
27  for (len_t i=0; i<m; i++)
28  sum += x[i + m*j];
29  }
30 
31  template <typename REAL>
32  static inline void col_mean(const len_t j, const len_t m, const REAL *x, REAL &mean)
33  {
34  mean = 0;
35  col_sum(j, m, x, mean);
36  mean /= (REAL) m;
37  }
38 
39  template <typename REAL>
40  static inline void col_mean_and_var(const len_t j, const len_t m, const REAL *x, REAL &mean, REAL &var)
41  {
42  mean = 0;
43  var = 0;
44 
45  for (len_t i=0; i<m; i++)
46  {
47  REAL dt = x[i + m*j] - mean;
48  mean += dt/((REAL) i+1);
49  var += dt * (x[i + m*j] - mean);
50  }
51 
52  var = sqrt(var / ((REAL) m-1));
53  }
54  }
55 
56 
57 
66  template <typename REAL>
67  void rowsums(const cpumat<REAL> &x, cpuvec<REAL> &s)
68  {
69  const REAL *x_d = x.data_ptr();
70  const len_t m = x.nrows();
71  const len_t n = x.ncols();
72 
73  s.resize(m);
74  s.fill_zero();
75  REAL *s_d = s.data_ptr();
76 
77  for (len_t j=0; j<n; j++)
78  {
79  #pragma omp for simd
80  for (len_t i=0; i<m; i++)
81  s_d[i] += x_d[i + m*j];
82  }
83  }
84 
85 
86 
95  template <typename REAL>
96  void rowmeans(const cpumat<REAL> &x, cpuvec<REAL> &s)
97  {
98  rowsums(x, s);
99 
100  const len_t m = x.nrows();
101  const len_t n = x.ncols();
102  REAL *s_d = s.data_ptr();
103 
104  #pragma omp for simd
105  for (len_t i=0; i<m; i++)
106  s_d[i] /= (REAL) n;
107  }
108 
109 
110 
119  template <typename REAL>
120  void colsums(const cpumat<REAL> &x, cpuvec<REAL> &s)
121  {
122  const REAL *x_d = x.data_ptr();
123  const len_t m = x.nrows();
124  const len_t n = x.ncols();
125 
126  s.resize(n);
127  s.fill_zero();
128  REAL *s_d = s.data_ptr();
129 
130  for (len_t j=0; j<n; j++)
131  internals::col_sum(j, m, x_d, s_d[j]);
132  }
133 
134 
135 
144  template <typename REAL>
145  void colmeans(const cpumat<REAL> &x, cpuvec<REAL> &s)
146  {
147  const REAL *x_d = x.data_ptr();
148  const len_t m = x.nrows();
149  const len_t n = x.ncols();
150 
151  s.resize(n);
152  s.fill_zero();
153  REAL *s_d = s.data_ptr();
154 
155  for (len_t j=0; j<n; j++)
156  internals::col_mean(j, m, x_d, s_d[j]);
157  }
158 
159 
160 
161  namespace internals
162  {
163  template <typename REAL>
164  static inline void center(cpumat<REAL> &x)
165  {
166  REAL *x_d = x.data_ptr();
167  const len_t m = x.nrows();
168  const len_t n = x.ncols();
169 
170  for (len_t j=0; j<n; j++)
171  {
172  REAL mean = 0;
173  internals::col_mean(j, m, x_d, mean);
174 
175  #pragma omp for simd
176  for (len_t i=0; i<m; i++)
177  x_d[i + m*j] -= mean;
178  }
179  }
180 
181  template <typename REAL>
182  static inline void scale(cpumat<REAL> &x)
183  {
184  REAL *x_d = x.data_ptr();
185  const len_t m = x.nrows();
186  const len_t n = x.ncols();
187 
188  for (len_t j=0; j<n; j++)
189  {
190  REAL mean = 0;
191  REAL var = 0;
192  internals::col_mean_and_var(j, m, x_d, mean, var);
193 
194  #pragma omp for simd
195  for (len_t i=0; i<m; i++)
196  x_d[i + m*j] /= var;
197  }
198  }
199 
200  template <typename REAL>
201  static inline void center_and_scale(cpumat<REAL> &x)
202  {
203  REAL *x_d = x.data_ptr();
204  const len_t m = x.nrows();
205  const len_t n = x.ncols();
206 
207  for (len_t j=0; j<n; j++)
208  {
209  REAL mean = 0;
210  REAL var = 0;
211  internals::col_mean_and_var(j, m, x_d, mean, var);
212 
213  #pragma omp for simd
214  for (len_t i=0; i<m; i++)
215  x_d[i + m*j] = (x_d[i + m*j] - mean) / var;
216  }
217  }
218  }
219 
220 
221 
231  template <typename REAL>
232  void scale(const bool rm_mean, const bool rm_sd, cpumat<REAL> &x)
233  {
234  if (rm_mean && rm_sd)
235  internals::center_and_scale(x);
236  else if (rm_mean)
237  internals::center(x);
238  else if (rm_sd)
239  internals::scale(x);
240  }
241 }
242 
243 
244 #endif
void rowsums(const cpumat< REAL > &x, cpuvec< REAL > &s)
Compute the row sums.
Definition: dimops.hh:67
Vector class for data held on a single CPU.
Definition: cpuvec.hh:27
void colmeans(const cpumat< REAL > &x, cpuvec< REAL > &s)
Compute the column means.
Definition: dimops.hh:145
void resize(len_t size)
Resize the internal object storage.
Definition: cpuvec.hh:177
REAL * data_ptr()
Pointer to the internal array.
Definition: unimat.hh:35
T * data_ptr()
Pointer to the internal array.
Definition: univec.hh:26
void rowmeans(const cpumat< REAL > &x, cpuvec< REAL > &s)
Compute the row means.
Definition: dimops.hh:96
len_t nrows() const
Number of rows.
Definition: unimat.hh:31
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:232
void colsums(const cpumat< REAL > &x, cpuvec< REAL > &s)
Compute the column sums.
Definition: dimops.hh:120
void fill_zero()
Set all values to zero.
Definition: cpuvec.hh:283
Row/column operations.
Definition: dimops.hh:17
Matrix class for data held on a single CPU.
Definition: cpumat.hh:32
len_t ncols() const
Number of columns.
Definition: unimat.hh:33