matpackV.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */

#include "matpackV.h"
#include "exceptions.h"

using std::runtime_error;

// Functions for ConstTensor5View:
// ------------------------------


bool ConstTensor5View::empty() const {
  return (nshelves() == 0 || nbooks() == 0 || npages() == 0 || nrows() == 0 ||
          ncols() == 0);
}

Index ConstTensor5View::nshelves() const { return msr.mextent; }

Index ConstTensor5View::nbooks() const { return mbr.mextent; }

Index ConstTensor5View::npages() const { return mpr.mextent; }

Index ConstTensor5View::nrows() const { return mrr.mextent; }

Index ConstTensor5View::ncols() const { return mcr.mextent; }

ConstTensor5View ConstTensor5View::operator()(const Range& s,
                                              const Range& b,
                                              const Range& p,
                                              const Range& r,
                                              const Range& c) const {
  return ConstTensor5View(mdata, msr, mbr, mpr, mrr, mcr, s, b, p, r, c);
}

ConstTensor4View ConstTensor5View::operator()(const Range& s,
                                              const Range& b,
                                              const Range& p,
                                              const Range& r,
                                              Index c) const {
  // Check that c is valid:
  assert(0 <= c);
  assert(c < mcr.mextent);

  return ConstTensor4View(
      mdata + mcr.mstart + c * mcr.mstride, msr, mbr, mpr, mrr, s, b, p, r);
}

ConstTensor4View ConstTensor5View::operator()(const Range& s,
                                              const Range& b,
                                              const Range& p,
                                              Index r,
                                              const Range& c) const {
  // Check that r is valid:
  assert(0 <= r);
  assert(r < mrr.mextent);

  return ConstTensor4View(
      mdata + mrr.mstart + r * mrr.mstride, msr, mbr, mpr, mcr, s, b, p, c);
}

ConstTensor4View ConstTensor5View::operator()(const Range& s,
                                              const Range& b,
                                              Index p,
                                              const Range& r,
                                              const Range& c) const {
  // Check that p is valid:
  assert(0 <= p);
  assert(p < mpr.mextent);

  return ConstTensor4View(
      mdata + mpr.mstart + p * mpr.mstride, msr, mbr, mrr, mcr, s, b, r, c);
}

ConstTensor4View ConstTensor5View::operator()(const Range& s,
                                              Index b,
                                              const Range& p,
                                              const Range& r,
                                              const Range& c) const {
  // Check that b is valid:
  assert(0 <= b);
  assert(b < mbr.mextent);

  return ConstTensor4View(
      mdata + mbr.mstart + b * mbr.mstride, msr, mpr, mrr, mcr, s, p, r, c);
}

ConstTensor4View ConstTensor5View::operator()(Index s,
                                              const Range& b,
                                              const Range& p,
                                              const Range& r,
                                              const Range& c) const {
  // Check that s is valid:
  assert(0 <= s);
  assert(s < msr.mextent);

  return ConstTensor4View(
      mdata + msr.mstart + s * msr.mstride, mbr, mpr, mrr, mcr, b, p, r, c);
}

ConstTensor3View ConstTensor5View::operator()(
    const Range& s, const Range& b, const Range& p, Index r, Index c) const {
  // Check that r and c is valid:
  assert(0 <= r);
  assert(0 <= c);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return ConstTensor3View(
      mdata + mrr.mstart + r * mrr.mstride + mcr.mstart + c * mcr.mstride,
      msr,
      mbr,
      mpr,
      s,
      b,
      p);
}

ConstTensor3View ConstTensor5View::operator()(
    const Range& s, const Range& b, Index p, const Range& r, Index c) const {
  // Check that p and c are valid:
  assert(0 <= p);
  assert(0 <= c);
  assert(p < mpr.mextent);
  assert(c < mcr.mextent);

  return ConstTensor3View(
      mdata + mpr.mstart + p * mpr.mstride + mcr.mstart + c * mcr.mstride,
      msr,
      mbr,
      mrr,
      s,
      b,
      r);
}

ConstTensor3View ConstTensor5View::operator()(
    const Range& s, const Range& b, Index p, Index r, const Range& c) const {
  // Check that p and r are valid:
  assert(0 <= p);
  assert(0 <= r);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);

  return ConstTensor3View(
      mdata + mpr.mstart + p * mpr.mstride + mrr.mstart + r * mrr.mstride,
      msr,
      mbr,
      mcr,
      s,
      b,
      c);
}

ConstTensor3View ConstTensor5View::operator()(
    const Range& s, Index b, const Range& p, Index r, const Range& c) const {
  // Check that b and r are valid:
  assert(0 <= b);
  assert(0 <= r);
  assert(b < mbr.mextent);
  assert(r < mrr.mextent);

  return ConstTensor3View(
      mdata + mbr.mstart + b * mbr.mstride + mrr.mstart + r * mrr.mstride,
      msr,
      mpr,
      mcr,
      s,
      p,
      c);
}

ConstTensor3View ConstTensor5View::operator()(
    const Range& s, Index b, const Range& p, const Range& r, Index c) const {
  // Check that b and c are valid:
  assert(0 <= b);
  assert(0 <= c);
  assert(b < mbr.mextent);
  assert(c < mcr.mextent);

  return ConstTensor3View(
      mdata + mbr.mstart + b * mbr.mstride + mcr.mstart + c * mcr.mstride,
      msr,
      mpr,
      mrr,
      s,
      p,
      r);
}

