hypercube.h

/* Copyright (C) 2012-2020 IBM Corp.
 * This program is Licensed under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance
 * with the License. You may obtain a copy of the License at
 *   http://www.apache.org/licenses/LICENSE-2.0
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License. See accompanying LICENSE file.
 */
#ifndef HELIB_HYPERCUBE_H
#define HELIB_HYPERCUBE_H
 
#include <helib/NumbTh.h>
 
namespace helib {
 
class PAlgebra; // forward declaration
 
class CubeSignature
{
private:
  NTL::Vec<long> dims;  // dims[i] is the size along the i'th dimension
  NTL::Vec<long> prods; // prods[i] = \prod_{j=i}^{n-1} dims[i]
 
public:
  CubeSignature() {} // a nullptr signature
 
  void initSignature(const long _dims[], long _ndims)
  {
 
    // can only initialize a nullptr signature
    assertEq(dims.length(),
             0l,
             "Can only initialize an un-initialized signature");
 
    assertTrue<InvalidArgument>(_ndims >= 0l, "Dimension count is negative");
 
    dims.SetLength(_ndims);
    prods.SetLength(_ndims + 1);
    prods[_ndims] = 1;
    for (long i = _ndims - 1; i >= 0; i--) {
      assertTrue<InvalidArgument>(
          _dims[i] > 0,
          "Bad (non-positive) cube signature dimension");
      dims[i] = _dims[i];
      prods[i] = prods[i + 1] * _dims[i];
    }
  }
  // VecType is either std::vector<intType> or NTL:Vec<intType>
  template <typename VecType>
  void initSignature(const VecType& _dims)
  {
    initSignature(_dims.data(), lsize(_dims));
  }
 
  CubeSignature(const long _dims[], long _ndims)
  {
    initSignature(_dims, _ndims);
  }
  CubeSignature(const NTL::Vec<long>& _dims) { initSignature(_dims); }
  CubeSignature(const std::vector<long>& _dims) { initSignature(_dims); }
 
  /* TODO: We could #include <initializer_list> and then do e.g.,
   * CubeSignature s {1,2,3};
  CubeSignature(std::initializer_list<long> _dims)
  {
     dims.SetLength(_dims.size());
     long i;
     std::initializer_list<long>::iterator it;
     for (i=0, it=_dims.begin(); it!=_dims.end(); ++i, ++it)
       dims[i] = *it;
 
     [...] // continue as above, initialize prods
  }
  */
 
  long getNumDims() const { return dims.length(); }
 
  long getSize() const { return ((getNumDims() > 0) ? prods[0] : 1); }
 
  long getDim(long d) const { return dims.at(d); }
 
  long getProd(long d) const { return prods.at(d); }
 
  long getProd(long from, long to) const
  {
    return prods.at(from) / prods.at(to);
  }
 
  long getCoord(long i, long d) const
  {
    assertInRange(i,
                  0l,
                  getSize(),
                  "Coordinate does not exist (index i out of range)");
 
    return (i % prods.at(d)) / prods.at(d + 1);
  }
 
  long addCoord(long i, long d, long offset) const
  {
    assertInRange(i,
                  0l,
                  getSize(),
                  "Coordinate does not exist (index i out of range)");
 
    offset = offset % dims.at(d);
    if (offset < 0)
      offset += dims.at(d);
 
    long i_d = getCoord(i, d);
    long i_d1 = (i_d + offset) % dims.at(d);
 
    long i1 = i + (i_d1 - i_d) * prods.at(d + 1);
 
    return i1;
  }
 
  template <typename VecType>
  bool incrementCoords(VecType& v) const
  {
    for (long i = getNumDims() - 1; i >= 0; --i) {
      if (i >= lsize(v))
        continue; // sanity check
 
      // increment current index, set all the ones after it to zero
      if (long(v[i]) < getDim(i) - 1) {
        v[i]++;
        for (long j = i + 1; j < lsize(v); j++)
          v[j] = 0;
        return true; // succeeded in incrementing the vector
      }
      // if buffer[i] >= getDim(i)-1, move to previous index i
    }
    return false; // cannot increment the vector anymore
  }
 
