Matrix.h
1 /* Copyright (C) 2020 IBM Corp.
2  * This program is Licensed under the Apache License, Version 2.0
3  * (the "License"); you may not use this file except in compliance
4  * with the License. You may obtain a copy of the License at
5  * http://www.apache.org/licenses/LICENSE-2.0
6  * Unless required by applicable law or agreed to in writing, software
7  * distributed under the License is distributed on an "AS IS" BASIS,
8  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
9  * See the License for the specific language governing permissions and
10  * limitations under the License. See accompanying LICENSE file.
11  */
12 #ifndef HELIB_MATRIX_H
13 #define HELIB_MATRIX_H
14 
15 #include <iostream>
16 #include <utility> // swap
17 #include <vector>
18 #include <numeric>
19 #include <array>
20 #include <string>
21 #include <algorithm>
22 #include <functional>
23 #include <type_traits>
24 #include <initializer_list>
25 
26 #include <NTL/BasicThreadPool.h>
27 
28 #include "assertions.h"
29 #include "helib.h"
30 #include "zeroValue.h"
31 
32 /*
33  * The code in this file should be considered very alpha.
34  * It is in flux, so not advised for public use, yet.
35  * These APIs differ from matmul; here we think of vectors of Ctxts
36  * not a Ctxt encoding a vector.
37  */
38 
39 namespace helib {
40 
41 template <std::size_t N>
43 {
44  std::array<std::size_t, N> lengths; // the dims.
45  std::array<std::size_t, N> strides;
46  std::vector<long> start = {0};
47  std::size_t size = 0;
48 
49  // Constructor used when taking a 'view'
50  template <typename Iter1, typename Iter2>
51  TensorSlice(const Iter1& firstLength,
52  const Iter1& lastLength,
53  const Iter2& firstStride,
54  const Iter2& lastStride,
55  const std::vector<long>& st) :
56  start(st)
57  {
58  std::copy(firstLength, lastLength, this->lengths.begin());
59  std::copy(firstStride, lastStride, this->strides.begin());
60  this->size = (std::accumulate(lengths.begin(),
61  lengths.end(),
62  1,
63  std::multiplies<std::size_t>()));
64  }
65 
66  // Constructor used when taking a 'view'
67  template <typename Iter1, typename Iter2>
68  TensorSlice(const Iter1& firstLength,
69  const Iter1& lastLength,
70  const Iter2& firstStride,
71  const Iter2& lastStride,
72  unsigned long pos) :
73  start(1, pos)
74  {
75  std::copy(firstLength, lastLength, this->lengths.begin());
76  std::copy(firstStride, lastStride, this->strides.begin());
77  this->size = (std::accumulate(lengths.begin(),
78  lengths.end(),
79  1,
80  std::multiplies<std::size_t>()));
81  }
82 
83  // Constructor for setting up a Tensor.
84  template <typename... Dims>
85  TensorSlice(Dims... dims) :
86  lengths{std::size_t(dims)...},
87  size(std::accumulate(lengths.begin(),
88  lengths.end(),
89  1,
90  std::multiplies<std::size_t>()))
91  {
92  // TODO add helib::assert or static checkhere w.r.t. N
93  this->strides.back() = 1;
94  if (N > 1) {
95  std::copy(lengths.begin() + 1, lengths.end(), strides.begin());
96  }
97  }
98 
99  template <typename... Dims>
100  std::size_t operator()(Dims... dims) const
101  {
102  static_assert(sizeof...(Dims) == N, "Wrong number of indices given.");
103 
104  std::array<std::size_t, N> args{std::size_t(dims)...};
105  // Check indices are within bounds.
106  for (long i = 0; i < long(N); ++i) {
107  if (args[i] >= this->lengths[i]) {
109  "Index given: " + std::to_string(args[i]) +
110  ". Max value is: " + std::to_string(this->lengths[i]));
111  }
112  }
113 
114  if (this->start.size() == 1) {
115  return std::inner_product(args.begin(),
116  args.end(),
117  this->strides.begin(),
118  this->start.at(0));
119  } else {
120  // FIXME - this will only work with Matrices
121  return std::inner_product(args.begin(),
122  args.end(),
123  this->strides.begin(),
124  this->start.at(args.at(1)) * strides.back());
125  }
126  }
127 
128  std::size_t order() const { return N; }
129 
130  bool operator==(const TensorSlice& rhs) const
131  {
132  if (this == &rhs)
133  return true;
134  else if (this->size == rhs.size && this->start == rhs.start &&
135  this->lengths == rhs.lengths && this->strides == rhs.strides)
136  return true;
137  else
138  return false;
139  }
140 
141  bool operator!=(const TensorSlice& rhs) const { return !(*this == rhs); }
142 };
143 
144 // N is the dimension of the tensor
145 // At the moment we care about N = 1 or 2 only.
146 
147 template <typename T, std::size_t N>
148 class Tensor
149 {
150 
151 private:
152  TensorSlice<N> subscripts;
153  std::shared_ptr<std::vector<T>> elements_ptr;
154  bool full_view = true;
155 
156 public:
157  // Construct with all entries equal to obj, often useful if T does not have a
158  // default constructor or if default-constructed Ts are invalid.
159  // This has to be disabled if dims is convertible to size_t for now,
160  // otherwise a call to e.g. Tensor(obj, 5) will use the below version
161  // instead of this. TODO: think of a way around this.
162  template <typename U = T,
163  typename... Dims,
164  typename std::enable_if_t<
165  !std::is_convertible<U, std::size_t>::value>* = nullptr>
166  Tensor(const T& obj, Dims... dims) :
167  subscripts{dims...},
168  elements_ptr(std::make_shared<std::vector<T>>(subscripts.size, obj))
169  {}
170 
171  template <typename... Dims>
172  Tensor(Dims... dims) :
173  subscripts{std::size_t(dims)...},
174  elements_ptr(std::make_shared<std::vector<T>>(subscripts.size))
175  {}
176 
177  // TODO generalise this
178  Tensor(std::initializer_list<std::vector<T>> lst) :
179  subscripts{lst.size(), lst.begin()->size()},
180  elements_ptr(
181  std::make_shared<std::vector<T>>(lst.size() * lst.begin()->size()))
182  {
183  int column_length = lst.begin()->size();
184  int cnt = 0;
185  for (const auto& v : lst) {
186  if (column_length != long(v.size()))
187  throw helib::LogicError(
188  "Column dimensions do not match on initializer list.");
189 
190  std::copy(v.begin(),
191  v.end(),
192  this->elements_ptr->begin() + (column_length * cnt++));
193  }
194  }
195 
197  const std::shared_ptr<std::vector<T>>& elems) :
198  subscripts(ts), elements_ptr(elems), full_view(false)
199  {}
200 
201  Tensor(const Tensor& other) = default;
202  Tensor(Tensor&& other) = default;
203 
204  Tensor& operator=(const Tensor& rhs) = default;
205  Tensor& operator=(Tensor&& rhs) = default;
206 
207  ~Tensor() = default;
208 
209  std::size_t order() const { return N; }
210 
211  template <typename... Args>
212  T& operator()(Args... args)
213  {
214  return this->elements_ptr->at(subscripts(args...));
215  }
216 
217  template <typename... Args>
218  const T& operator()(Args... args) const
219  {
220  return this->elements_ptr->at(subscripts(args...));
221  }
222 
223  std::size_t size() const { return this->subscripts.size; }
224 
225  std::size_t dims(int i) const { return this->subscripts.lengths.at(i); }
226 
227  bool fullView() const { return this->full_view; }
228 
229  bool operator==(const Tensor& rhs) const
230  {
231  if (this == &rhs) {
232  return true;
233  } else if (this->subscripts != rhs.subscripts) {
234  return false;
235  } else if (this->fullView() && rhs.fullView() &&
236  *elements_ptr == *(rhs.elements_ptr)) {
237  return true;
238  } else {
239  for (size_t i = 0; i < dims(0); ++i)
240  for (size_t j = 0; j < dims(1); ++j)
241  if (this->operator()(i, j) != rhs(i, j)) {
242  return false;
243  }
244  return true;
245  }
246  }
247 
248  bool operator!=(const Tensor& rhs) const { return !(*this == rhs); }
249 
250  Tensor<T, N - 1> row(std::size_t i) const
251  {
252  TensorSlice<N - 1> ts(this->subscripts.lengths.begin() + 1,
253  this->subscripts.lengths.end(),
254  this->subscripts.strides.begin() + 1,
255  this->subscripts.strides.end(),
256  /*row X strides*/ i * this->subscripts.strides.at(0));
257  return Tensor<T, N - 1>(ts, this->elements_ptr);
258  }
259 
260  Tensor<T, N - 1> column(std::size_t j) const
261  {
262  // TODO - can only exist with N > 1
263  TensorSlice<N - 1> ts(this->subscripts.lengths.begin(),
264  this->subscripts.lengths.end() - 1,
265  this->subscripts.strides.begin(),
266  this->subscripts.strides.end() - 1,
267  /*col*/ j);
268  return Tensor<T, N - 1>(ts, this->elements_ptr);
269  }
270 
271  // FIXME - returns several columns, same order as this. Only Matrices.
272  Tensor<T, N> columns(const std::vector<long>& js) const
273  {
274  // Check the js are valid
275  for (const auto& j : js)
276  assertInRange<LogicError>(
277  j,
278  0l,
279  static_cast<long>(this->dims(1)),
280  "Index for column does not exist. Given index " + std::to_string(j) +
281  ". Expected index in " + "range [0, " +
282  std::to_string(this->dims(1)) + ").");
283 
284  // FIXME: lengths - works only for matrices at mo.
285  std::vector<std::size_t> lengths = {this->dims(0), js.size()};
286 
287  std::vector<long> offsets(js);
288  for (std::size_t i = 0; i < offsets.size(); ++i) {
289  offsets[i] -= i;
290  }
291 
292  TensorSlice<N> ts(lengths.begin(),
293  lengths.end(),
294  this->subscripts.strides.begin(),
295  this->subscripts.strides.end(),
296  offsets);
297 
298  return Tensor<T, N>(ts, this->elements_ptr);
299  }
300 
301  template <typename T2>
303  std::function<T&(T&, const T2&)> operation)
304  {
305  // rhs is not of the same type thus need the dims.
306  std::array<std::size_t, N> rhs_subscripts;
307  for (std::size_t i = 0; i < N; ++i) {
308  rhs_subscripts[i] = rhs.dims(i);
309  }
310 
311  // Sanity Check: Are they the same dimensions?
312  if (!std::equal(this->subscripts.lengths.begin(),
313  this->subscripts.lengths.end(),
314  rhs_subscripts.begin())) {
315  throw helib::LogicError("Matrix dimensions do not match.");
316  }
317 
318  // TODO For now, we do not allow views to operate on views to the same data.
319  if (static_cast<const void*>(&this->data()) ==
320  static_cast<const void*>(&rhs.data())) {
321  throw helib::LogicError("Views point to same underlying data.");
322  }
323 
324  // Optimisation if they have full view of underlying memmory.
325  if (this->full_view && rhs.fullView()) {
326  const std::vector<T2>& rhs_v = rhs.data();
327  for (std::size_t i = 0; i < this->elements_ptr->size(); ++i) {
328  operation((*this->elements_ptr)[i], rhs_v[i]);
329  }
330  } else {
331  // TODO - again will only work for Matrices.
332  for (std::size_t j = 0; j < this->dims(1); ++j) {
333  for (std::size_t i = 0; i < this->dims(0); ++i) {
334  operation(this->operator()(i, j), rhs(i, j));
335  }
336  }
337  }
338 
339  return *this;
340  }
341 
342  template <typename T2>
344  {
345  return entrywiseOperation<T2>(
346  rhs,
347  [](auto& lhs, const auto& rhs) -> decltype(auto) {
348  return lhs += rhs;
349  });
350  }
351 
352  template <typename T2>
354  {
355  return entrywiseOperation<T2>(
356  rhs,
357  [](auto& lhs, const auto& rhs) -> decltype(auto) {
358  return lhs -= rhs;
359  });
360  }
361 
362  template <typename T2>
364  {
365  return entrywiseOperation<T2>(
366  rhs,
367  [](auto& lhs, const auto& rhs) -> decltype(auto) {
368  return lhs *= rhs;
369  });
370  }
371 
372  Tensor<T, N>& apply(std::function<void(T& x)> fn)
373  {
374  // Optimisation if they have full view of underlying memory.
375  if (this->full_view) {
376  NTL_EXEC_RANGE(long(this->elements_ptr->size()), first, last)
377  for (long i = first; i < last; ++i)
378  fn((*elements_ptr)[i]);
379  NTL_EXEC_RANGE_END
380 
381  } else {
382  // TODO - again will only work for Matrices.
383  NTL_EXEC_RANGE(this->dims(1), first, last)
384  for (long j = first; j < last; ++j)
385  for (std::size_t i = 0; i < this->dims(0); ++i)
386  fn(this->operator()(i, j));
387  NTL_EXEC_RANGE_END
388  }
389  return *this;
390  }
391 
392  // Matrix special
394  {
395  if (fullView()) {
396  // In this case we can perform a transpose copy-free using swaps
397  // First express the desired permutation as a list of numbers
398  std::vector<int> permutation(size());
399  std::iota(permutation.begin(), permutation.end(), 0);
400  for (int& num : permutation)
401  num = (num % subscripts.lengths[1]) * subscripts.lengths[0] +
402  num / subscripts.lengths[1];
403  // Now write the permutation as a list of disjoint cycles
404  std::vector<std::vector<int>> cycles;
405  std::vector<bool> seen(size(), false);
406  int num_processed = 0;
407  int current_pos = 0;
408  while (num_processed < long(size())) {
409  // Build the cycle starting at current_pos
410  std::vector<int> cycle = {current_pos};
411  seen[current_pos] = true;
412  while (permutation.at(cycle.back()) != cycle.front()) {
413  seen[permutation.at(cycle.back())] = true;
414  cycle.push_back(permutation.at(cycle.back()));
415  }
416  num_processed += cycle.size();
417  cycles.push_back(std::move(cycle));
418  // Move current_pos to the next non-seen value
419  while (current_pos < long(size()) && seen[current_pos])
420  ++current_pos;
421  }
422  // Now turn this product of disjoint cycles into a list of pairs
423  std::vector<std::pair<int, int>> swaps;
424  for (const auto& cycle : cycles)
425  if (cycle.size() >= 2)
426  for (int i = cycle.size() - 1; i > 0; --i)
427  swaps.emplace_back(cycle[i], cycle[i - 1]);
428  // Now do the swaps
429  for (const auto& swap : swaps)
430  std::swap(elements_ptr->at(swap.first), elements_ptr->at(swap.second));
431  } else {
432  auto new_elements =
433  std::make_shared<std::vector<T>>(size(), data().front());
434  int count = size();
435  int i = 0;
436  int j = 0;
437  int new_i = 0;
438  int new_j = 0;
439  while (count--) {
440  // Put the new element in the correct place
441  new_elements->at(new_i * subscripts.lengths[0] + new_j) = operator()(i,
442  j);
443 
444  // Increment everything.
445  // (i,j) wraps around as per the current configuration.
446  j = (j + 1) % subscripts.lengths[1];
447  if (j == 0)
448  i = (i + 1) % subscripts.lengths[0];
449 
450  // (new_i,new_j) increments and wraps around as per the new, transposed
451  // configuration.
452  new_i = (new_i + 1) % subscripts.lengths[1];
453  if (new_i == 0)
454  new_j = (new_j + 1) % subscripts.lengths[0];
455  }
456  elements_ptr = new_elements;
457  full_view = true;
458  }
459  std::reverse(subscripts.lengths.begin(), subscripts.lengths.end());
460  subscripts.strides.front() = subscripts.lengths.back();
461  subscripts.start = {0};
462  return *this;
463  }
464 
465  const std::vector<T>& data() const { return *this->elements_ptr; }
466 };
467 
468 // Matrix special - Different types
469 template <typename T,
470  typename T2,
471  typename std::enable_if_t<
472  std::is_convertible<T, std::size_t>::value>* = nullptr>
473 inline Tensor<T, 2> operator*(const Tensor<T, 2>& M1, const Tensor<T2, 2>& M2)
474 {
475  HELIB_NTIMER_START(MatrixMultiplicationConv);
476  // rows from M1 x cols from M2
477  Tensor<T, 2> R(M1.dims(0), M2.dims(1));
478 
479  // Sanity check: M1 cols should match M2 rows
480  if (M1.dims(1) != M2.dims(0)) {
481  throw helib::LogicError("Matrix inner dimensions do not match.");
482  }
483 
484  // TODO add NTL thread pool.
485  // TODO swap += for some binary sum.
486  NTL_EXEC_RANGE(M1.dims(0), first, last)
487  for (long i = first; i < last; ++i)
488  // for (std::size_t i = 0; i < M1.dims(0); ++i)
489  for (std::size_t j = 0; j < M2.dims(1); ++j)
490  for (std::size_t k = 0; k < M2.dims(0); ++k) {
491  T R_tmp = M1(i, k);
492  R_tmp *= M2(k, j);
493  R(i, j) += R_tmp;
494  }
495  NTL_EXEC_RANGE_END
496 
497  HELIB_NTIMER_STOP(MatrixMultiplicationConv);
498  return R;
499 }
500 
501 // Matrix special - Different types
502 template <typename T,
503  typename T2,
504  typename std::enable_if_t<
505  !std::is_convertible<T, std::size_t>::value>* = nullptr>
506 inline Tensor<T, 2> operator*(const Tensor<T, 2>& M1, const Tensor<T2, 2>& M2)
507 {
508  HELIB_NTIMER_START(MatrixMultiplicationNotConv);
509  // rows from M1 x cols from M2
510  Tensor<T, 2> R(helib::zeroValue(M1(0, 0)), M1.dims(0), M2.dims(1));
511 
512  // Sanity check: M1 cols should match M2 rows
513  if (M1.dims(1) != M2.dims(0)) {
514  throw helib::LogicError(
515  "The number of columns in left matrix (" + std::to_string(M1.dims(1)) +
516  ") do not match the number of rows of the right matrix (" +
517  std::to_string(M2.dims(0)) + ").");
518  }
519 
520  // TODO swap += for some binary sum.
521  NTL_EXEC_RANGE(M1.dims(0), first, last)
522  for (long i = first; i < last; ++i)
523  // For reference before NTL threads: for (std::size_t i = 0; i < M1.dims(0);
524  // ++i)
525  for (std::size_t j = 0; j < M2.dims(1); ++j)
526  for (std::size_t k = 0; k < M2.dims(0); ++k) {
527  T R_tmp = M1(i, k);
528  R_tmp *= M2(k, j);
529  R(i, j) += R_tmp;
530  }
531  NTL_EXEC_RANGE_END
532 
533  HELIB_NTIMER_STOP(MatrixMultiplicationNotConv);
534  return R;
535 }
536 
537 // These are the alias we will actually use.
538 template <typename T>
540 
541 template <typename T>
543 
544 // TODO Columns Tuple
545 template <typename T>
547 {
548 
549 private:
550  std::vector<Matrix<T>> columns;
551 
552 public:
553  MatrixView(const std::initializer_list<Matrix<T>> lst) : columns(lst) {}
554 };
555 
556 template <typename T>
557 void printMatrix(const Matrix<T>& M, std::ostream& out = std::cout)
558 {
559  for (std::size_t i = 0; i < M.dims(0); ++i) {
560  for (std::size_t j = 0; j < M.dims(1); ++j)
561  out << M(i, j) << " ";
562  out << "\n";
563  }
564 }
565 
566 } // namespace helib
567 
568 #endif // ifndef HELIB_MATRIX_H
bool operator==(const TensorSlice &rhs) const
Definition: Matrix.h:130
~Tensor()=default
Tensor(Dims... dims)
Definition: Matrix.h:172
Tensor< T, N > & hadamard(const Tensor< T2, N > &rhs)
Definition: Matrix.h:363
std::size_t size() const
Definition: Matrix.h:223
TensorSlice(const Iter1 &firstLength, const Iter1 &lastLength, const Iter2 &firstStride, const Iter2 &lastStride, unsigned long pos)
Definition: Matrix.h:68
Tensor< T, 2 > operator*(const Tensor< T, 2 > &M1, const Tensor< T2, 2 > &M2)
Definition: Matrix.h:473
std::size_t size
Definition: Matrix.h:47
Inherits from Exception and std::logic_error.
Definition: exceptions.h:68
std::size_t dims(int i) const
Definition: Matrix.h:225
bool operator!=(const TensorSlice &rhs) const
Definition: Matrix.h:141
Tensor< T, N > & operator-=(const Tensor< T2, N > &rhs)
Definition: Matrix.h:353
const T & operator()(Args... args) const
Definition: Matrix.h:218
std::vector< long > start
Definition: Matrix.h:46
Definition: Matrix.h:547
std::array< std::size_t, N > strides
Definition: Matrix.h:45
std::array< std::size_t, N > lengths
Definition: Matrix.h:44
Tensor(const T &obj, Dims... dims)
Definition: Matrix.h:166
T zeroValue(const T &x)
Given an object x return a zero object of the same type.
Definition: zeroValue.h:29
Tensor(const TensorSlice< N > &ts, const std::shared_ptr< std::vector< T >> &elems)
Definition: Matrix.h:196
TensorSlice(const Iter1 &firstLength, const Iter1 &lastLength, const Iter2 &firstStride, const Iter2 &lastStride, const std::vector< long > &st)
Definition: Matrix.h:51
TensorSlice(Dims... dims)
Definition: Matrix.h:85
bool fullView() const
Definition: Matrix.h:227
Definition: Matrix.h:43
MatrixView(const std::initializer_list< Matrix< T >> lst)
Definition: Matrix.h:553
Tensor< T, N > & entrywiseOperation(const Tensor< T2, N > &rhs, std::function< T &(T &, const T2 &)> operation)
Definition: Matrix.h:302
Definition: apiAttributes.h:21
std::size_t operator()(Dims... dims) const
Definition: Matrix.h:100
const std::vector< T > & data() const
Definition: Matrix.h:465
Tensor< T, N - 1 > column(std::size_t j) const
Definition: Matrix.h:260
Tensor< T, 2 > & transpose()
Definition: Matrix.h:393
Definition: Matrix.h:149
std::size_t order() const
Definition: Matrix.h:128
bool operator!=(const Tensor &rhs) const
Definition: Matrix.h:248
Tensor(std::initializer_list< std::vector< T >> lst)
Definition: Matrix.h:178
Tensor & operator=(const Tensor &rhs)=default
Tensor< T, N - 1 > row(std::size_t i) const
Definition: Matrix.h:250
Inherits from Exception and std::out_of_range.
Definition: exceptions.h:86
Tensor(Tensor &&other)=default
bool operator==(const Tensor &rhs) const
Definition: Matrix.h:229
void swap(cloned_ptr< X, Cloner > &x, cloned_ptr< X, Cloner > &y)
Definition: clonedPtr.h:170
void printMatrix(const Matrix< T > &M, std::ostream &out=std::cout)
Definition: Matrix.h:557
Tensor & operator=(Tensor &&rhs)=default
Tensor(const Tensor &other)=default
Tensor< T, N > columns(const std::vector< long > &js) const
Definition: Matrix.h:272
std::size_t order() const
Definition: Matrix.h:209
Tensor< T, N > & operator+=(const Tensor< T2, N > &rhs)
Definition: Matrix.h:343
Tensor< T, N > & apply(std::function< void(T &x)> fn)
Definition: Matrix.h:372
T & operator()(Args... args)
Definition: Matrix.h:212