ConstTensor3View ConstTensor5View::operator()(
    const Range& s, Index b, Index p, const Range& r, const Range& c) const {
  // Check that b and p are valid:
  assert(0 <= b);
  assert(0 <= p);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);

  return ConstTensor3View(
      mdata + mbr.mstart + b * mbr.mstride + mpr.mstart + p * mpr.mstride,
      msr,
      mrr,
      mcr,
      s,
      r,
      c);
}

ConstTensor3View ConstTensor5View::operator()(
    Index s, const Range& b, Index p, const Range& r, const Range& c) const {
  // Check that s and p are valid:
  assert(0 <= s);
  assert(0 <= p);
  assert(s < msr.mextent);
  assert(p < mpr.mextent);

  return ConstTensor3View(
      mdata + msr.mstart + s * msr.mstride + mpr.mstart + p * mpr.mstride,
      mbr,
      mrr,
      mcr,
      b,
      r,
      c);
}

ConstTensor3View ConstTensor5View::operator()(
    Index s, const Range& b, const Range& p, Index r, const Range& c) const {
  // Check that s and r are valid:
  assert(0 <= s);
  assert(0 <= r);
  assert(s < msr.mextent);
  assert(r < mrr.mextent);

  return ConstTensor3View(
      mdata + msr.mstart + s * msr.mstride + mrr.mstart + r * mrr.mstride,
      mbr,
      mpr,
      mcr,
      b,
      p,
      c);
}

ConstTensor3View ConstTensor5View::operator()(
    Index s, const Range& b, const Range& p, const Range& r, Index c) const {
  // Check that s and c are valid:
  assert(0 <= s);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(c < mcr.mextent);

  return ConstTensor3View(
      mdata + msr.mstart + s * msr.mstride + mcr.mstart + c * mcr.mstride,
      mbr,
      mpr,
      mrr,
      b,
      p,
      r);
}

ConstTensor3View ConstTensor5View::operator()(
    Index s, Index b, const Range& p, const Range& r, const Range& c) const {
  // Check that s and b are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);

  return ConstTensor3View(
      mdata + msr.mstart + s * msr.mstride + mbr.mstart + b * mbr.mstride,
      mpr,
      mrr,
      mcr,
      p,
      r,
      c);
}

ConstMatrixView ConstTensor5View::operator()(
    const Range& s, const Range& b, Index p, Index r, Index c) const {
  // Check that p, r and c are valid:
  assert(0 <= p);
  assert(0 <= r);
  assert(0 <= c);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return ConstMatrixView(mdata + mpr.mstart + p * mpr.mstride + mrr.mstart +
                             r * mrr.mstride + mcr.mstart + c * mcr.mstride,
                         msr,
                         mbr,
                         s,
                         b);
}

ConstMatrixView ConstTensor5View::operator()(
    const Range& s, Index b, const Range& p, Index r, Index c) const {
  // Check that b, r and c are valid:
  assert(0 <= b);
  assert(0 <= r);
  assert(0 <= c);
  assert(b < mbr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return ConstMatrixView(mdata + mbr.mstart + b * mbr.mstride + mrr.mstart +
                             r * mpr.mstride + mcr.mstart + c * mcr.mstride,
                         msr,
                         mpr,
                         s,
                         p);
}

ConstMatrixView ConstTensor5View::operator()(
    const Range& s, Index b, Index p, const Range& r, Index c) const {
  // Check that b, p and c are valid:
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= c);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(c < mcr.mextent);

  return ConstMatrixView(mdata + mbr.mstart + b * mbr.mstride + mpr.mstart +
                             p * mpr.mstride + mcr.mstart + c * mcr.mstride,
                         msr,
                         mrr,
                         s,
                         r);
}

ConstMatrixView ConstTensor5View::operator()(
    const Range& s, Index b, Index p, Index r, const Range& c) const {
  // Check that b, p and r are valid:
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= r);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);

  return ConstMatrixView(mdata + mbr.mstart + b * mbr.mstride + mpr.mstart +
                             p * mpr.mstride + mrr.mstart + r * mrr.mstride,
                         msr,
                         mcr,
                         s,
                         c);
}

ConstMatrixView ConstTensor5View::operator()(
    Index s, const Range& b, Index p, Index r, const Range& c) const {
  // Check that s, p and r are valid:
  assert(0 <= s);
  assert(0 <= p);
  assert(0 <= r);
  assert(s < msr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);

  return ConstMatrixView(mdata + msr.mstart + s * msr.mstride + mpr.mstart +
                             p * mpr.mstride + mrr.mstart + r * mrr.mstride,
                         mbr,
                         mcr,
                         b,
                         c);
}

ConstMatrixView ConstTensor5View::operator()(
    Index s, const Range& b, Index p, const Range& r, Index c) const {
  // Check that s, p and c are valid:
  assert(0 <= s);
  assert(0 <= p);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(p < mpr.mextent);
  assert(c < mcr.mextent);

  return ConstMatrixView(mdata + msr.mstart + s * msr.mstride + mpr.mstart +
                             p * mpr.mstride + mcr.mstart + c * mcr.mstride,
                         mbr,
                         mrr,
                         b,
                         r);
}

ConstMatrixView ConstTensor5View::operator()(
    Index s, const Range& b, const Range& p, Index r, Index c) const {
  // Check that s, r and c are valid:
  assert(0 <= s);
  assert(0 <= r);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return ConstMatrixView(mdata + msr.mstart + s * msr.mstride + mrr.mstart +
                             r * mrr.mstride + mcr.mstart + c * mcr.mstride,
                         mbr,
                         mpr,
                         b,
                         p);
}