  template <typename VecType>
  void getAllCoords(VecType& v, long i) const
  {
    assertInRange(i,
                  0l,
                  getSize(),
                  "Coordinate does not exist (index i out of range)");
 
    resize(v, getNumDims()); // resize(*), lsize(*) defined in NumbTh.h
    for (long j = getNumDims() - 1; j >= 0; --j) {
      v[j] = i % getDim(j);
      i = (i - v[j]) / getDim(j);
    }
  }
 
  template <typename VecType>
  long assembleCoords(VecType& v) const
  {
    assertEq(lsize(v),
             getNumDims(),
             "Vector size is different to the number of dimensions");
 
    long idx = 0;
    for (long i = 0; i < getNumDims(); i++) {
      idx += v[i] * prods[i + 1];
    }
    return idx;
  }
 
  long numSlices(long d = 1) const { return getProd(0, d); }
 
  long sliceSize(long d = 1) const { return getProd(d); }
 
  long numCols() const { return getProd(1); }
 
  std::pair<long, long> breakIndexByDim(long idx, long dim) const;
 
  long assembleIndexByDim(std::pair<long, long> idx, long dim) const;
 
  friend std::ostream& operator<<(std::ostream& s, const CubeSignature& sig);
};
 
inline std::ostream& operator<<(std::ostream& s, const CubeSignature& sig)
{
  return s << sig.dims;
}
 
template <typename T>
class HyperCube
{
private:
  const CubeSignature& sig;
  NTL::Vec<T> data;
 
  HyperCube(); // disable default constructor
 
public:
  // Default copy-constructor
  HyperCube(const HyperCube& other) = default;
 
  HyperCube(const CubeSignature& _sig) : sig(_sig)
  {
    data.FixLength(sig.getSize());
  }
 
  // use default copy constructor
 
  HyperCube& operator=(const HyperCube<T>& other)
  {
    assertEq(&this->sig,
             &other.sig,
             "Cannot assign HyperCubes with different signatures");
    data = other.data;
    return *this;
  }
 
  bool operator==(const HyperCube<T>& other) const
  {
    assertEq(&this->sig,
             &other.sig,
             "Cannot compare HyperCubes with different signatures");
    return data == other.data;
  }
 
  bool operator!=(const HyperCube<T>& other) const { return !(*this == other); }
 
  const CubeSignature& getSig() const { return sig; }
 
  NTL::Vec<T>& getData() { return data; }
 
  const NTL::Vec<T>& getData() const { return data; }
 
  long getSize() const { return sig.getSize(); }
 
  long getNumDims() const { return sig.getNumDims(); }
 
  long getDim(long d) const { return sig.getDim(d); }
 
  long getProd(long d) const { return sig.getProd(d); }
 
  long getProd(long from, long to) const { return sig.getProd(from, to); }
 
  long getCoord(long i, long d) const { return sig.getCoord(i, d); }
 
  long addCoord(long i, long d, long offset) const
  {
    return sig.addCoord(i, d, offset);
  }
 
  long numSlices(long d = 1) const { return getProd(0, d); }
 
  long sliceSize(long d = 1) const { return getProd(d); }
 
  long numCols() const { return getProd(1); }
 
  T& at(long i) { return data.at(i); }
 
  T& operator[](long i) { return data[i]; }
 
  const T& at(long i) const { return data.at(i); }
 
  const T& operator[](long i) const { return data[i]; }
 
  void rotate1D(long i, long k);
 
  void shift1D(long i, long k);
};
 
 
template <typename T>
class ConstCubeSlice
{
private:
  const NTL::Vec<T>* data;
  const CubeSignature* sig;
  long dimOffset;  // # of "missing dimensions" is this slice vs. the full cube
  long sizeOffset; // index in the cube of the first element in this slice
 
  ConstCubeSlice(); // disabled
 
public:
  explicit ConstCubeSlice(const HyperCube<T>& _cube)
  {
    data = &_cube.getData();
    sig = &_cube.getSig();
    dimOffset = 0;
    sizeOffset = 0;
  }
 
