tensor.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_TENSOR_H
24 #define O2SCL_TENSOR_H
25 
26 /** \file tensor.h
27  \brief File defining \ref o2scl::tensor and rank-specific children
28 */
29 #include <iostream>
30 #include <cstdlib>
31 #include <string>
32 #include <fstream>
33 #include <sstream>
34 
35 #include <boost/numeric/ublas/vector.hpp>
36 #include <boost/numeric/ublas/vector_proxy.hpp>
37 #include <boost/numeric/ublas/matrix.hpp>
38 #include <boost/numeric/ublas/matrix_proxy.hpp>
39 
40 #include <gsl/gsl_matrix.h>
41 #include <gsl/gsl_ieee_utils.h>
42 
43 #include <o2scl/err_hnd.h>
44 #include <o2scl/interp.h>
45 
46 #ifndef DOXYGEN_NO_O2NS
47 namespace o2scl {
48 #endif
49 
50  /** \brief Tensor class with arbitrary dimensions
51 
52  The elements of a tensor are typically specified as a list of
53  <tt>size_t</tt> numbers with length equal to the tensor rank.
54  For a rank-4 tensor named \c t, the element
55  <tt>t[1][2][0][3]</tt> can be obtained with something similar to
56  \code
57  size_t ix[4]={1,2,0,3};
58  double x=t.get(ix);
59  \endcode
60 
61  Empty tensors have zero rank.
62 
63  The type <tt>vec_t</tt> can be any vector type with
64  <tt>operator[]</tt>, <tt>size()</tt> and <tt>resize()</tt>
65  methods. The type <tt>vec_size_t</tt> can be any integer-like
66  vector type with <tt>operator[]</tt>, <tt>size()</tt> and
67  <tt>resize()</tt> methods.
68 
69  For I/O with tensors, see \ref o2scl_hdf::hdf_file::setd_ten()
70  and \ref o2scl_hdf::hdf_file::getd_ten() . See also
71  the discussion in the sections \ref tensor_subsect and
72  \ref vec_io_cont_subsect of the user's guide.
73 
74  The storage pattern is a generalization of row-major order.
75  In the case of a 4-rank tensor, the location of a generic
76  element is
77  \f[
78  \left( \left( i_0 s_1 + i_1 \right) s_2 + i_2 \right) s_3 + i_3 \, .
79  \f]
80  In this case the distance between two elements \f$(i_0,i_1,
81  i_2,i_3)\f$ and \f$ (i_0+1,i_1,i_2,i_3) \f$ is \f$ s_1 s_2 s_3
82  \f$,the distance between two elements \f$(i_0,i_1,i_2,
83  i_3)\f$ and \f$ (i_0,i_1+1,i_2,i_3) \f$ is \f$ s_2 s_3 \f$, and
84  the distance between two elements \f$(i_0,i_1,i_2,i_3)\f$ and
85  \f$ (i_0,i_1,i_2,i_3+1) \f$ is just unity.
86 
87  \note Slices of tensors are subsets obtained from fixing the
88  index of several dimensions, while letting other dimensions
89  vary. For a slice with one dimension not fixed, see \ref
90  vector_slice(). The \ref o2scl::tensor::vector_slice() function
91  should clearly work for uBlas vectors, and seems to work with
92  std::vector objects also, but latter use has not been fully
93  tested.
94 
95  \future Create an operator[] for tensor and not just tensor1?
96 
97  \future Could implement arithmetic operators + and - and some
98  different products.
99 
100  \future Implement copies to and from vector
101  and matrices
102 
103  \future Implement tensor contractions, i.e. tensor
104  = tensor * tensor
105 
106  \future Consider making a template type to replace double,
107  e.g. <tt>value_t</tt>.
108 
109  \future Could be interesting to write an iterator for this class.
110 
111  \future Create an is_valid() function which checks sizes
112  \comment
113  A tensor is defined to be empty if it's rank is zero. A tensor's
114  rank should be zero if and only if no memory has been allocated
115  for it.
116  \endcomment
117 
118  */
119  template<class vec_t=std::vector<double>,
120  class vec_size_t=std::vector<size_t> > class tensor {
121 
122  public:
123 
124 #ifndef DOXYGEN_INTERNAL
125 
126  protected:
127 
128  /// The data
129  vec_t data;
130 
131  /// A rank-sized vector of the sizes of each dimension
132  vec_size_t size;
133 
134  /// Rank
135  size_t rk;
136 
137 #endif
138 
139  public:
140 
141  /// Create an empty tensor with zero rank
142  tensor() {
143  rk=0;
144  }
145 
146  /** \brief Create a tensor of rank \c rank with sizes given in \c dim
147 
148  The parameter \c dim must be a pointer to a vector of sizes with
149  length \c rank. If the user requests any of the sizes to be
150  zero, this constructor will call the error handler, create an
151  empty tensor, and will allocate no memory.
152  */
153  template<class size_vec_t>
154  tensor(size_t rank, const size_vec_t &dim) {
155  if (rank==0) {
156  rk=0;
157  } else {
158  rk=rank;
159  for(size_t i=0;i<rk;i++) {
160  if (dim[i]==0) {
161  rk=0;
162  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
163  "rank for index "+szttos(i)+
164  " in tensor::tensor(size_t,size_vec_t)").c_str(),
165  exc_einval);
166  }
167  }
168  size.resize(rk);
169  size_t tot=1;
170  for(size_t i=0;i<rk;i++) {
171  size[i]=dim[i];
172  tot*=size[i];
173  }
174  data.resize(tot);
175  }
176  }
177 
178  ~tensor() {
179  }
180 
181  /// \name Clear method
182  //@{
183  /// Clear the tensor of all data and free allocated memory
184  void clear() {
185  rk=0;
186  data.resize(0);
187  size.resize(0);
188  return;
189  }
190  //@}
191 
192  /// \name Set functions
193  //@{
194  /// Set the element indexed by \c index to value \c val
195  template<class size_vec_t>
196  void set(const size_vec_t &index, double val) {
197 #if O2SCL_NO_RANGE_CHECK
198 #else
199  if (rk==0) {
200  O2SCL_ERR("Empty tensor in tensor::set().",exc_einval);
201  }
202  if (index[0]>=size[0]) {
203  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
204  " greater than or equal to size[0]="+
205  szttos(size[0])+" in tensor::set().").c_str(),
206  exc_eindex);
207  }
208 #endif
209  size_t ix=index[0];
210  for(size_t i=1;i<rk;i++) {
211 #if O2SCL_NO_RANGE_CHECK
212 #else
213  if (index[i]>=size[i]) {
214  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
215  szttos(index[i])+" greater than or equal to size "+
216  szttos(size[i])+" in tensor::set().").c_str(),
217  exc_eindex);
218  }
219 #endif
220  ix*=size[i];
221  ix+=index[i];
222  }
223  data[ix]=val;
224  return;
225  }
226 
227  /// Set all elements in a tensor to some fixed value
228  void set_all(double x) {
229  for(size_t i=0;i<total_size();i++) data[i]=x;
230  return;
231  }
232  //@}
233 
234  /// \name Get functions
235  //@{
236  /// Get the element indexed by \c index
237  template<class size_vec_t> double &get(const size_vec_t &index) {
238 #if O2SCL_NO_RANGE_CHECK
239 #else
240  if (rk==0) {
241  O2SCL_ERR("Empty tensor in tensor::get().",exc_einval);
242  }
243  if (index[0]>=size[0]) {
244  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
245  " greater than or equal to size[0]="+
246  szttos(size[0])+" in tensor::get().").c_str(),
247  exc_eindex);
248  }
249 #endif
250  size_t ix=index[0];
251  for(size_t i=1;i<rk;i++) {
252 #if O2SCL_NO_RANGE_CHECK
253 #else
254  if (index[i]>=size[i]) {
255  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
256  szttos(index[i])+" greater than or equal to size "+
257  szttos(size[i])+" in tensor::get().").c_str(),
258  exc_eindex);
259  }
260 #endif
261  ix*=size[i];
262  ix+=index[i];
263  }
264  return data[ix];
265  }
266 
267  /// Get a const reference to the element indexed by \c index
268  template<class size_vec_t>
269  double const &get(const size_vec_t &index) const {
270 #if O2SCL_NO_RANGE_CHECK
271 #else
272  if (rk==0) {
273  O2SCL_ERR("Empty tensor in tensor::get() (const).",exc_einval);
274  }
275  if (index[0]>=size[0]) {
276  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
277  " greater than or equal to size[0]="+
278  szttos(size[0])+" in tensor::get() (const).").c_str(),
279  exc_eindex);
280  }
281 #endif
282  size_t ix=index[0];
283  for(size_t i=1;i<rk;i++) {
284 #if O2SCL_NO_RANGE_CHECK
285 #else
286  if (index[i]>=size[i]) {
287  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
288  szttos(index[i])+" greater than or equal to size "+
289  szttos(size[i])+" in tensor::get() (const).").c_str(),
290  exc_eindex);
291  }
292 #endif
293  ix*=size[i];
294  ix+=index[i];
295  }
296  return data[ix];
297  }
298  //@}
299 
300  typedef boost::numeric::ublas::vector_slice<
301  boost::numeric::ublas::vector<double> > ubvector_slice;
302  typedef boost::numeric::ublas::slice slice;
303 
304  /// \name Slice function
305  //@{
306  /** \brief Fix all but one index to create a vector
307 
308  This fixes all of the indices to the values given in \c index
309  except for the index number \c ix, and returns the
310  corresponding vector, whose length is equal to the size of the
311  tensor in that index. The value <tt>index[ix]</tt> is ignored.
312 
313  For example, for a rank 3 tensor allocated with
314  \code
315  tensor t;
316  size_t dim[3]={3,4,5};
317  t.resize(3,dim);
318  \endcode
319  the following code
320  \code
321  size_t index[3]={1,0,3};
322  ubvector_view v=t.vector_slice(1,index);
323  \endcode
324  Gives a vector \c v of length 4 which refers to the values
325  <tt>t(1,0,3)</tt>, <tt>t(1,1,3)</tt>, <tt>t(1,2,3)</tt>, and
326  <tt>t(1,3,3)</tt>.
327  */
328  template<class size_vec_t>
329  ubvector_slice vector_slice(size_t ix, const size_vec_t &index) {
330  if (ix+1>rk) {
331  O2SCL_ERR((((std::string)"Specified index ")+szttos(ix)+
332  " greater than or equal to rank "+szttos(rk)+
333  " in tensor::vector_slice()").c_str(),
334  exc_eindex);
335  }
336  size_t start;
337  if (ix==0) start=0;
338  else start=index[0];
339  for(size_t i=1;i<rk;i++) {
340  start*=size[i];
341  if (i!=ix) start+=index[i];
342  }
343  size_t stride=1;
344  for(size_t i=ix+1;i<rk;i++) stride*=size[i];
345  return ubvector_slice(data,slice(start,stride,size[ix]));
346  }
347  //@}
348 
349 #ifdef O2SCL_NEVER_DEFINED
350 
351  /** \brief Fix all but two indices to create a matrix
352 
353  One could use ublas::make_matrix_from_pointer() and
354  then ublas matrix slicing to perform this operation, but
355  only if the template types are built on ublas objects.
356  */
357  template<class size_vec_t> ubmatrix_array matrix_slice() {
358  }
359 
360 #endif
361 
362  /// \name Resize method
363  //@{
364  /** \brief Resize the tensor to rank \c rank with sizes
365  given in \c dim
366 
367  The parameter \c dim must be a vector of sizes with a length
368  equal to \c rank. This resize method is always destructive.
369 
370  If the user requests any of the sizes to be zero, this
371  function will call the error handler.
372  */
373  template<class size_vec_t>
374  void resize(size_t rank, const size_vec_t &dim) {
375  for(size_t i=0;i<rank;i++) {
376  if (dim[i]==0) {
377  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
378  "rank for index "+szttos(i)+" in tensor::"+
379  "resize().").c_str(),exc_einval);
380  }
381  }
382  rk=rank;
383  size.resize(rk);
384  if (rk==0) {
385  data.resize(0);
386  } else {
387  size_t tot=1;
388  for(size_t i=0;i<rk;i++) {
389  size[i]=dim[i];
390  tot*=size[i];
391  }
392  data.resize(tot);
393  }
394  return;
395  }
396  //@}
397 
398  /// \name Size functions
399  //@{
400  /// Return the rank of the tensor
401  size_t get_rank() const { return rk; }
402 
403  /// Returns the size of the ith index
404  size_t get_size(size_t i) const {
405  if (i>=rk) {
406  O2SCL_ERR((((std::string)"Specified index ")+szttos(i)+
407  " greater than or equal to rank "+szttos(rk)+
408  " in tensor::get_size()").c_str(),
409  exc_einval);
410  }
411  return size[i];
412  }
413 
414  /// Return the full vector of sizes
415  const vec_size_t &get_size_arr() const {
416  return size;
417  }
418 
419  /// Return the full data vector
420  const vec_t &get_data() const {
421  return data;
422  }
423 
424  /** \brief Returns the size of the tensor (the product of
425  the sizes over every index)
426  */
427  size_t total_size() const {
428  if (rk==0) return 0;
429  size_t tot=1;
430  for(size_t i=0;i<rk;i++) tot*=size[i];
431  return tot;
432  }
433  //@}
434 
435  /// \name Index manipulation
436  //@{
437  /// Pack the indices into a single vector index
438  template<class size_vec_t>
439  size_t pack_indices(const size_vec_t &index) {
440  if (rk==0) {
441  O2SCL_ERR("Empty tensor in tensor::pack_indices().",exc_einval);
442  return 0;
443  }
444  if (index[0]>=size[0]) {
445  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
446  " greater than or equal to size[0]="+
447  szttos(size[0])+" in tensor::pack_indices().").c_str(),
448  exc_eindex);
449  }
450  size_t ix=index[0];
451  for(size_t i=1;i<rk;i++) {
452  if (index[i]>=size[i]) {
453  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
454  szttos(index[i])+" greater than or equal to size "+
455  szttos(size[i])+" in tensor::pack_indices().").c_str(),
456  exc_eindex);
457  }
458  ix*=size[i];
459  ix+=index[i];
460  }
461  return ix;
462  }
463 
464  /// Unpack the single vector index into indices
465  template<class size_vec_t>
466  void unpack_indices(size_t ix, size_vec_t &index) {
467  if (ix>=total_size()) {
468  O2SCL_ERR((((std::string)"Value of index ")+szttos(ix)+
469  " greater than or equal to total size"+
470  szttos(total_size())+
471  " in tensor::unpack_indices().").c_str(),
472  exc_eindex);
473  return;
474  }
475  size_t ix2, sub_size;
476  for(size_t i=0;i<rk;i++) {
477  if (i==rk-1) {
478  index[i]=ix;
479  } else {
480  sub_size=1;
481  for(size_t j=i+1;j<rk;j++) sub_size*=size[j];
482  index[i]=ix/sub_size;
483  // (Remember we're doing integer arithmetic here.)
484  ix-=sub_size*(ix/sub_size);
485  }
486  }
487  return;
488  }
489  //@}
490 
491  /// \name Minimum and maximum
492  //@{
493  /** \brief Compute the minimum value in the tensor
494  */
495  double min_value() {
496  return o2scl::vector_min_value<vec_t,double>(total_size(),data);
497  }
498 
499  /** \brief Compute the index of the minimum value in the tensor
500  */
501  void min_index(vec_size_t &index) {
502  size_t ix=o2scl::vector_min_index<vec_t,double>(total_size(),data);
503  unpack_indices(ix,index);
504  return;
505  }
506 
507  /** \brief Compute the index of the minimum value in the tensor
508  and return the minimum
509  */
510  void min(vec_size_t &index, double &val) {
511  size_t ix;
512  o2scl::vector_min<vec_t,double>(total_size(),data,ix,val);
513  unpack_indices(ix,index);
514  return;
515  }
516 
517  /** \brief Compute the maximum value in the tensor
518  */
519  double max_value() {
520  return o2scl::vector_max_value<vec_t,double>(total_size(),data);
521  }
522 
523  /** \brief Compute the index of the maximum value in the tensor
524  */
525  void max_index(vec_size_t &index) {
526  size_t ix=o2scl::vector_max_index<vec_t,double>(total_size(),data);
527  unpack_indices(ix,index);
528  return;
529  }
530 
531  /** \brief Compute the index and value of the maximum value in the tensor
532  and return the maximum
533  */
534  void max(vec_size_t &index, double &val) {
535  size_t ix;
536  o2scl::vector_max<vec_t,double>(total_size(),data,ix,val);
537  unpack_indices(ix,index);
538  return;
539  }
540 
541  /** \brief Compute the minimum and maximum values in the tensor
542  */
543  void minmax_value(double &min, double &max) {
544  return o2scl::vector_minmax_value<vec_t,double>(total_size(),data,min,max);
545  }
546 
547  /** \brief Compute the indices of the minimum and maximum values in the tensor
548  */
549  void minmax_index(vec_size_t &index_min, vec_size_t &index_max) {
550  size_t ix_min, ix_max;
551  o2scl::vector_minmax_index<vec_t,double>
552  (total_size(),data,ix_min,ix_max);
553  unpack_indices(ix_min,index_min);
554  unpack_indices(ix_max,index_max);
555  return;
556  }
557 
558  /** \brief Compute the indices and values of the maximum and minimum
559  in the tensor
560  */
561  void minmax(vec_size_t &index, size_t &index_min, double &min,
562  size_t &index_max, double &max) {
563  size_t ix_min, ix_max;
564  o2scl::vector_minmax<vec_t,double>(total_size(),data,ix_min,min,
565  ix_max,max);
566  unpack_indices(ix_min,index_min);
567  unpack_indices(ix_max,index_max);
568  return;
569  }
570  //@}
571 
572  };
573 
574  /** \brief Rank 1 tensor
575  */
576  template<class vec_t=std::vector<double>,
577  class vec_size_t=std::vector<size_t> > class tensor1 :
578  public tensor<vec_t,vec_size_t> {
579 
580  public:
581 
582  /// Create an empty tensor
583  tensor1() : tensor<vec_t,vec_size_t>() {}
584 
585  /// Create a rank 1 tensory of size \c sz
586  tensor1(size_t sz) : tensor<vec_t,vec_size_t>() {
587  vec_size_t sizex(1);
588  sizex[0]=sz;
589  this->resize(1,sizex);
590  }
591 
592  /// Get the element indexed by \c ix
593  double &get(size_t ix) {
594  return tensor<vec_t,vec_size_t>::get(&ix);
595  }
596 
597  /// Get the element indexed by \c ix
598  const double &get(size_t ix) const {
599  return tensor<vec_t,vec_size_t>::get(&ix);
600  }
601 
602  /// Set the element indexed by \c index to value \c val
603  void set(size_t index, double val)
604  { tensor<vec_t,vec_size_t>::set(&index,val); }
605 
606  /** \brief Set the element indexed by \c index to value \c val
607 
608  (We have to explicitly provide this version since the set()
609  function is overloaded in this child of \ref tensor.)
610  */
611  template<class size_vec_t>
612  void set(const size_vec_t &index, double val) {
614  }
615 
616  /// Get an element using array-like indexing
617  double &operator[](size_t ix) { return this->data[ix]; }
618 
619  /// Get an element using array-like indexing (const version)
620  const double &operator[](size_t ix) const { return this->data[ix]; }
621 
622  /// Get an element using operator()
623  double &operator()(size_t ix) { return this->data[ix]; }
624 
625  /// Get an element using operator() (const version)
626  const double &operator()(size_t ix) const { return this->data[ix]; }
627  };
628 
629  /** \brief Rank 2 tensor
630  */
631  template<class vec_t=std::vector<double>,
632  class vec_size_t=std::vector<size_t> > class tensor2 :
633  public tensor<vec_t,vec_size_t> {
634 
635  public:
636 
637  /// Create an empty tensor
638  tensor2() : tensor<vec_t,vec_size_t>() {}
639 
640  /// Create a rank 2 tensor of size \c (sz,sz2)
641  tensor2(size_t sz, size_t sz2) : tensor<vec_t,vec_size_t>() {
642  this->rk=2;
643  this->size.resize(2);
644  this->size[0]=sz;
645  this->size[1]=sz2;
646  size_t tot=sz*sz2;
647  this->data.resize(tot);
648  }
649 
650  /// Get the element indexed by \c (ix1,ix2)
651  double &get(size_t ix1, size_t ix2) {
652  size_t sz[2]={ix1,ix2};
653  return tensor<vec_t,vec_size_t>::get(sz);
654  }
655 
656  /// Get the element indexed by \c (ix1,ix2)
657  const double &get(size_t ix1, size_t ix2) const {
658  size_t sz[2]={ix1,ix2};
659  return tensor<vec_t,vec_size_t>::get(sz);
660  }
661 
662  /// Set the element indexed by \c (ix1,ix2) to value \c val
663  void set(size_t ix1, size_t ix2, double val) {
664  size_t sz[2]={ix1,ix2};
666  return;
667  }
668 
669  /** \brief Set the element indexed by \c index to value \c val
670 
671  (We have to explicitly provide this version since the set()
672  function is overloaded in this child of \ref tensor.)
673  */
674  template<class size_vec_t>
675  void set(const size_vec_t &index, double val) {
677  return;
678  }
679 
680  /// Get the element indexed by \c (ix1,ix2)
681  double &operator()(size_t ix, size_t iy)
682  { return this->data[ix*this->size[1]+iy]; }
683 
684  /// Get the element indexed by \c (ix1,ix2) (const version)
685  const double &operator()(size_t ix, size_t iy) const
686  { return this->data[ix*this->size[1]+iy]; }
687  };
688 
689  /** \brief Rank 3 tensor
690  */
691  template<class vec_t=std::vector<double>,
692  class vec_size_t=std::vector<size_t> > class tensor3 :
693  public tensor<vec_t,vec_size_t> {
694 
695  public:
696 
697  /// Create an empty tensor
698  tensor3() : tensor<vec_t,vec_size_t>() {}
699 
700  /// Create a rank 3 tensor of size \c (sz,sz2,sz3)
701  tensor3(size_t sz, size_t sz2, size_t sz3) :
702  tensor<vec_t,vec_size_t>() {
703  this->rk=3;
704  this->size.resize(3);
705  this->size[0]=sz;
706  this->size[1]=sz2;
707  this->size[2]=sz3;
708  size_t tot=sz*sz2*sz3;
709  this->data.resize(tot);
710  }
711 
712  /// Get the element indexed by \c (ix1,ix2,ix3)
713  double &get(size_t ix1, size_t ix2, size_t ix3) {
714  size_t sz[3]={ix1,ix2,ix3};
715  return tensor<vec_t,vec_size_t>::get(sz);
716  }
717 
718  /// Get the element indexed by \c (ix1,ix2,ix3)
719  const double &get(size_t ix1, size_t ix2, size_t ix3) const {
720  size_t sz[3]={ix1,ix2,ix3};
721  return tensor<vec_t,vec_size_t>::get(sz);
722  }
723 
724  /// Set the element indexed by \c (ix1,ix2,ix3) to value \c val
725  void set(size_t ix1, size_t ix2, size_t ix3, double val) {
726  size_t sz[3]={ix1,ix2,ix3};
728  return;
729  }
730 
731  /** \brief Set the element indexed by \c index to value \c val
732 
733  (We have to explicitly provide this version since the set()
734  function is overloaded in this child of \ref tensor.)
735  */
736  template<class size_vec_t>
737  void set(const size_vec_t &index, double val) {
739  return;
740  }
741 
742  };
743 
744  /** \brief Rank 4 tensor
745  */
746  template<class vec_t=std::vector<double>,
747  class vec_size_t=std::vector<size_t> > class tensor4 :
748  public tensor<vec_t,vec_size_t> {
749 
750  public:
751 
752  /// Create an empty tensor
753  tensor4() : tensor<vec_t,vec_size_t>() {}
754 
755  /// Create a rank 4 tensor of size \c (sz,sz2,sz3,sz4)
756  tensor4(size_t sz, size_t sz2, size_t sz3, size_t sz4) :
757  tensor<vec_t,vec_size_t>() {
758  this->rk=4;
759  this->size.resize(4);
760  this->size[0]=sz;
761  this->size[1]=sz2;
762  this->size[2]=sz3;
763  this->size[3]=sz4;
764  size_t tot=sz*sz2*sz3*sz4;
765  this->data.resize(tot);
766  }
767 
768  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
769  double &get(size_t ix1, size_t ix2, size_t ix3, size_t ix4) {
770  size_t sz[4]={ix1,ix2,ix3,ix4};
771  return tensor<vec_t,vec_size_t>::get(sz);
772  }
773 
774  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
775  const double &get(size_t ix1, size_t ix2, size_t ix3,
776  size_t ix4) const {
777  size_t sz[4]={ix1,ix2,ix3,ix4};
778  return tensor<vec_t,vec_size_t>::get(sz);
779  }
780 
781  /// Set the element indexed by \c (ix1,ix2,ix3,ix4) to value \c val
782  void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4,
783  double val) {
784  size_t sz[4]={ix1,ix2,ix3,ix4};
786  return;
787  }
788 
789  /** \brief Set the element indexed by \c index to value \c val
790 
791  (We have to explicitly provide this version since the set()
792  function is overloaded in this child of \ref tensor.)
793  */
794  template<class size_vec_t>
795  void set(const size_vec_t &index, double val) {
797  return;
798  }
799  };
800 
801 #ifndef DOXYGEN_NO_O2NS
802 }
803 #endif
804 
805 #endif
806 
807 
808 
tensor(size_t rank, const size_vec_t &dim)
Create a tensor of rank rank with sizes given in dim.
Definition: tensor.h:154
size_t pack_indices(const size_vec_t &index)
Pack the indices into a single vector index.
Definition: tensor.h:439
double max_value()
Compute the maximum value in the tensor.
Definition: tensor.h:519
Rank 1 tensor.
Definition: tensor.h:577
const double & operator()(size_t ix, size_t iy) const
Get the element indexed by (ix1,ix2) (const version)
Definition: tensor.h:685
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
Rank 3 tensor.
Definition: tensor.h:692
void max(vec_size_t &index, double &val)
Compute the index and value of the maximum value in the tensor and return the maximum.
Definition: tensor.h:534
void min(vec_size_t &index, double &val)
Compute the index of the minimum value in the tensor and return the minimum.
Definition: tensor.h:510
const double & operator[](size_t ix) const
Get an element using array-like indexing (const version)
Definition: tensor.h:620
void set(const size_vec_t &index, double val)
Set the element indexed by index to value val.
Definition: tensor.h:196
invalid argument supplied by user
Definition: err_hnd.h:59
double & operator()(size_t ix, size_t iy)
Get the element indexed by (ix1,ix2)
Definition: tensor.h:681
void unpack_indices(size_t ix, size_vec_t &index)
Unpack the single vector index into indices.
Definition: tensor.h:466
void max_index(vec_size_t &index)
Compute the index of the maximum value in the tensor.
Definition: tensor.h:525
void resize(size_t rank, const size_vec_t &dim)
Resize the tensor to rank rank with sizes given in dim.
Definition: tensor.h:374
void minmax(vec_size_t &index, size_t &index_min, double &min, size_t &index_max, double &max)
Compute the indices and values of the maximum and minimum in the tensor.
Definition: tensor.h:561
vec_t data
The data.
Definition: tensor.h:129
size_t get_size(size_t i) const
Returns the size of the ith index.
Definition: tensor.h:404
void min_index(vec_size_t &index)
Compute the index of the minimum value in the tensor.
Definition: tensor.h:501
size_t get_rank() const
Return the rank of the tensor.
Definition: tensor.h:401
tensor2()
Create an empty tensor.
Definition: tensor.h:638
vec_size_t size
A rank-sized vector of the sizes of each dimension.
Definition: tensor.h:132
const vec_t & get_data() const
Return the full data vector.
Definition: tensor.h:420
tensor3()
Create an empty tensor.
Definition: tensor.h:698
void minmax_value(double &min, double &max)
Compute the minimum and maximum values in the tensor.
Definition: tensor.h:543
size_t rk
Rank.
Definition: tensor.h:135
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void set_all(double x)
Set all elements in a tensor to some fixed value.
Definition: tensor.h:228
tensor1()
Create an empty tensor.
Definition: tensor.h:583
void minmax_index(vec_size_t &index_min, vec_size_t &index_max)
Compute the indices of the minimum and maximum values in the tensor.
Definition: tensor.h:549
tensor1(size_t sz)
Create a rank 1 tensory of size sz.
Definition: tensor.h:586
tensor4()
Create an empty tensor.
Definition: tensor.h:753
ubvector_slice vector_slice(size_t ix, const size_vec_t &index)
Fix all but one index to create a vector.
Definition: tensor.h:329
size_t total_size() const
Returns the size of the tensor (the product of the sizes over every index)
Definition: tensor.h:427
tensor()
Create an empty tensor with zero rank.
Definition: tensor.h:142
const vec_size_t & get_size_arr() const
Return the full vector of sizes.
Definition: tensor.h:415
double & get(const size_vec_t &index)
Get the element indexed by index.
Definition: tensor.h:237
Rank 2 tensor.
Definition: tensor.h:632
Rank 4 tensor.
Definition: tensor.h:747
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor.h:184
Invalid index for array or matrix.
Definition: err_hnd.h:123
tensor3(size_t sz, size_t sz2, size_t sz3)
Create a rank 3 tensor of size (sz,sz2,sz3)
Definition: tensor.h:701
double & operator[](size_t ix)
Get an element using array-like indexing.
Definition: tensor.h:617
std::string szttos(size_t x)
Convert a size_t to a string.
tensor2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
Definition: tensor.h:641
const double & operator()(size_t ix) const
Get an element using operator() (const version)
Definition: tensor.h:626
tensor4(size_t sz, size_t sz2, size_t sz3, size_t sz4)
Create a rank 4 tensor of size (sz,sz2,sz3,sz4)
Definition: tensor.h:756
double & operator()(size_t ix)
Get an element using operator()
Definition: tensor.h:623
Tensor class with arbitrary dimensions.
Definition: tensor.h:120
double min_value()
Compute the minimum value in the tensor.
Definition: tensor.h:495

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).