ConstMatrixView ConstTensor5View::operator()(
    Index s, Index b, const Range& p, const Range& r, Index c) const {
  // Check that s, b and c are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(c < mcr.mextent);

  return ConstMatrixView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                             b * mbr.mstride + mcr.mstart + c * mcr.mstride,
                         mpr,
                         mrr,
                         p,
                         r);
}

ConstMatrixView ConstTensor5View::operator()(
    Index s, Index b, const Range& p, Index r, const Range& c) const {
  // Check that s, b and r are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= r);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(r < mrr.mextent);

  return ConstMatrixView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                             b * mbr.mstride + mrr.mstart + r * mrr.mstride,
                         mpr,
                         mcr,
                         p,
                         c);
}

ConstMatrixView ConstTensor5View::operator()(
    Index s, Index b, Index p, const Range& r, const Range& c) const {
  // Check that s, b and p are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= p);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);

  return ConstMatrixView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                             b * mbr.mstride + mpr.mstart + p * mpr.mstride,
                         mrr,
                         mcr,
                         r,
                         c);
}

ConstVectorView ConstTensor5View::operator()(
    const Range& s, Index b, Index p, Index r, Index c) const {
  // Check that b, p, r and c are valid:
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= r);
  assert(0 <= c);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return ConstVectorView(mdata + mbr.mstart + b * mbr.mstride + mpr.mstart +
                             p * mpr.mstride + mrr.mstart + r * mrr.mstride +
                             mcr.mstart + c * mcr.mstride,
                         msr,
                         s);
}

ConstVectorView ConstTensor5View::operator()(
    Index s, const Range& b, Index p, Index r, Index c) const {
  // Check that s, p, r and c are valid:
  assert(0 <= s);
  assert(0 <= p);
  assert(0 <= r);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return ConstVectorView(mdata + msr.mstart + s * msr.mstride + mpr.mstart +
                             p * mpr.mstride + mrr.mstart + r * mrr.mstride +
                             mcr.mstart + c * mcr.mstride,
                         mbr,
                         b);
}

ConstVectorView ConstTensor5View::operator()(
    Index s, Index b, const Range& p, Index r, Index c) const {
  // Check that s, b, r and c are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= r);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return ConstVectorView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                             b * mbr.mstride + mrr.mstart + r * mrr.mstride +
                             mcr.mstart + c * mcr.mstride,
                         mpr,
                         p);
}

ConstVectorView ConstTensor5View::operator()(
    Index s, Index b, Index p, const Range& r, Index c) const {
  // Check that s, b, p and c are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(c < mcr.mextent);

  return ConstVectorView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                             b * mbr.mstride + mpr.mstart + p * mpr.mstride +
                             mcr.mstart + c * mcr.mstride,
                         mrr,
                         r);
}

ConstVectorView ConstTensor5View::operator()(
    Index s, Index b, Index p, Index r, const Range& c) const {
  // Check that s, b, p and r are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= r);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);

  return ConstVectorView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                             b * mbr.mstride + mpr.mstart + p * mpr.mstride +
                             mrr.mstart + r * mrr.mstride,
                         mcr,
                         c);
}

Numeric* Tensor5View::get_c_array() {
  if (msr.mstart != 0 ||
      msr.mstride != mbr.mextent * mpr.mextent * mrr.mextent * mcr.mextent ||
      mbr.mstart != 0 ||
      mbr.mstride != mpr.mextent * mrr.mextent * mcr.mextent ||
      mpr.mstart != 0 || mpr.mstride != mrr.mextent * mcr.mextent ||
      mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
      mcr.mstride != 1)
    throw std::runtime_error(
        "A Tensor5View can only be converted to a plain C-array if it's pointing to a continuous block of data");

  return mdata;
}

const Numeric* Tensor5View::get_c_array() const {
  if (msr.mstart != 0 ||
      msr.mstride != mbr.mextent * mpr.mextent * mrr.mextent * mcr.mextent ||
      mbr.mstart != 0 ||
      mbr.mstride != mpr.mextent * mrr.mextent * mcr.mextent ||
      mpr.mstart != 0 || mpr.mstride != mrr.mextent * mcr.mextent ||
      mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
      mcr.mstride != 1)
    throw std::runtime_error(
        "A Tensor5View can only be converted to a plain C-array if it's pointing to a continuous block of data");

  return mdata;
}

ConstIterator5D ConstTensor5View::begin() const {
  return ConstIterator5D(
      ConstTensor4View(mdata + msr.mstart, mbr, mpr, mrr, mcr), msr.mstride);
}

ConstIterator5D ConstTensor5View::end() const {
  return ConstIterator5D(
      ConstTensor4View(
          mdata + msr.mstart + (msr.mextent) * msr.mstride, mbr, mpr, mrr, mcr),
      msr.mstride);
}

ConstTensor5View::ConstTensor5View(const ConstTensor4View& a)
    : msr(0, 1, a.mbr.mextent * a.mpr.mextent * a.mrr.mextent * a.mcr.mextent),
      mbr(a.mbr),
      mpr(a.mpr),
      mrr(a.mrr),
      mcr(a.mcr),
      mdata(a.mdata) {
  // Nothing to do here.
}