  ConstCubeSlice(const NTL::Vec<T>& _data, const CubeSignature& _sig)
  {
    assertEq<InvalidArgument>(_data.length(),
                              _sig.getSize(),
                              "Data and signature sizes are different");
    data = &_data;
    sig = &_sig;
    dimOffset = 0;
    sizeOffset = 0;
  }
 
  ConstCubeSlice(const ConstCubeSlice& bigger, long i, long _dimOffset = 1);
  ConstCubeSlice(const HyperCube<T>& _cube, long i, long _dimOffset = 1);
 
  // use default copy constructor and assignment operators,
  // which means shallow copy
 
  // the following mimic the corresponding methods
  // in the HyperCube class, restricted to the slice
 
  const CubeSignature& getSig() const { return *sig; }
 
  long getSize() const { return sig->getProd(dimOffset); }
 
  long getNumDims() const { return sig->getNumDims() - dimOffset; }
 
  long getDim(long d) const { return sig->getDim(d + dimOffset); }
 
  long getProd(long d) const { return sig->getProd(d + dimOffset); }
 
  long getProd(long from, long to) const
  {
    return sig->getProd(from + dimOffset, to + dimOffset);
  }
 
  long getCoord(long i, long d) const
  {
    assertInRange(i,
                  0l,
                  getSize(),
                  "Coordinate does not exist (index i out of range)");
    return sig->getCoord(i + sizeOffset, d + dimOffset);
  }
 
  long addCoord(long i, long d, long offset) const
  {
    assertInRange(i,
                  0l,
                  getSize(),
                  "Coordinate does not exist (index i out of range)");
    return sig->addCoord(i + sizeOffset, d + dimOffset, offset);
  }
 
  long numSlices(long d = 1) const { return getProd(0, d); }
 
  long sliceSize(long d = 1) const { return getProd(d); }
 
  long numCols() const { return getProd(1); }
 
  const T& at(long i) const
  {
    assertInRange(i,
                  0l,
                  getSize(),
                  "Coordinate does not exist (index i out of range)");
    return (*data)[i + sizeOffset];
  }
 
  const T& operator[](long i) const { return (*data)[i + sizeOffset]; }
};
 
template <typename T>
class CubeSlice : public ConstCubeSlice<T>
{
private:
  CubeSlice(); // disabled
public:
  // initialize the slice to the full cube
  explicit CubeSlice(HyperCube<T>& _cube) : ConstCubeSlice<T>(_cube) {}
 
  CubeSlice(NTL::Vec<T>& _data, const CubeSignature& _sig) :
      ConstCubeSlice<T>(_data, _sig)
  {}
 
  // initialize the slice to point to the i-th subcube
  // of the cube pointed to by bigger
  CubeSlice(const CubeSlice<T>& bigger, long i, long _dimOffset = 1) :
      ConstCubeSlice<T>(bigger, i, _dimOffset)
  {}
 
  CubeSlice(HyperCube<T>& _cube, long i, long _dimOffset = 1) :
      ConstCubeSlice<T>(_cube, i, _dimOffset)
  {}
 
  // deep copy of a slice: copies other into this
  void copy(const ConstCubeSlice<T>& other) const;
 
  // reference to element at position i, with bounds check
  T& at(long i) const { return const_cast<T&>(this->ConstCubeSlice<T>::at(i)); }
 
  // reference to element at position i, without bounds check
  T& operator[](long i) const
  {
    return const_cast<T&>(this->ConstCubeSlice<T>::operator[](i));
  }
};
 
template <typename T>
void getHyperColumn(NTL::Vec<T>& v, const ConstCubeSlice<T>& s, long pos);
 
template <typename T>
void setHyperColumn(const NTL::Vec<T>& v, const CubeSlice<T>& s, long pos);
 
template <typename T>
void setHyperColumn(const NTL::Vec<T>& v,
                    const CubeSlice<T>& s,
                    long pos,
                    const T& val);
 
template <typename T>
void print3D(const HyperCube<T>& c);
 
} // namespace helib
 
#endif // ifndef HELIB_HYPERCUBE_H