ARTS  2.3.1285(git:92a29ea9-dirty)
matpackIII.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
25 #include "matpackIII.h"
26 #include "exceptions.h"
27 
28 using std::runtime_error;
29 
30 // Functions for ConstTensor3View:
31 // ------------------------------
32 
34 
39  return (npages() == 0 || nrows() == 0 || ncols() == 0);
40 }
41 
46  const Range& r,
47  const Range& c) const {
48  return ConstTensor3View(mdata, mpr, mrr, mcr, p, r, c);
49 }
50 
54  const Range& r,
55  Index c) const {
56  // Check that c is valid:
57  assert(0 <= c);
58  assert(c < mcr.mextent);
59 
60  return ConstMatrixView(mdata + mcr.mstart + c * mcr.mstride, mpr, mrr, p, r);
61 }
62 
66  Index r,
67  const Range& c) const {
68  // Check that r is valid:
69  assert(0 <= r);
70  assert(r < mrr.mextent);
71 
72  return ConstMatrixView(mdata + mrr.mstart + r * mrr.mstride, mpr, mcr, p, c);
73 }
74 
78  const Range& r,
79  const Range& c) const {
80  // Check that p is valid:
81  assert(0 <= p);
82  assert(p < mpr.mextent);
83 
84  return ConstMatrixView(mdata + mpr.mstart + p * mpr.mstride, mrr, mcr, r, c);
85 }
86 
90  Index r,
91  const Range& c) const {
92  // Check that p and r are valid:
93  assert(0 <= p);
94  assert(0 <= r);
95  assert(p < mpr.mextent);
96  assert(r < mrr.mextent);
97 
98  return ConstVectorView(
99  mdata + mpr.mstart + p * mpr.mstride + mrr.mstart + r * mrr.mstride,
100  mcr,
101  c);
102 }
103 
107  const Range& r,
108  Index c) const {
109  // Check that p and c are valid:
110  assert(0 <= p);
111  assert(0 <= c);
112  assert(p < mpr.mextent);
113  assert(c < mcr.mextent);
114 
115  return ConstVectorView(
116  mdata + mpr.mstart + p * mpr.mstride + mcr.mstart + c * mcr.mstride,
117  mrr,
118  r);
119 }
120 
124  Index r,
125  Index c) const {
126  // Check that r and c are valid:
127  assert(0 <= r);
128  assert(0 <= c);
129  assert(r < mrr.mextent);
130  assert(c < mcr.mextent);
131 
132  return ConstVectorView(
133  mdata + mrr.mstart + r * mrr.mstride + mcr.mstart + c * mcr.mstride,
134  mpr,
135  p);
136 }
137 
141  mpr.mstride);
142 }
143 
146  return ConstIterator3D(
148  mdata + mpr.mstart + (mpr.mextent) * mpr.mstride, mrr, mcr),
149  mpr.mstride);
150 }
151 
154  : mpr(0, 1, a.mrr.mextent * a.mcr.mextent),
155  mrr(a.mrr),
156  mcr(a.mcr),
157  mdata(a.mdata) {
158  // Nothing to do here.
159 }
160 
166  const Range& pr,
167  const Range& rr,
168  const Range& cr)
169  : mpr(pr), mrr(rr), mcr(cr), mdata(data) {
170  // Nothing to do here.
171 }
172 
181  const Range& pp,
182  const Range& pr,
183  const Range& pc,
184  const Range& np,
185  const Range& nr,
186  const Range& nc)
187  : mpr(pp, np), mrr(pr, nr), mcr(pc, nc), mdata(data) {
188  // Nothing to do here.
189 }
190 
194 std::ostream& operator<<(std::ostream& os, const ConstTensor3View& v) {
195  // Page iterators:
196  ConstIterator3D ip = v.begin();
197  const ConstIterator3D end_page = v.end();
198 
199  if (ip != end_page) {
200  os << *ip;
201  ++ip;
202  }
203 
204  for (; ip != end_page; ++ip) {
205  os << "\n\n";
206  os << *ip;
207  }
208 
209  return os;
210 }
211 
212 // Functions for Tensor3View:
213 // -------------------------
214 
219  const Range& r,
220  const Range& c) {
221  return Tensor3View(mdata, mpr, mrr, mcr, p, r, c);
222 }
223 
227  // Check that c is valid:
228  assert(0 <= c);
229  assert(c < mcr.mextent);
230 
231  return MatrixView(mdata + mcr.mstart + c * mcr.mstride, mpr, mrr, p, r);
232 }
233 
237  // Check that r is valid:
238  assert(0 <= r);
239  assert(r < mrr.mextent);
240 
241  return MatrixView(mdata + mrr.mstart + r * mrr.mstride, mpr, mcr, p, c);
242 }
243 
247  // Check that p is valid:
248  assert(0 <= p);
249  assert(p < mpr.mextent);
250 
251  return MatrixView(mdata + mpr.mstart + p * mpr.mstride, mrr, mcr, r, c);
252 }
253 
257  // Check that p and r are valid:
258  assert(0 <= p);
259  assert(0 <= r);
260  assert(p < mpr.mextent);
261  assert(r < mrr.mextent);
262 
263  return VectorView(
264  mdata + mpr.mstart + p * mpr.mstride + mrr.mstart + r * mrr.mstride,
265  mcr,
266  c);
267 }
268 
272  // Check that p and c are valid:
273  assert(0 <= p);
274  assert(0 <= c);
275  assert(p < mpr.mextent);
276  assert(c < mcr.mextent);
277 
278  return VectorView(
279  mdata + mpr.mstart + p * mpr.mstride + mcr.mstart + c * mcr.mstride,
280  mrr,
281  r);
282 }
283 
287  // Check that r and r are valid:
288  assert(0 <= r);
289  assert(0 <= c);
290  assert(r < mrr.mextent);
291  assert(c < mcr.mextent);
292 
293  return VectorView(
294  mdata + mrr.mstart + r * mrr.mstride + mcr.mstart + c * mcr.mstride,
295  mpr,
296  p);
297 }
298 
306  if (mpr.mstart != 0 || mpr.mstride != mrr.mextent * mcr.mextent ||
307  mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
308  mcr.mstride != 1)
309  throw std::runtime_error(
310  "A Tensor3View can only be converted to a plain C-array if it's pointing to a continuous block of data");
311 
312  return mdata;
313 }
314 
322  if (mpr.mstart != 0 || mpr.mstride != mrr.mextent * mcr.mextent ||
323  mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
324  mcr.mstride != 1)
325  throw std::runtime_error(
326  "A Tensor3View can only be converted to a plain C-array if it's pointing to a continuous block of data");
327 
328  return mdata;
329 }
330 
334 }
335 
338  return Iterator3D(
340  mpr.mstride);
341 }
342 
348  // Check that sizes are compatible:
349  assert(mpr.mextent == m.mpr.mextent);
350  assert(mrr.mextent == m.mrr.mextent);
351  assert(mcr.mextent == m.mcr.mextent);
352 
353  copy(m.begin(), m.end(), begin());
354  return *this;
355 }
356 
363  // Check that sizes are compatible:
364  assert(mpr.mextent == m.mpr.mextent);
365  assert(mrr.mextent == m.mrr.mextent);
366  assert(mcr.mextent == m.mcr.mextent);
367 
368  copy(m.begin(), m.end(), begin());
369  return *this;
370 }
371 
376  // Check that sizes are compatible:
377  assert(mpr.mextent == m.mpr.mextent);
378  assert(mrr.mextent == m.mrr.mextent);
379  assert(mcr.mextent == m.mcr.mextent);
380 
381  copy(m.begin(), m.end(), begin());
382  return *this;
383 }
384 
388  copy(x, begin(), end());
389  return *this;
390 }
391 
392 // Some little helper functions:
393 //------------------------------
394 
395 Numeric add(Numeric x, Numeric y) { return x + y; }
396 
399  const Iterator3D ep = end();
400  for (Iterator3D p = begin(); p != ep; ++p) {
401  *p *= x;
402  }
403  return *this;
404 }
405 
408  const Iterator3D ep = end();
409  for (Iterator3D p = begin(); p != ep; ++p) {
410  *p /= x;
411  }
412  return *this;
413 }
414 
417  const Iterator3D ep = end();
418  for (Iterator3D p = begin(); p != ep; ++p) {
419  *p += x;
420  }
421  return *this;
422 }
423 
426  const Iterator3D ep = end();
427  for (Iterator3D p = begin(); p != ep; ++p) {
428  *p -= x;
429  }
430  return *this;
431 }
432 
435  assert(npages() == x.npages());
436  assert(nrows() == x.nrows());
437  assert(ncols() == x.ncols());
438  ConstIterator3D xp = x.begin();
439  Iterator3D p = begin();
440  const Iterator3D ep = end();
441  for (; p != ep; ++p, ++xp) {
442  *p *= *xp;
443  }
444  return *this;
445 }
446 
449  assert(npages() == x.npages());
450  assert(nrows() == x.nrows());
451  assert(ncols() == x.ncols());
452  ConstIterator3D xp = x.begin();
453  Iterator3D p = begin();
454  const Iterator3D ep = end();
455  for (; p != ep; ++p, ++xp) {
456  *p /= *xp;
457  }
458  return *this;
459 }
460 
463  assert(npages() == x.npages());
464  assert(nrows() == x.nrows());
465  assert(ncols() == x.ncols());
466  ConstIterator3D xp = x.begin();
467  Iterator3D p = begin();
468  const Iterator3D ep = end();
469  for (; p != ep; ++p, ++xp) {
470  *p += *xp;
471  }
472  return *this;
473 }
474 
477  assert(npages() == x.npages());
478  assert(nrows() == x.nrows());
479  assert(ncols() == x.ncols());
480  ConstIterator3D xp = x.begin();
481  Iterator3D p = begin();
482  const Iterator3D ep = end();
483  for (; p != ep; ++p, ++xp) {
484  *p -= *xp;
485  }
486  return *this;
487 }
488 
492  a.mdata, Range(0, 1, a.mrr.mextent * a.mcr.mextent), a.mrr, a.mcr) {
493  // Nothing to do here.
494 }
495 
500  const Range& pr,
501  const Range& rr,
502  const Range& cr)
503  : ConstTensor3View(data, pr, rr, cr) {
504  // Nothing to do here.
505 }
506 
524  const Range& pp,
525  const Range& pr,
526  const Range& pc,
527  const Range& np,
528  const Range& nr,
529  const Range& nc)
530  : ConstTensor3View(data, pp, pr, pc, np, nr, nc) {
531  // Nothing to do here.
532 }
533 
538 void copy(ConstIterator3D origin,
539  const ConstIterator3D& end,
540  Iterator3D target) {
541  for (; origin != end; ++origin, ++target) {
542  // We use the copy function for the next smaller rank of tensor
543  // recursively:
544  copy(origin->begin(), origin->end(), target->begin());
545  }
546 }
547 
549 void copy(Numeric x, Iterator3D target, const Iterator3D& end) {
550  for (; target != end; ++target) {
551  // We use the copy function for the next smaller rank of tensor
552  // recursively:
553  copy(x, target->begin(), target->end());
554  }
555 }
556 
557 // Functions for Tensor3:
558 // ---------------------
559 
563  : Tensor3View(new Numeric[p * r * c],
564  Range(0, p, r * c),
565  Range(0, r, c),
566  Range(0, c)) {
567  // Nothing to do here.
568 }
569 
572  : Tensor3View(new Numeric[p * r * c],
573  Range(0, p, r * c),
574  Range(0, r, c),
575  Range(0, c)) {
576  // Here we can access the raw memory directly, for slightly
577  // increased efficiency:
578  std::fill_n(mdata, p * r * c, fill);
579 }
580 
584  : Tensor3View(new Numeric[m.npages() * m.nrows() * m.ncols()],
585  Range(0, m.npages(), m.nrows() * m.ncols()),
586  Range(0, m.nrows(), m.ncols()),
587  Range(0, m.ncols())) {
588  copy(m.begin(), m.end(), begin());
589 }
590 
594  : Tensor3View(new Numeric[m.npages() * m.nrows() * m.ncols()],
595  Range(0, m.npages(), m.nrows() * m.ncols()),
596  Range(0, m.nrows(), m.ncols()),
597  Range(0, m.ncols())) {
598  // There is a catch here: If m is an empty tensor, then it will have
599  // dimensions of size 0. But these are used to initialize the stride
600  // for higher dimensions! Thus, this method has to be consistent
601  // with the behaviour of Range::Range. For now, Range::Range allows
602  // also stride 0.
603  std::memcpy(mdata, m.mdata, npages() * nrows() * ncols() * sizeof(Numeric));
604 }
605 
607 
631  if (this != &x) {
632  resize(x.npages(), x.nrows(), x.ncols());
633  std::memcpy(mdata, x.mdata, npages() * nrows() * ncols() * sizeof(Numeric));
634  }
635  return *this;
636 }
637 
640  if (this != &x) {
641  delete[] mdata;
642  mdata = x.mdata;
643  mpr = x.mpr;
644  mrr = x.mrr;
645  mcr = x.mcr;
646  x.mpr = Range(0, 0);
647  x.mrr = Range(0, 0);
648  x.mcr = Range(0, 0);
649  x.mdata = nullptr;
650  }
651  return *this;
652 }
653 
657  std::fill_n(mdata, npages() * nrows() * ncols(), x);
658  return *this;
659 }
660 
665  assert(0 <= p);
666  assert(0 <= r);
667  assert(0 <= c);
668 
669  if (mpr.mextent != p || mrr.mextent != r || mcr.mextent != c) {
670  delete[] mdata;
671  mdata = new Numeric[p * r * c];
672 
673  mpr.mstart = 0;
674  mpr.mextent = p;
675  mpr.mstride = r * c;
676 
677  mrr.mstart = 0;
678  mrr.mextent = r;
679  mrr.mstride = c;
680 
681  mcr.mstart = 0;
682  mcr.mextent = c;
683  mcr.mstride = 1;
684  }
685 }
686 
688 void swap(Tensor3& t1, Tensor3& t2) {
689  std::swap(t1.mpr, t2.mpr);
690  std::swap(t1.mrr, t2.mrr);
691  std::swap(t1.mcr, t2.mcr);
692  std::swap(t1.mdata, t2.mdata);
693 }
694 
698  // cout << "Destroying a Tensor3:\n"
699  // << *this << "\n........................................\n";
700  delete[] mdata;
701 }
702 
718 void transform(Tensor3View y, double (&my_func)(double), ConstTensor3View x) {
719  // Check dimensions:
720  assert(y.npages() == x.npages());
721  assert(y.nrows() == x.nrows());
722  assert(y.ncols() == x.ncols());
723 
724  const ConstIterator3D xe = x.end();
725  ConstIterator3D xi = x.begin();
726  Iterator3D yi = y.begin();
727  for (; xi != xe; ++xi, ++yi) {
728  // Use the transform function of lower dimensional tensors
729  // recursively:
730  transform(*yi, my_func, *xi);
731  }
732 }
733 
736  const ConstIterator3D xe = x.end();
737  ConstIterator3D xi = x.begin();
738 
739  // Initial value for max:
740  Numeric themax = max(*xi);
741  ++xi;
742 
743  for (; xi != xe; ++xi) {
744  // Use the max function of lower dimensional tensors
745  // recursively:
746  Numeric maxi = max(*xi);
747  if (maxi > themax) themax = maxi;
748  }
749 
750  return themax;
751 }
752 
755  const ConstIterator3D xe = x.end();
756  ConstIterator3D xi = x.begin();
757 
758  // Initial value for min:
759  Numeric themin = min(*xi);
760  ++xi;
761 
762  for (; xi != xe; ++xi) {
763  // Use the min function of lower dimensional tensors
764  // recursively:
765  Numeric mini = min(*xi);
766  if (mini < themin) themin = mini;
767  }
768 
769  return themin;
770 }
771 
773 // Helper function for debugging
774 #ifndef NDEBUG
775 
791  return tv(p, r, c);
792 }
793 
794 #endif
795 
798 
814  assert(A.npages() == B.nelem());
815  assert(A.nrows() == C.nrows());
816  assert(A.ncols() == C.ncols());
817 
818  for (Index ii = 0; ii < B.nelem(); ii++) {
819  A(ii, joker, joker) = C;
820  A(ii, joker, joker) *= B[ii];
821  }
822 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Iterator2D begin()
Return const iterator to first row.
Definition: matpackI.cc:618
The VectorView class.
Definition: matpackI.h:610
Tensor3View operator()(const Range &p, const Range &r, const Range &c)
Index operator for subrange.
Definition: matpackIII.cc:218
Index mstride
The stride.
Definition: matpackI.h:355
void swap(Tensor3 &t1, Tensor3 &t2)
Swaps two objects.
Definition: matpackIII.cc:688
Index mstart
The start index.
Definition: matpackI.h:346
Tensor3()=default
Range mcr
The column range of mdata that is actually used.
Definition: matpackIII.h:225
The MatrixView class.
Definition: matpackI.h:1093
The range class.
Definition: matpackI.h:160
void copy(ConstIterator3D origin, const ConstIterator3D &end, Iterator3D target)
Copy data between begin and end to target.
Definition: matpackIII.cc:538
bool empty() const
Check if variable is empty.
Definition: matpackIII.cc:38
Iterator2D end()
Return iterator behind last row.
Definition: matpackI.cc:623
Range mpr
The page range of mdata that is actually used.
Definition: matpackIII.h:221
Numeric max(const ConstTensor3View &x)
Max function, tensor version.
Definition: matpackIII.cc:735
Iterator3D end()
Return iterator behind last page.
Definition: matpackIII.cc:337
ConstTensor3View()=default
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
Tensor3View & operator-=(Numeric x)
Subtraction of scalar.
Definition: matpackIII.cc:425
G0 G2 FVC Y DV Numeric Numeric Numeric Zeeman LowerQuantumNumbers void * data
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
virtual ~Tensor3()
Destructor for Tensor3.
Definition: matpackIII.cc:697
ConstIterator2D begin() const
Return const iterator to first row.
Definition: matpackI.cc:469
Tensor3View()=default
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
The Tensor3 class.
Definition: matpackIII.h:339
ConstIterator3D end() const
Return const iterator behind last page.
Definition: matpackIII.cc:145
Tensor3View & operator/=(Numeric x)
Division by scalar.
Definition: matpackIII.cc:407
friend class Tensor3View
Definition: matpackIII.h:193
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
Iterator3D begin()
Return iterator to first page.
Definition: matpackIII.cc:332
The Tensor3View class.
Definition: matpackIII.h:239
Tensor3View & operator*=(Numeric x)
Multiplication by scalar.
Definition: matpackIII.cc:398
The declarations of all the exception classes.
ConstIterator2D end() const
Return const iterator behind last row.
Definition: matpackI.cc:474
Numeric debug_tensor3view_get_elem(Tensor3View &tv, Index p, Index r, Index c)
Helper function to access tensor elements.
Definition: matpackIII.cc:790
void transform(Tensor3View y, double(&my_func)(double), ConstTensor3View x)
A generic transform function for tensors, which can be used to implement mathematical functions opera...
Definition: matpackIII.cc:718
const Joker joker
friend void swap(Tensor3 &t1, Tensor3 &t2)
Swaps two objects.
Definition: matpackIII.cc:688
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Index mextent
The number of elements.
Definition: matpackI.h:353
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:664
Numeric min(const ConstTensor3View &x)
Min function, tensor version.
Definition: matpackIII.cc:754
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Tensor3 & operator=(const Tensor3 &x)
Assignment operator from another tensor.
Definition: matpackIII.cc:630
ConstIterator3D begin() const
Return const iterator to first page.
Definition: matpackIII.cc:139
A constant view of a Tensor3.
Definition: matpackIII.h:132
A constant view of a Vector.
Definition: matpackI.h:476
Tensor3View & operator+=(Numeric x)
Addition of scalar.
Definition: matpackIII.cc:416
ConstTensor3View operator()(const Range &p, const Range &r, const Range &c) const
Const index operator for subrange.
Definition: matpackIII.cc:45
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackIII.h:227
A constant view of a Matrix.
Definition: matpackI.h:982
Implementation of Tensors of Rank 3.
Definition: matpackIII.h:34
Tensor3View & operator=(const ConstTensor3View &v)
Assignment operator.
Definition: matpackIII.cc:347
Range mrr
The row range of mdata that is actually used.
Definition: matpackIII.h:223
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackIII.cc:321
std::ostream & operator<<(std::ostream &os, const ConstTensor3View &v)
Output operator.
Definition: matpackIII.cc:194
Const version of Iterator3D.
Definition: matpackIII.h:76
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
void mult(Tensor3View A, const ConstVectorView B, const ConstMatrixView C)
mult Tensor3
Definition: matpackIII.cc:813
Numeric add(Numeric x, Numeric y)
Definition: matpackIII.cc:395