ConstTensor5View::ConstTensor5View(Numeric* data,
                                   const Range& sr,
                                   const Range& br,
                                   const Range& pr,
                                   const Range& rr,
                                   const Range& cr)
    : msr(sr), mbr(br), mpr(pr), mrr(rr), mcr(cr), mdata(data) {
  // Nothing to do here.
}

ConstTensor5View::ConstTensor5View(Numeric* data,
                                   const Range& ps,
                                   const Range& pb,
                                   const Range& pp,
                                   const Range& pr,
                                   const Range& pc,
                                   const Range& ns,
                                   const Range& nb,
                                   const Range& np,
                                   const Range& nr,
                                   const Range& nc)
    : msr(ps, ns),
      mbr(pb, nb),
      mpr(pp, np),
      mrr(pr, nr),
      mcr(pc, nc),
      mdata(data) {
  // Nothing to do here.
}

std::ostream& operator<<(std::ostream& os, const ConstTensor5View& v) {
  // Page iterators:
  ConstIterator5D is = v.begin();
  const ConstIterator5D end_shelf = v.end();

  if (is != end_shelf) {
    os << *is;
    ++is;
  }

  for (; is != end_shelf; ++is) {
    os << "\n\n";
    os << *is;
  }

  return os;
}

// Functions for Tensor5View:
// -------------------------

Tensor5View Tensor5View::operator()(const Range& s,
                                    const Range& b,
                                    const Range& p,
                                    const Range& r,
                                    const Range& c) {
  return Tensor5View(mdata, msr, mbr, mpr, mrr, mcr, s, b, p, r, c);
}

Tensor4View Tensor5View::operator()(
    const Range& s, const Range& b, const Range& p, const Range& r, Index c) {
  // Check that c is valid:
  assert(0 <= c);
  assert(c < mcr.mextent);

  return Tensor4View(
      mdata + mcr.mstart + c * mcr.mstride, msr, mbr, mpr, mrr, s, b, p, r);
}

Tensor4View Tensor5View::operator()(
    const Range& s, const Range& b, const Range& p, Index r, const Range& c) {
  // Check that r is valid:
  assert(0 <= r);
  assert(r < mrr.mextent);

  return Tensor4View(
      mdata + mrr.mstart + r * mrr.mstride, msr, mbr, mpr, mcr, s, b, p, c);
}

Tensor4View Tensor5View::operator()(
    const Range& s, const Range& b, Index p, const Range& r, const Range& c) {
  // Check that p is valid:
  assert(0 <= p);
  assert(p < mpr.mextent);

  return Tensor4View(
      mdata + mpr.mstart + p * mpr.mstride, msr, mbr, mrr, mcr, s, b, r, c);
}

Tensor4View Tensor5View::operator()(
    const Range& s, Index b, const Range& p, const Range& r, const Range& c) {
  // Check that b is valid:
  assert(0 <= b);
  assert(b < mbr.mextent);

  return Tensor4View(
      mdata + mbr.mstart + b * mbr.mstride, msr, mpr, mrr, mcr, s, p, r, c);
}

Tensor4View Tensor5View::operator()(
    Index s, const Range& b, const Range& p, const Range& r, const Range& c) {
  // Check that s is valid:
  assert(0 <= s);
  assert(s < msr.mextent);

  return Tensor4View(
      mdata + msr.mstart + s * msr.mstride, mbr, mpr, mrr, mcr, b, p, r, c);
}

Tensor3View Tensor5View::operator()(
    const Range& s, const Range& b, const Range& p, Index r, Index c) {
  // Check that r and c is valid:
  assert(0 <= r);
  assert(0 <= c);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return Tensor3View(
      mdata + mrr.mstart + r * mrr.mstride + mcr.mstart + c * mcr.mstride,
      msr,
      mbr,
      mpr,
      s,
      b,
      p);
}

Tensor3View Tensor5View::operator()(
    const Range& s, const Range& b, Index p, const Range& r, Index c) {
  // Check that p and c are valid:
  assert(0 <= p);
  assert(0 <= c);
  assert(p < mpr.mextent);
  assert(c < mcr.mextent);

  return Tensor3View(
      mdata + mpr.mstart + p * mpr.mstride + mcr.mstart + c * mcr.mstride,
      msr,
      mbr,
      mrr,
      s,
      b,
      r);
}

Tensor3View Tensor5View::operator()(
    const Range& s, const Range& b, Index p, Index r, const Range& c) {
  // Check that p and r are valid:
  assert(0 <= p);
  assert(0 <= r);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);

  return Tensor3View(
      mdata + mpr.mstart + p * mpr.mstride + mrr.mstart + r * mrr.mstride,
      msr,
      mbr,
      mcr,
      s,
      b,
      c);
}

Tensor3View Tensor5View::operator()(
    const Range& s, Index b, const Range& p, Index r, const Range& c) {
  // Check that b and r are valid:
  assert(0 <= b);
  assert(0 <= r);
  assert(b < mbr.mextent);
  assert(r < mrr.mextent);

  return Tensor3View(
      mdata + mbr.mstart + b * mbr.mstride + mrr.mstart + r * mrr.mstride,
      msr,
      mpr,
      mcr,
      s,
      p,
      c);
}

Tensor3View Tensor5View::operator()(
    const Range& s, Index b, const Range& p, const Range& r, Index c) {
  // Check that b and c are valid:
  assert(0 <= b);
  assert(0 <= c);
  assert(b < mbr.mextent);
  assert(c < mcr.mextent);

  return Tensor3View(
      mdata + mbr.mstart + b * mbr.mstride + mcr.mstart + c * mcr.mstride,
      msr,
      mpr,
      mrr,
      s,
      p,
      r);
}

Tensor3View Tensor5View::operator()(
    const Range& s, Index b, Index p, const Range& r, const Range& c) {
  // Check that b and p are valid:
  assert(0 <= b);
  assert(0 <= p);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);

  return Tensor3View(
      mdata + mbr.mstart + b * mbr.mstride + mpr.mstart + p * mpr.mstride,
      msr,
      mrr,
      mcr,
      s,
      r,
      c);
}

Tensor3View Tensor5View::operator()(
    Index s, const Range& b, Index p, const Range& r, const Range& c) {
  // Check that s and p are valid:
  assert(0 <= s);
  assert(0 <= p);
  assert(s < msr.mextent);
  assert(p < mpr.mextent);

  return Tensor3View(
      mdata + msr.mstart + s * msr.mstride + mpr.mstart + p * mpr.mstride,
      mbr,
      mrr,
      mcr,
      b,
      r,
      c);
}

Tensor3View Tensor5View::operator()(
    Index s, const Range& b, const Range& p, Index r, const Range& c) {
  // Check that s and r are valid:
  assert(0 <= s);
  assert(0 <= r);
  assert(s < msr.mextent);
  assert(r < mrr.mextent);

  return Tensor3View(
      mdata + msr.mstart + s * msr.mstride + mrr.mstart + r * mrr.mstride,
      mbr,
      mpr,
      mcr,
      b,
      p,
      c);
}

Tensor3View Tensor5View::operator()(
    Index s, const Range& b, const Range& p, const Range& r, Index c) {
  // Check that s and c are valid:
  assert(0 <= s);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(c < mcr.mextent);

  return Tensor3View(
      mdata + msr.mstart + s * msr.mstride + mcr.mstart + c * mcr.mstride,
      mbr,
      mpr,
      mrr,
      b,
      p,
      r);
}

Tensor3View Tensor5View::operator()(
    Index s, Index b, const Range& p, const Range& r, const Range& c) {
  // Check that s and b are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);

  return Tensor3View(
      mdata + msr.mstart + s * msr.mstride + mbr.mstart + b * mbr.mstride,
      mpr,
      mrr,
      mcr,
      p,
      r,
      c);
}

MatrixView Tensor5View::operator()(
    const Range& s, const Range& b, Index p, Index r, Index c) {
  // Check that p, r and c are valid:
  assert(0 <= p);
  assert(0 <= r);
  assert(0 <= c);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return MatrixView(mdata + mpr.mstart + p * mpr.mstride + mrr.mstart +
                        r * mrr.mstride + mcr.mstart + c * mcr.mstride,
                    msr,
                    mbr,
                    s,
                    b);
}

MatrixView Tensor5View::operator()(
    const Range& s, Index b, const Range& p, Index r, Index c) {
  // Check that b, r and c are valid:
  assert(0 <= b);
  assert(0 <= r);
  assert(0 <= c);
  assert(b < mbr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return MatrixView(mdata + mbr.mstart + b * mbr.mstride + mrr.mstart +
                        r * mpr.mstride + mcr.mstart + c * mcr.mstride,
                    msr,
                    mpr,
                    s,
                    p);
}

MatrixView Tensor5View::operator()(
    const Range& s, Index b, Index p, const Range& r, Index c) {
  // Check that b, p and c are valid:
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= c);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(c < mcr.mextent);

  return MatrixView(mdata + mbr.mstart + b * mbr.mstride + mpr.mstart +
                        p * mpr.mstride + mcr.mstart + c * mcr.mstride,
                    msr,
                    mrr,
                    s,
                    r);
}

MatrixView Tensor5View::operator()(
    const Range& s, Index b, Index p, Index r, const Range& c) {
  // Check that b, p and r are valid:
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= r);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);

  return MatrixView(mdata + mbr.mstart + b * mbr.mstride + mpr.mstart +
                        p * mpr.mstride + mrr.mstart + r * mrr.mstride,
                    msr,
                    mcr,
                    s,
                    c);
}

MatrixView Tensor5View::operator()(
    Index s, const Range& b, Index p, Index r, const Range& c) {
  // Check that s, p and r are valid:
  assert(0 <= s);
  assert(0 <= p);
  assert(0 <= r);
  assert(s < msr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);

  return MatrixView(mdata + msr.mstart + s * msr.mstride + mpr.mstart +
                        p * mpr.mstride + mrr.mstart + r * mrr.mstride,
                    mbr,
                    mcr,
                    b,
                    c);
}

MatrixView Tensor5View::operator()(
    Index s, const Range& b, Index p, const Range& r, Index c) {
  // Check that s, p and c are valid:
  assert(0 <= s);
  assert(0 <= p);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(p < mpr.mextent);
  assert(c < mcr.mextent);

  return MatrixView(mdata + msr.mstart + s * msr.mstride + mpr.mstart +
                        p * mpr.mstride + mcr.mstart + c * mcr.mstride,
                    mbr,
                    mrr,
                    b,
                    r);
}

MatrixView Tensor5View::operator()(
    Index s, const Range& b, const Range& p, Index r, Index c) {
  // Check that s, r and c are valid:
  assert(0 <= s);
  assert(0 <= r);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return MatrixView(mdata + msr.mstart + s * msr.mstride + mrr.mstart +
                        r * mrr.mstride + mcr.mstart + c * mcr.mstride,
                    mbr,
                    mpr,
                    b,
                    p);
}

MatrixView Tensor5View::operator()(
    Index s, Index b, const Range& p, const Range& r, Index c) {
  // Check that s, b and c are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(c < mcr.mextent);

  return MatrixView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                        b * mbr.mstride + mcr.mstart + c * mcr.mstride,
                    mpr,
                    mrr,
                    p,
                    r);
}

MatrixView Tensor5View::operator()(
    Index s, Index b, const Range& p, Index r, const Range& c) {
  // Check that s, b and r are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= r);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(r < mrr.mextent);

  return MatrixView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                        b * mbr.mstride + mrr.mstart + r * mrr.mstride,
                    mpr,
                    mcr,
                    p,
                    c);
}

MatrixView Tensor5View::operator()(
    Index s, Index b, Index p, const Range& r, const Range& c) {
  // Check that s, b and p are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= p);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);

  return MatrixView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                        b * mbr.mstride + mpr.mstart + p * mpr.mstride,
                    mrr,
                    mcr,
                    r,
                    c);
}

VectorView Tensor5View::operator()(
    const Range& s, Index b, Index p, Index r, Index c) {
  // Check that b, p, r and c are valid:
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= r);
  assert(0 <= c);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return VectorView(mdata + mbr.mstart + b * mbr.mstride + mpr.mstart +
                        p * mpr.mstride + mrr.mstart + r * mrr.mstride +
                        mcr.mstart + c * mcr.mstride,
                    msr,
                    s);
}

VectorView Tensor5View::operator()(
    Index s, const Range& b, Index p, Index r, Index c) {
  // Check that s, p, r and c are valid:
  assert(0 <= s);
  assert(0 <= p);
  assert(0 <= r);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return VectorView(mdata + msr.mstart + s * msr.mstride + mpr.mstart +
                        p * mpr.mstride + mrr.mstart + r * mrr.mstride +
                        mcr.mstart + c * mcr.mstride,
                    mbr,
                    b);
}

VectorView Tensor5View::operator()(
    Index s, Index b, const Range& p, Index r, Index c) {
  // Check that s, b, r and c are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= r);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(r < mrr.mextent);
  assert(c < mcr.mextent);

  return VectorView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                        b * mbr.mstride + mrr.mstart + r * mrr.mstride +
                        mcr.mstart + c * mcr.mstride,
                    mpr,
                    p);
}

VectorView Tensor5View::operator()(
    Index s, Index b, Index p, const Range& r, Index c) {
  // Check that s, b, p and c are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= c);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(c < mcr.mextent);

  return VectorView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                        b * mbr.mstride + mpr.mstart + p * mpr.mstride +
                        mcr.mstart + c * mcr.mstride,
                    mrr,
                    r);
}

VectorView Tensor5View::operator()(
    Index s, Index b, Index p, Index r, const Range& c) {
  // Check that s, b, p and r are valid:
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= r);
  assert(s < msr.mextent);
  assert(b < mbr.mextent);
  assert(p < mpr.mextent);
  assert(r < mrr.mextent);

  return VectorView(mdata + msr.mstart + s * msr.mstride + mbr.mstart +
                        b * mbr.mstride + mpr.mstart + p * mpr.mstride +
                        mrr.mstart + r * mrr.mstride,
                    mcr,
                    c);
}

Iterator5D Tensor5View::begin() {
  return Iterator5D(Tensor4View(mdata + msr.mstart, mbr, mpr, mrr, mcr),
                    msr.mstride);
}

Iterator5D Tensor5View::end() {
  return Iterator5D(
      Tensor4View(
          mdata + msr.mstart + (msr.mextent) * msr.mstride, mbr, mpr, mrr, mcr),
      msr.mstride);
}

Tensor5View& Tensor5View::operator=(const ConstTensor5View& m) {
  // Check that sizes are compatible:
  assert(msr.mextent == m.msr.mextent);
  assert(mbr.mextent == m.mbr.mextent);
  assert(mpr.mextent == m.mpr.mextent);
  assert(mrr.mextent == m.mrr.mextent);
  assert(mcr.mextent == m.mcr.mextent);

  copy(m.begin(), m.end(), begin());
  return *this;
}

Tensor5View& Tensor5View::operator=(const Tensor5View& m) {
  // Check that sizes are compatible:
  assert(msr.mextent == m.msr.mextent);
  assert(mbr.mextent == m.mbr.mextent);
  assert(mpr.mextent == m.mpr.mextent);
  assert(mrr.mextent == m.mrr.mextent);
  assert(mcr.mextent == m.mcr.mextent);

  copy(m.begin(), m.end(), begin());
  return *this;
}

Tensor5View& Tensor5View::operator=(const Tensor5& m) {
  // Check that sizes are compatible:
  assert(msr.mextent == m.msr.mextent);
  assert(mbr.mextent == m.mbr.mextent);
  assert(mpr.mextent == m.mpr.mextent);
  assert(mrr.mextent == m.mrr.mextent);
  assert(mcr.mextent == m.mcr.mextent);

  copy(m.begin(), m.end(), begin());
  return *this;
}

Tensor5View& Tensor5View::operator=(Numeric x) {
  copy(x, begin(), end());
  return *this;
}

// Some little helper functions:
//------------------------------

Tensor5View& Tensor5View::operator*=(Numeric x) {
  const Iterator5D es = end();
  for (Iterator5D s = begin(); s != es; ++s) {
    *s *= x;
  }
  return *this;
}

Tensor5View& Tensor5View::operator/=(Numeric x) {
  const Iterator5D es = end();
  for (Iterator5D s = begin(); s != es; ++s) {
    *s /= x;
  }
  return *this;
}

Tensor5View& Tensor5View::operator+=(Numeric x) {
  const Iterator5D es = end();
  for (Iterator5D s = begin(); s != es; ++s) {
    *s += x;
  }
  return *this;
}

Tensor5View& Tensor5View::operator-=(Numeric x) {
  const Iterator5D es = end();
  for (Iterator5D s = begin(); s != es; ++s) {
    *s -= x;
  }
  return *this;
}

Tensor5View& Tensor5View::operator*=(const ConstTensor5View& x) {
  assert(nshelves() == x.nshelves());
  assert(nbooks() == x.nbooks());
  assert(npages() == x.npages());
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstIterator5D xs = x.begin();
  Iterator5D s = begin();
  const Iterator5D es = end();
  for (; s != es; ++s, ++xs) {
    *s *= *xs;
  }
  return *this;
}

Tensor5View& Tensor5View::operator/=(const ConstTensor5View& x) {
  assert(nshelves() == x.nshelves());
  assert(nbooks() == x.nbooks());
  assert(npages() == x.npages());
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstIterator5D xs = x.begin();
  Iterator5D s = begin();
  const Iterator5D es = end();
  for (; s != es; ++s, ++xs) {
    *s /= *xs;
  }
  return *this;
}

Tensor5View& Tensor5View::operator+=(const ConstTensor5View& x) {
  assert(nshelves() == x.nshelves());
  assert(nbooks() == x.nbooks());
  assert(npages() == x.npages());
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstIterator5D xs = x.begin();
  Iterator5D s = begin();
  const Iterator5D es = end();
  for (; s != es; ++s, ++xs) {
    *s += *xs;
  }
  return *this;
}

Tensor5View& Tensor5View::operator-=(const ConstTensor5View& x) {
  assert(nshelves() == x.nshelves());
  assert(nbooks() == x.nbooks());
  assert(npages() == x.npages());
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstIterator5D xs = x.begin();
  Iterator5D s = begin();
  const Iterator5D es = end();
  for (; s != es; ++s, ++xs) {
    *s -= *xs;
  }
  return *this;
}

Tensor5View::Tensor5View(const Tensor4View& a)
    : ConstTensor5View(
          a.mdata,
          Range(0,
                1,
                a.mbr.mextent * a.mpr.mextent * a.mrr.mextent * a.mcr.mextent),
          a.mbr,
          a.mpr,
          a.mrr,
          a.mcr) {
  // Nothing to do here.
}

Tensor5View::Tensor5View(Numeric* data,
                         const Range& sr,
                         const Range& br,
                         const Range& pr,
                         const Range& rr,
                         const Range& cr)
    : ConstTensor5View(data, sr, br, pr, rr, cr) {
  // Nothing to do here.
}

Tensor5View::Tensor5View(Numeric* data,
                         const Range& ps,
                         const Range& pb,
                         const Range& pp,
                         const Range& pr,
                         const Range& pc,
                         const Range& ns,
                         const Range& nb,
                         const Range& np,
                         const Range& nr,
                         const Range& nc)
    : ConstTensor5View(data, ps, pb, pp, pr, pc, ns, nb, np, nr, nc) {
  // Nothing to do here.
}

void copy(ConstIterator5D origin,
          const ConstIterator5D& end,
          Iterator5D target) {
  for (; origin != end; ++origin, ++target) {
    // We use the copy function for the next smaller rank of tensor
    // recursively:
    copy(origin->begin(), origin->end(), target->begin());
  }
}

void copy(Numeric x, Iterator5D target, const Iterator5D& end) {
  for (; target != end; ++target) {
    // We use the copy function for the next smaller rank of tensor
    // recursively:
    copy(x, target->begin(), target->end());
  }
}

// Functions for Tensor5:
// ---------------------

Tensor5::Tensor5(Index s, Index b, Index p, Index r, Index c)
    : Tensor5View(new Numeric[s * b * p * r * c],
                  Range(0, s, b * p * r * c),
                  Range(0, b, p * r * c),
                  Range(0, p, r * c),
                  Range(0, r, c),
                  Range(0, c)) {
  // Nothing to do here.
}

Tensor5::Tensor5(Index s, Index b, Index p, Index r, Index c, Numeric fill)
    : Tensor5View(new Numeric[s * b * p * r * c],
                  Range(0, s, b * p * r * c),
                  Range(0, b, p * r * c),
                  Range(0, p, r * c),
                  Range(0, r, c),
                  Range(0, c)) {
  // Here we can access the raw memory directly, for slightly
  // increased efficiency:
  std::fill_n(mdata, s * b * p * r * c, fill);
}

Tensor5::Tensor5(const ConstTensor5View& m)
    : Tensor5View(
          new Numeric[m.nshelves() * m.nbooks() * m.npages() * m.nrows() *
                      m.ncols()],
          Range(
              0, m.nshelves(), m.nbooks() * m.npages() * m.nrows() * m.ncols()),
          Range(0, m.nbooks(), m.npages() * m.nrows() * m.ncols()),
          Range(0, m.npages(), m.nrows() * m.ncols()),
          Range(0, m.nrows(), m.ncols()),
          Range(0, m.ncols())) {
  copy(m.begin(), m.end(), begin());
}

Tensor5::Tensor5(const Tensor5& m)
    : Tensor5View(
          new Numeric[m.nshelves() * m.nbooks() * m.npages() * m.nrows() *
                      m.ncols()],
          Range(
              0, m.nshelves(), m.nbooks() * m.npages() * m.nrows() * m.ncols()),
          Range(0, m.nbooks(), m.npages() * m.nrows() * m.ncols()),
          Range(0, m.npages(), m.nrows() * m.ncols()),
          Range(0, m.nrows(), m.ncols()),
          Range(0, m.ncols())) {
  // There is a catch here: If m is an empty tensor, then it will have
  // dimensions of size 0. But these are used to initialize the stride
  // for higher dimensions! Thus, this method has to be consistent
  // with the behaviour of Range::Range. For now, Range::Range allows
  // also stride 0.
  std::memcpy(
      mdata,
      m.mdata,
      nshelves() * nbooks() * npages() * nrows() * ncols() * sizeof(Numeric));
}


Tensor5& Tensor5::operator=(const Tensor5& x) {
  if (this != &x) {
    resize(x.nshelves(), x.nbooks(), x.npages(), x.nrows(), x.ncols());
    std::memcpy(
        mdata,
        x.mdata,
        nshelves() * nbooks() * npages() * nrows() * ncols() * sizeof(Numeric));
  }
  return *this;
}

Tensor5& Tensor5::operator=(Tensor5&& x) noexcept {
  if (this != &x) {
    delete[] mdata;
    mdata = x.mdata;
    msr = x.msr;
    mbr = x.mbr;
    mpr = x.mpr;
    mrr = x.mrr;
    mcr = x.mcr;
    x.msr = Range(0, 0);
    x.mbr = Range(0, 0);
    x.mpr = Range(0, 0);
    x.mrr = Range(0, 0);
    x.mcr = Range(0, 0);
    x.mdata = nullptr;
  }
  return *this;
}

Tensor5& Tensor5::operator=(Numeric x) {
  std::fill_n(mdata, nshelves() * nbooks() * npages() * nrows() * ncols(), x);
  return *this;
}

void Tensor5::resize(Index s, Index b, Index p, Index r, Index c) {
  assert(0 <= s);
  assert(0 <= b);
  assert(0 <= p);
  assert(0 <= r);
  assert(0 <= c);

  if (msr.mextent != s || mbr.mextent != b || mpr.mextent != p ||
      mrr.mextent != r || mcr.mextent != c) {
    delete[] mdata;
    mdata = new Numeric[s * b * p * r * c];

    msr.mstart = 0;
    msr.mextent = s;
    msr.mstride = b * p * r * c;

    mbr.mstart = 0;
    mbr.mextent = b;
    mbr.mstride = p * r * c;

    mpr.mstart = 0;
    mpr.mextent = p;
    mpr.mstride = r * c;

    mrr.mstart = 0;
    mrr.mextent = r;
    mrr.mstride = c;

    mcr.mstart = 0;
    mcr.mextent = c;
    mcr.mstride = 1;
  }
}

void swap(Tensor5& t1, Tensor5& t2) {
  std::swap(t1.msr, t2.msr);
  std::swap(t1.mbr, t2.mbr);
  std::swap(t1.mpr, t2.mpr);
  std::swap(t1.mrr, t2.mrr);
  std::swap(t1.mcr, t2.mcr);
  std::swap(t1.mdata, t2.mdata);
}

Tensor5::~Tensor5() {
  //   cout << "Destroying a Tensor5:\n"
  //        << *this << "\n........................................\n";
  delete[] mdata;
}

void transform(Tensor5View y, double (&my_func)(double), ConstTensor5View x) {
  // Check dimensions:
  assert(y.nshelves() == x.nshelves());
  assert(y.nbooks() == x.nbooks());
  assert(y.npages() == x.npages());
  assert(y.nrows() == x.nrows());
  assert(y.ncols() == x.ncols());

  const ConstIterator5D xe = x.end();
  ConstIterator5D xi = x.begin();
  Iterator5D yi = y.begin();
  for (; xi != xe; ++xi, ++yi) {
    // Use the transform function of lower dimensional tensors
    // recursively:
    transform(*yi, my_func, *xi);
  }
}

Numeric max(const ConstTensor5View& x) {
  const ConstIterator5D xe = x.end();
  ConstIterator5D xi = x.begin();

  // Initial value for max:
  Numeric themax = max(*xi);
  ++xi;

  for (; xi != xe; ++xi) {
    // Use the max function of lower dimensional tensors
    // recursively:
    Numeric maxi = max(*xi);
    if (maxi > themax) themax = maxi;
  }

  return themax;
}

Numeric min(const ConstTensor5View& x) {
  const ConstIterator5D xe = x.end();
  ConstIterator5D xi = x.begin();

  // Initial value for min:
  Numeric themin = min(*xi);
  ++xi;

  for (; xi != xe; ++xi) {
    // Use the min function of lower dimensional tensors
    // recursively:
    Numeric mini = min(*xi);
    if (mini < themin) themin = mini;
  }

  return themin;
}

// Helper function for debugging
#ifndef NDEBUG

Numeric debug_tensor5view_get_elem(
    Tensor5View& tv, Index s, Index b, Index p, Index r, Index c) {
  return tv(s, b, p, r, c);
}

#endif