tensor.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 #include <o2scl/table3d.h>
46 #include <o2scl/misc.h>
47 
48 #ifndef DOXYGEN_NO_O2NS
49 namespace o2scl {
50 #endif
51 
52  /** \brief Index specification
53  */
54  class index_spec {
55 
56  public:
57 
58  /// Type of specification
59  size_t type;
60  /// First argument
61  size_t ix1;
62  /// Second argument
63  size_t ix2;
64  /// Third argument
65  size_t ix3;
66  /// First double argument
67  double val1;
68  /// Second double argument
69  double val2;
70  /// Third double argument
71  double val3;
72 
73  /// \name Possible values for type
74  //@{
75  /// Empty specification
76  static const size_t empty=0;
77  /// Retain an index
78  static const size_t index=1;
79  /// Fix the value of an index
80  static const size_t fixed=2;
81  /// Sum over an index
82  static const size_t sum=3;
83  /// Perform a trace (sum over two indices)
84  static const size_t trace=4;
85  /// Reverse an index
86  static const size_t reverse=5;
87  /// Choose a new range for an index
88  static const size_t range=6;
89  /// Interpolate a value to fix an index
90  static const size_t interp=7;
91  /// Interpolate a value to set a new grid (fixed bin number)
92  static const size_t grid=8;
93  /// Interpolate a value to set a new grid (fixed bin width)
94  static const size_t gridw=9;
95  //@}
96 
97  /// Default constructor
99  type=0;
100  ix1=0;
101  ix2=0;
102  ix3=0;
103  val1=0.0;
104  val2=0.0;
105  val3=0.0;
106  }
107 
108  /// Explicit full constructor
109  index_spec(size_t typ, size_t i1, size_t i2=0, size_t i3=0,
110  double v1=0.0, double v2=0.0, double v3=0.0) {
111  type=typ;
112  ix1=i1;
113  ix2=i2;
114  ix3=i3;
115  val1=v1;
116  val2=v2;
117  val3=v3;
118  }
119 
120  };
121 
122  /// \name Tensor index functions in src/base/tensor.h
123  //@{
124  /** \brief Choose an index
125  */
126  index_spec ix_index(size_t ix);
127 
128  /** \brief Fix index \c ix to value \c ix2
129  */
130  index_spec ix_fixed(size_t ix, size_t ix2);
131 
132  /** \brief Sum over index \c ix
133  */
134  index_spec ix_sum(size_t ix);
135 
136  /** \brief Perform a trace over indices \c ix and \c ix2
137  */
138  index_spec ix_trace(size_t ix, size_t ix2);
139 
140  /** \brief Reverse index \c ix
141  */
142  index_spec ix_reverse(size_t ix);
143 
144  /** \brief Index covers a range of values
145  */
146  index_spec ix_range(size_t ix, size_t start, size_t end);
147 
148  /** \brief Interpolate value \c v into index \c ix
149  (for \ref o2scl::tensor_grid only)
150  */
151  index_spec ix_interp(size_t ix, double v);
152 
153  /** \brief Interpolate grid with fixed number of bins into index \c ix
154  (for \ref o2scl::tensor_grid only)
155  */
156  index_spec ix_grid(size_t ix, double begin, double end, size_t n_bins,
157  bool log=false);
158 
159  /** \brief Interpolate grid with fixed bin width into index \c ix
160  (for \ref o2scl::tensor_grid only)
161  */
162  index_spec ix_gridw(size_t ix, double begin, double end, double width,
163  bool log=false);
164  //@}
165 
166  /** \brief Tensor class with arbitrary dimensions
167 
168  The elements of a tensor are typically specified as a list of
169  <tt>size_t</tt> numbers with length equal to the tensor rank.
170  For a rank-4 tensor named \c t, the element
171  <tt>t[1][2][0][3]</tt> can be obtained with something similar to
172  \code
173  size_t ix[4]={1,2,0,3};
174  double x=t.get(ix);
175  \endcode
176 
177  Empty tensors have zero rank.
178 
179  The type <tt>vec_t</tt> can be any vector type with
180  <tt>operator[]</tt>, <tt>size()</tt> and <tt>resize()</tt>
181  methods. The type <tt>vec_size_t</tt> can be any integer-like
182  vector type with <tt>operator[]</tt>, <tt>size()</tt> and
183  <tt>resize()</tt> methods.
184 
185  For I/O with tensors, see \ref o2scl_hdf::hdf_file::setd_ten()
186  and \ref o2scl_hdf::hdf_file::getd_ten() . See also
187  the discussion in the sections \ref tensor_subsect and
188  \ref vec_io_cont_subsect of the user's guide.
189 
190  The storage pattern is a generalization of row-major order.
191  In the case of a 4-rank tensor, the location of a generic
192  element is
193  \f[
194  \left( \left( i_0 s_1 + i_1 \right) s_2 + i_2 \right) s_3 + i_3 \, .
195  \f]
196  In this case the distance between two elements \f$(i_0,i_1,
197  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
198  \f$, the distance between two elements \f$(i_0,i_1,i_2, i_3)\f$
199  and \f$ (i_0,i_1+1,i_2,i_3) \f$ is \f$ s_2 s_3 \f$, and the
200  elements \f$(i_0,i_1,i_2,i_3)\f$ and \f$ (i_0,i_1,i_2,i_3+1) \f$
201  are adjacent.
202 
203  \note Slices of tensors are subsets obtained from fixing the
204  index of several dimensions, while letting other dimensions
205  vary. For a slice with one dimension not fixed, see \ref
206  vector_slice(). The \ref o2scl::tensor::vector_slice() function
207  should clearly work for uBlas vectors, and seems to work with
208  std::vector objects also, but latter use has not been fully
209  tested.
210 
211  \future Create an operator[] for tensor and not just tensor1?
212 
213  \future Could implement arithmetic operators + and - and some
214  different products.
215 
216  \future Implement copies to and from vector
217  and matrices
218 
219  \future Implement tensor contractions, i.e. tensor
220  = tensor * tensor
221 
222  \future Could be interesting to write an iterator for this class.
223 
224  */
225  template<class data_t=double, class vec_t=std::vector<data_t>,
226  class vec_size_t=std::vector<size_t> > class tensor {
227 
228  public:
229 
230 #ifndef DOXYGEN_INTERNAL
231 
232  protected:
233 
234  /// The data
235  vec_t data;
236 
237  /// A rank-sized vector of the sizes of each dimension
238  vec_size_t size;
239 
240  /// Rank
241  size_t rk;
242 
243 #endif
244 
245  public:
246 
247  /// Create an empty tensor with zero rank
248  tensor() {
249  rk=0;
250  }
251 
252  /** \brief Create a tensor of rank \c rank with sizes given in \c dim
253 
254  The parameter \c dim must be a pointer to a vector of sizes with
255  length \c rank. If the user requests any of the sizes to be
256  zero, this constructor will call the error handler, create an
257  empty tensor, and will allocate no memory.
258  */
259  template<class size_vec_t>
260  tensor(size_t rank, const size_vec_t &dim) {
261  if (rank==0) {
262  rk=0;
263  } else {
264  rk=rank;
265  for(size_t i=0;i<rk;i++) {
266  if (dim[i]==0) {
267  rk=0;
268  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
269  "rank for index "+szttos(i)+
270  " in tensor::tensor(size_t,size_vec_t)").c_str(),
271  exc_einval);
272  }
273  }
274  size.resize(rk);
275  size_t tot=1;
276  for(size_t i=0;i<rk;i++) {
277  size[i]=dim[i];
278  tot*=size[i];
279  }
280  data.resize(tot);
281  }
282  }
283 
284  ~tensor() {
285  }
286 
287  /// \name Method to check for valid object
288  //@{
289  /** \brief Check that the \ref o2scl::tensor object is valid
290  */
291  void is_valid() const {
292  if (rk==0) {
293  if (data.size()!=0) {
294  O2SCL_ERR2("Rank is zero but the data vector has non-zero size ",
295  "in tensor::is_valid().",o2scl::exc_esanity);
296  }
297  }
298 
299  if (rk!=size.size()) {
300  O2SCL_ERR2("Rank does not match size vector size ",
301  "in tensor::is_valid().",o2scl::exc_esanity);
302  }
303 
304  if (rk>0) {
305  size_t tot=1;
306  for(size_t i=0;i<rk;i++) tot*=size[i];
307  if (tot==0) {
308  O2SCL_ERR2("One entry in the size vector is zero ",
309  "in tensor::is_valid().",o2scl::exc_esanity);
310  }
311  if (tot!=data.size()) {
312  O2SCL_ERR2("Product of size vector entries does not match data ",
313  "vector size in tensor::is_valid().",o2scl::exc_esanity);
314  }
315  }
316 
317  return;
318  }
319  //@}
320 
321  /// \name Copy constructors
322  //@{
323  /** \brief Copy using <tt>operator()</tt>
324  */
327  rk=t.rk;
328  data=t.data;
329  size=t.size;
330  }
331 
332  /** \brief Copy using <tt>operator=()</tt>
333  */
336  if (this!=&t) {
337  rk=t.rk;
338  data=t.data;
339  size=t.size;
340  }
341  return *this;
342  }
343  //@}
344 
345  /// \name Clear method
346  //@{
347  /// Clear the tensor of all data and free allocated memory
348  void clear() {
349  rk=0;
350  data.resize(0);
351  size.resize(0);
352  return;
353  }
354  //@}
355 
356  /// \name Set functions
357  //@{
358  /// Set the element indexed by \c index to value \c val
359  template<class size_vec_t>
360  void set(const size_vec_t &index, data_t val) {
361 #if O2SCL_NO_RANGE_CHECK
362 #else
363  if (rk==0) {
364  O2SCL_ERR("Empty tensor in tensor::set().",exc_einval);
365  }
366  if (index[0]>=size[0]) {
367  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
368  " greater than or equal to size[0]="+
369  szttos(size[0])+" in tensor::set().").c_str(),
370  exc_eindex);
371  }
372 #endif
373  size_t ix=index[0];
374  for(size_t i=1;i<rk;i++) {
375 #if O2SCL_NO_RANGE_CHECK
376 #else
377  if (index[i]>=size[i]) {
378  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
379  szttos(index[i])+" greater than or equal to size "+
380  szttos(size[i])+" in tensor::set().").c_str(),
381  exc_eindex);
382  }
383 #endif
384  ix*=size[i];
385  ix+=index[i];
386  }
387  data[ix]=val;
388  return;
389  }
390 
391  /// Set all elements in a tensor to some fixed value
392  void set_all(data_t x) {
393  for(size_t i=0;i<total_size();i++) data[i]=x;
394  return;
395  }
396 
397  /** \brief Swap the data vector
398  */
399  void swap_data(vec_t &dat) {
400  if (data.size()!=dat.size()) {
401  O2SCL_ERR2("Size of new vector does not equal tensor size in ",
402  "tensor::swap_data().",o2scl::exc_einval);
403  }
404  std::swap(dat,data);
405  return;
406  }
407  //@}
408 
409  /// \name Get functions
410  //@{
411  /// Get the element indexed by \c index
412  template<class size_vec_t> data_t &get(const size_vec_t &index) {
413 #if O2SCL_NO_RANGE_CHECK
414 #else
415  if (rk==0) {
416  O2SCL_ERR("Empty tensor in tensor::get().",exc_einval);
417  }
418  if (index[0]>=size[0]) {
419  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
420  " greater than or equal to size[0]="+
421  szttos(size[0])+" in tensor::get().").c_str(),
422  exc_eindex);
423  }
424 #endif
425  size_t ix=index[0];
426  for(size_t i=1;i<rk;i++) {
427 #if O2SCL_NO_RANGE_CHECK
428 #else
429  if (index[i]>=size[i]) {
430  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
431  szttos(index[i])+" greater than or equal to size "+
432  szttos(size[i])+" in tensor::get().").c_str(),
433  exc_eindex);
434  }
435 #endif
436  ix*=size[i];
437  ix+=index[i];
438  }
439  return data[ix];
440  }
441 
442  /// Get a const reference to the element indexed by \c index
443  template<class size_vec_t>
444  data_t const &get(const size_vec_t &index) const {
445 #if O2SCL_NO_RANGE_CHECK
446 #else
447  if (rk==0) {
448  O2SCL_ERR("Empty tensor in tensor::get() (const).",exc_einval);
449  }
450  if (index[0]>=size[0]) {
451  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
452  " greater than or equal to size[0]="+
453  szttos(size[0])+" in tensor::get() (const).").c_str(),
454  exc_eindex);
455  }
456 #endif
457  size_t ix=index[0];
458  for(size_t i=1;i<rk;i++) {
459 #if O2SCL_NO_RANGE_CHECK
460 #else
461  if (index[i]>=size[i]) {
462  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
463  szttos(index[i])+" greater than or equal to size "+
464  szttos(size[i])+" in tensor::get() (const).").c_str(),
465  exc_eindex);
466  }
467 #endif
468  ix*=size[i];
469  ix+=index[i];
470  }
471  return data[ix];
472  }
473  //@}
474 
475  typedef boost::numeric::ublas::vector_slice<
476  boost::numeric::ublas::vector<data_t> > ubvector_slice;
477  typedef boost::numeric::ublas::slice slice;
478 
479  /// \name Slice function
480  //@{
481  /** \brief Fix all but one index to create a vector
482 
483  This fixes all of the indices to the values given in \c index
484  except for the index number \c ix, and returns the
485  corresponding vector, whose length is equal to the size of the
486  tensor in that index. The value <tt>index[ix]</tt> is ignored.
487 
488  For example, for a rank 3 tensor allocated with
489  \code
490  tensor t;
491  size_t dim[3]={3,4,5};
492  t.resize(3,dim);
493  \endcode
494  the following code
495  \code
496  size_t index[3]={1,0,3};
497  ubvector_view v=t.vector_slice(1,index);
498  \endcode
499  Gives a vector \c v of length 4 which refers to the values
500  <tt>t(1,0,3)</tt>, <tt>t(1,1,3)</tt>, <tt>t(1,2,3)</tt>, and
501  <tt>t(1,3,3)</tt>.
502  */
503  template<class size_vec_t>
504  ubvector_slice vector_slice(size_t ix, const size_vec_t &index) {
505  if (ix+1>rk) {
506  O2SCL_ERR((((std::string)"Specified index ")+szttos(ix)+
507  " greater than or equal to rank "+szttos(rk)+
508  " in tensor::vector_slice()").c_str(),
509  exc_eindex);
510  }
511  size_t start;
512  if (ix==0) start=0;
513  else start=index[0];
514  for(size_t i=1;i<rk;i++) {
515  start*=size[i];
516  if (i!=ix) start+=index[i];
517  }
518  size_t stride=1;
519  for(size_t i=ix+1;i<rk;i++) stride*=size[i];
520  return ubvector_slice(data,slice(start,stride,size[ix]));
521  }
522  //@}
523 
524 #ifdef O2SCL_NEVER_DEFINED
525 
526  /** \brief Fix all but two indices to create a matrix
527 
528  One could use ublas::make_matrix_from_pointer() and
529  then ublas matrix slicing to perform this operation, but
530  only if the template types are built on ublas objects.
531  */
532  template<class size_vec_t> ubmatrix_array matrix_slice() {
533  }
534 
535 #endif
536 
537  /// \name Resize method
538  //@{
539  /** \brief Resize the tensor to rank \c rank with sizes
540  given in \c dim
541 
542  The parameter \c dim must be a vector of sizes with a length
543  equal to \c rank. This resize method is always destructive.
544 
545  If the user requests any of the sizes to be zero, this
546  function will call the error handler.
547  */
548  template<class size_vec_t>
549  void resize(size_t rank, const size_vec_t &dim) {
550  for(size_t i=0;i<rank;i++) {
551  if (dim[i]==0) {
552  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
553  "rank for index "+szttos(i)+" in tensor::"+
554  "resize().").c_str(),exc_einval);
555  }
556  }
557  rk=rank;
558  size.resize(rk);
559  if (rk==0) {
560  data.resize(0);
561  } else {
562  size_t tot=1;
563  for(size_t i=0;i<rk;i++) {
564  size[i]=dim[i];
565  tot*=size[i];
566  }
567  data.resize(tot);
568  }
569  return;
570  }
571  //@}
572 
573  /// \name Size functions
574  //@{
575  /// Return the rank of the tensor
576  size_t get_rank() const { return rk; }
577 
578  /// Returns the size of the ith index
579  size_t get_size(size_t i) const {
580  if (i>=rk) {
581  O2SCL_ERR((((std::string)"Specified index ")+szttos(i)+
582  " greater than or equal to rank "+szttos(rk)+
583  " in tensor::get_size()").c_str(),
584  exc_einval);
585  }
586  return size[i];
587  }
588 
589  /// Return the full vector of sizes
590  const vec_size_t &get_size_arr() const {
591  return size;
592  }
593 
594  /// Return the full data vector
595  const vec_t &get_data() const {
596  return data;
597  }
598 
599  /** \brief Returns the size of the tensor (the product of
600  the sizes over every index)
601  */
602  size_t total_size() const {
603  if (rk==0) return 0;
604  size_t tot=1;
605  for(size_t i=0;i<rk;i++) tot*=size[i];
606  return tot;
607  }
608  //@}
609 
610  /// \name Index manipulation
611  //@{
612  /// Pack the indices into a single vector index
613  template<class size_vec_t>
614  size_t pack_indices(const size_vec_t &index) {
615  if (rk==0) {
616  O2SCL_ERR("Empty tensor in tensor::pack_indices().",exc_einval);
617  return 0;
618  }
619  if (index[0]>=size[0]) {
620  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
621  " greater than or equal to size[0]="+
622  szttos(size[0])+" in tensor::pack_indices().").c_str(),
623  exc_eindex);
624  }
625  size_t ix=index[0];
626  for(size_t i=1;i<rk;i++) {
627  if (index[i]>=size[i]) {
628  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
629  szttos(index[i])+" greater than or equal to size "+
630  szttos(size[i])+" in tensor::pack_indices().").c_str(),
631  exc_eindex);
632  }
633  ix*=size[i];
634  ix+=index[i];
635  }
636  return ix;
637  }
638 
639  /// Unpack the single vector index into indices
640  template<class size_vec_t>
641  void unpack_index(size_t ix, size_vec_t &index) {
642  if (ix>=total_size()) {
643  O2SCL_ERR((((std::string)"Value of index ")+szttos(ix)+
644  " greater than or equal to total size"+
645  szttos(total_size())+
646  " in tensor::unpack_index().").c_str(),
647  exc_eindex);
648  return;
649  }
650  size_t ix2, sub_size;
651  for(size_t i=0;i<rk;i++) {
652  if (i==rk-1) {
653  index[i]=ix;
654  } else {
655  sub_size=1;
656  for(size_t j=i+1;j<rk;j++) sub_size*=size[j];
657  index[i]=ix/sub_size;
658  // (Remember we're doing integer arithmetic here.)
659  ix-=sub_size*(ix/sub_size);
660  }
661  }
662  return;
663  }
664  //@}
665 
666  /// \name Minimum, maximum, and sum
667  //@{
668  /** \brief Compute the minimum value in the tensor
669  */
670  data_t min_value() {
671  return o2scl::vector_min_value<vec_t,data_t>(total_size(),data);
672  }
673 
674  /** \brief Compute the index of the minimum value in the tensor
675  */
676  void min_index(vec_size_t &index) {
677  size_t ix=o2scl::vector_min_index<vec_t,data_t>(total_size(),data);
678  unpack_index(ix,index);
679  return;
680  }
681 
682  /** \brief Compute the index of the minimum value in the tensor
683  and return the minimum
684  */
685  void min(vec_size_t &index, data_t &val) {
686  size_t ix;
687  o2scl::vector_min<vec_t,data_t>(total_size(),data,ix,val);
688  unpack_index(ix,index);
689  return;
690  }
691 
692  /** \brief Compute the maximum value in the tensor
693  */
694  data_t max_value() {
695  return o2scl::vector_max_value<vec_t,data_t>(total_size(),data);
696  }
697 
698  /** \brief Compute the index of the maximum value in the tensor
699  */
700  void max_index(vec_size_t &index) {
701  size_t ix=o2scl::vector_max_index<vec_t,data_t>(total_size(),data);
702  unpack_index(ix,index);
703  return;
704  }
705 
706  /** \brief Compute the index and value of the maximum value in the tensor
707  and return the maximum
708  */
709  void max(vec_size_t &index, data_t &val) {
710  size_t ix;
711  o2scl::vector_max<vec_t,data_t>(total_size(),data,ix,val);
712  unpack_index(ix,index);
713  return;
714  }
715 
716  /** \brief Compute the minimum and maximum values in the tensor
717  */
718  void minmax_value(data_t &min, data_t &max) {
719  return o2scl::vector_minmax_value<vec_t,data_t>(total_size(),data,min,max);
720  }
721 
722  /** \brief Compute the indices of the minimum and maximum values in the tensor
723  */
724  void minmax_index(vec_size_t &index_min, vec_size_t &index_max) {
725  size_t ix_min, ix_max;
726  o2scl::vector_minmax_index<vec_t,data_t>
727  (total_size(),data,ix_min,ix_max);
728  unpack_index(ix_min,index_min);
729  unpack_index(ix_max,index_max);
730  return;
731  }
732 
733  /** \brief Compute the indices and values of the maximum and minimum
734  in the tensor
735  */
736  void minmax(vec_size_t &index, size_t &index_min, data_t &min,
737  size_t &index_max, data_t &max) {
738  size_t ix_min, ix_max;
739  o2scl::vector_minmax<vec_t,data_t>(total_size(),data,ix_min,min,
740  ix_max,max);
741  unpack_index(ix_min,index_min);
742  unpack_index(ix_max,index_max);
743  return;
744  }
745 
746  /** \brief Return the sum over every element in the tensor
747  */
748  double total_sum() const {
749  if (rk==0) return 0.0;
750  double tot=0.0;
751  for(size_t i=0;i<data.size();i++) {
752  tot+=data[i];
753  }
754  return tot;
755  }
756  //@}
757 
758  /// \name Slicing and converting to table3d objects
759  //@{
760  /** \brief Convert to a \ref o2scl::table3d object by
761  summing over all but two indices
762  */
764  (size_t ix_x, size_t ix_y, table3d &tab, std::string x_name="x",
765  std::string y_name="y", std::string slice_name="z") {
766 
767  // Get current table3d grid
768  size_t nx, ny;
769  tab.get_size(nx,ny);
770 
771  if (nx==0 && ny==0) {
772 
773  if (x_name.length()==0) x_name="x";
774  if (y_name.length()==0) y_name="y";
775 
776  // If there's no grid, then create a grid in the table3d
777  // object which just enumerates the indices
778  std::vector<double> grid_x(size[ix_x]), grid_y(size[ix_y]);
779  for(size_t i=0;i<size[ix_x];i++) {
780  grid_x[i]=((double)i);
781  }
782  for(size_t i=0;i<size[ix_y];i++) {
783  grid_y[i]=((double)i);
784  }
785  tab.set_xy("x",grid_x.size(),grid_x,
786  "y",grid_y.size(),grid_y);
787  // Now that the grid is set, get nx and ny
788  tab.get_size(nx,ny);
789  }
790 
791  // Check that the grids are commensurate
792  if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
793  O2SCL_ERR2("Grids not commensurate in ",
794  "tensor_grid::convert_table3d_sum().",exc_einval);
795  }
796 
797  tab.set_slice_all(slice_name,0.0);
798 
799  std::vector<size_t> ix;
800  for(size_t i=0;i<this->total_size();i++) {
801  this->unpack_index(i,ix);
802  tab.set(ix[ix_x],ix[ix_y],slice_name,
803  tab.get(ix[ix_x],ix[ix_y],slice_name)+
804  this->data[i]);
805  }
806 
807  return;
808  }
809  //@}
810 
811  /** \brief Rearrange, sum and copy current tensor to a new tensor
812 
813  \future Return a scalar if possible as a rank 1 tensor with
814  1 element.
815  */
816  tensor<data_t> rearrange_and_copy(std::vector<index_spec> spec,
817  int verbose=0, bool err_on_fail=true) {
818 
819  // Old rank and new rank (computed later)
820  size_t rank_old=this->rk;
821  size_t rank_new=0;
822 
823  // Size of the new indices
824  std::vector<size_t> size_new;
825 
826  // Reorganize the index specifications by both the old
827  // and new indexing system
828  std::vector<index_spec> spec_old(rank_old);
829  std::vector<index_spec> spec_new;
830 
831  // Number of elements to sum over
832  size_t n_sum_loop=1;
833 
834  // Size of sums
835  std::vector<size_t> sum_sizes;
836 
837  // Collect the statistics on the transformation
838  for(size_t i=0;i<spec.size();i++) {
839  if (spec[i].type==index_spec::index ||
840  spec[i].type==index_spec::reverse) {
841  if (spec[i].ix1>=rank_old) {
842  if (err_on_fail) {
843  O2SCL_ERR2("Index too large (index,reverse) in ",
844  "tensor::rearrange_and_copy().",o2scl::exc_einval);
845  } else {
846  if (verbose>0) {
847  std::cout << "Index " << spec[i].ix1
848  << " too large (index,reverse) in "
849  << "tensor in tensor::rearrange_and_copy()."
850  << std::endl;
851  }
852  return tensor<data_t>();
853  }
854  }
855  size_new.push_back(this->size[spec[i].ix1]);
856  // Use ix1 to store the destination index (which is
857  // at this point equal to rank_new)
858  spec_old[spec[i].ix1]=index_spec(spec[i].type,rank_new);
859  spec_new.push_back(index_spec(spec[i].type,spec[i].ix1));
860  rank_new++;
861  } else if (spec[i].type==index_spec::range) {
862  if (verbose>2) {
863  std::cout << "In range " << spec[i].ix1 << " "
864  << spec[i].ix2 << " " << spec[i].ix3 << std::endl;
865  }
866  if (spec[i].ix1>=rank_old ||
867  spec[i].ix2>=this->size[spec[i].ix1] ||
868  spec[i].ix3>=this->size[spec[i].ix1]) {
869  if (err_on_fail) {
870  O2SCL_ERR2("Index too large (range) in ",
871  "tensor::rearrange_and_copy().",o2scl::exc_einval);
872  } else {
873  if (verbose>0) {
874  std::cout << "Index " << spec[i].ix1 << " "
875  << spec[i].ix2 << " " << spec[i].ix3
876  << " too large (range) in "
877  << "tensor in tensor::rearrange_and_copy()."
878  << std::endl;
879  }
880  return tensor<data_t>();
881  }
882  }
883  if (spec[i].ix3>spec[i].ix2) {
884  size_new.push_back(spec[i].ix3-spec[i].ix2+1);
885  } else {
886  size_new.push_back(spec[i].ix2-spec[i].ix3+1);
887  }
888  // Use ix1 to store the destination index (which is
889  // at this point equal to rank_new)
890  spec_old[spec[i].ix1]=
891  index_spec(spec[i].type,rank_new,spec[i].ix2,spec[i].ix3);
892  spec_new.push_back
893  (index_spec(spec[i].type,spec[i].ix1,
894  spec[i].ix2,spec[i].ix3));
895  rank_new++;
896  if (verbose>2) {
897  std::cout << "Out range " << size_new[size_new.size()-1]
898  << std::endl;
899  }
900  } else if (spec[i].type==index_spec::trace) {
901  if (spec[i].ix1>=rank_old || spec[i].ix2>=rank_old) {
902  if (err_on_fail) {
903  O2SCL_ERR2("Index too large (trace) in ",
904  "tensor::rearrange_and_copy().",o2scl::exc_einval);
905  } else {
906  if (verbose>0) {
907  std::cout << "Indices " << spec[i].ix1 << " or "
908  << spec[i].ix2 << " too large (trace) in "
909  << "tensor in tensor::rearrange_and_copy()."
910  << std::endl;
911  }
912  return tensor<data_t>();
913  }
914  }
915  if (size[spec[i].ix1]<size[spec[i].ix2]) {
916  n_sum_loop*=size[spec[i].ix1];
917  sum_sizes.push_back(size[spec[i].ix1]);
918  } else {
919  n_sum_loop*=size[spec[i].ix2];
920  sum_sizes.push_back(size[spec[i].ix2]);
921  }
922  // We set the values of ix1 and ix2 so that ix2
923  // always refers to the other index being traced over
924  spec_old[spec[i].ix1]=index_spec(spec[i].type,
925  spec[i].ix1,spec[i].ix2);
926  spec_old[spec[i].ix2]=index_spec(spec[i].type,
927  spec[i].ix2,spec[i].ix1);
928  } else if (spec[i].type==index_spec::sum) {
929  if (spec[i].ix1>=rank_old) {
930  if (err_on_fail) {
931  O2SCL_ERR2("Index too large (sum) in ",
932  "tensor::rearrange_and_copy().",o2scl::exc_einval);
933  } else {
934  if (verbose>0) {
935  std::cout << "Index " << spec[i].ix1
936  << " too large (sum) in "
937  << "tensor in tensor::rearrange_and_copy()."
938  << std::endl;
939  }
940  return tensor<data_t>();
941  }
942  }
943  n_sum_loop*=size[spec[i].ix1];
944  sum_sizes.push_back(size[spec[i].ix1]);
945  spec_old[spec[i].ix1]=index_spec(spec[i].type,
946  spec[i].ix1,spec[i].ix2,0);
947  } else if (spec[i].type==index_spec::fixed) {
948  if (spec[i].ix1>=rank_old ||
949  spec[i].ix2>=this->size[spec[i].ix1]) {
950  if (err_on_fail) {
951  O2SCL_ERR2("Index too large (fixed) in ",
952  "tensor::rearrange_and_copy().",o2scl::exc_einval);
953  } else {
954  if (verbose>0) {
955  std::cout << "Index too large (fixed) in "
956  << "tensor in tensor::rearrange_and_copy()."
957  << std::endl;
958  }
959  return tensor<data_t>();
960  }
961  }
962  // Use ix1 to store the destination index (which is
963  // at this point equal to rank_new)
964  spec_old[spec[i].ix1]=index_spec(spec[i].type,
965  rank_new,spec[i].ix2);
966  } else {
967  if (err_on_fail) {
968  O2SCL_ERR2("Index specification type not allowed in ",
969  "tensor::rearrange_and_copy().",o2scl::exc_einval);
970  } else {
971  if (verbose>0) {
972  std::cout << "Index specification type not allowed in "
973  << "tensor::rearrange_and_copy()." << std::endl;
974  }
975  return tensor<data_t>();
976  }
977  }
978  }
979  size_t n_sums=sum_sizes.size();
980 
981  // Call the error handler if the input is invalid
982  if (rank_new==0) {
983  if (err_on_fail) {
984  O2SCL_ERR2("Zero new indices in ",
985  "tensor::rearrange_and_copy().",o2scl::exc_einval);
986  } else {
987  if (verbose>0) {
988  std::cout << "Zero new indices in "
989  << "tensor::rearrange_and_copy()." << std::endl;
990  }
991  return tensor<data_t>();
992  }
993  }
994 
995  for(size_t i=0;i<rank_old;i++) {
996  if (spec_old[i].type==index_spec::empty) {
997  if (err_on_fail) {
998  O2SCL_ERR2("Not all indices accounted for in ",
999  "tensor::rearrange_and_copy().",o2scl::exc_einval);
1000  } else {
1001  if (verbose>0) {
1002  std::cout << "Index " << i << " not accounted for in "
1003  << "tensor::rearrange_and_copy()." << std::endl;
1004  }
1005  return tensor<data_t>();
1006  }
1007  }
1008  }
1009 
1010  // Verbose output if necessary
1011  if (verbose>0) {
1012  std::cout << "Using a " << rank_old << " rank tensor to create a new "
1013  << rank_new << " rank tensor." << std::endl;
1014  }
1015  if (verbose>1) {
1016  for(size_t i=0;i<rank_old;i++) {
1017  std::cout << "Old index " << i;
1018  if (spec_old[i].type==index_spec::index) {
1019  std::cout << " is being remapped to new index " << spec_old[i].ix1
1020  << "." << std::endl;
1021  } else if (spec_old[i].type==index_spec::range) {
1022  std::cout << " is being remapped to new index " << spec_old[i].ix1
1023  << " with a range from " << spec_old[i].ix2
1024  << " to " << spec_old[i].ix3 << "." << std::endl;
1025  } else if (spec_old[i].type==index_spec::reverse) {
1026  std::cout << " is being reversed and remapped to new index "
1027  << spec_old[i].ix1 << "." << std::endl;
1028  } else if (spec_old[i].type==index_spec::trace) {
1029  std::cout << " is being traced with index "
1030  << spec_old[i].ix2 << "." << std::endl;
1031  } else if (spec_old[i].type==index_spec::sum) {
1032  std::cout << " is being summed." << std::endl;
1033  } else if (spec_old[i].type==index_spec::fixed) {
1034  std::cout << " is being fixed to " << spec_old[i].ix2
1035  << "." << std::endl;
1036  }
1037  }
1038  for(size_t i=0;i<rank_new;i++) {
1039  std::cout << "New index " << i;
1040  if (spec_new[i].type==index_spec::index) {
1041  std::cout << " was remapped from old index " << spec_new[i].ix1
1042  << "." << std::endl;
1043  } else if (spec_new[i].type==index_spec::range) {
1044  std::cout << " was remapped from old index " << spec_new[i].ix1
1045  << " using range from " << spec_new[i].ix2 << " to "
1046  << spec_new[i].ix3 << "." << std::endl;
1047  } else if (spec_new[i].type==index_spec::reverse) {
1048  std::cout << " was reversed and remapped from old index "
1049  << spec_new[i].ix1 << "." << std::endl;
1050  }
1051  }
1052  }
1053 
1054  // Create the new tensor object
1055  tensor<data_t> t_new(rank_new,size_new);
1056 
1057  // Index arrays
1058  std::vector<size_t> ix_new(rank_new);
1059  std::vector<size_t> ix_old(rank_old);
1060  std::vector<size_t> sum_ix(n_sums);
1061 
1062  // Loop over the new tensor object
1063  for(size_t i=0;i<t_new.total_size();i++) {
1064 
1065  // Find the location in the new tensor object
1066  t_new.unpack_index(i,ix_new);
1067 
1068  // Determine the location in the old tensor object
1069  for(size_t j=0;j<rank_old;j++) {
1070  if (spec_old[j].type==index_spec::index) {
1071  ix_old[j]=ix_new[spec_old[j].ix1];
1072  } else if (spec_old[j].type==index_spec::range) {
1073  if (spec_old[j].ix2<spec_old[j].ix3) {
1074  ix_old[j]=ix_new[spec_old[j].ix1]+spec_old[j].ix2;
1075  } else {
1076  ix_old[j]=spec_old[j].ix2-ix_new[spec_old[j].ix1];
1077  }
1078  } else if (spec_old[j].type==index_spec::reverse) {
1079  ix_old[j]=get_size(j)-1-ix_new[spec_old[j].ix1];
1080  } else if (spec_old[j].type==index_spec::fixed) {
1081  ix_old[j]=spec_old[j].ix2;
1082  }
1083  }
1084 
1085  data_t val=0;
1086 
1087  for(size_t j=0;j<n_sum_loop;j++) {
1088 
1089  // This code is similar to tensor::unpack_index(), it unpacks
1090  // the index j to the indices which we are summing over.
1091  size_t j2=j, sub_size;
1092  for(size_t k=0;k<n_sums;k++) {
1093  if (k==n_sums-1) {
1094  sum_ix[k]=j2;
1095  } else {
1096  sub_size=1;
1097  for(size_t kk=k+1;kk<n_sums;kk++) sub_size*=sum_sizes[kk];
1098  sum_ix[k]=j2/sub_size;
1099  j2-=sub_size*(j2/sub_size);
1100  }
1101  }
1102  if (verbose>2) {
1103  std::cout << "n_sum_loop: " << n_sum_loop << " n_sums: "
1104  << n_sums << " sum_sizes: ";
1105  vector_out(std::cout,sum_sizes,true);
1106  std::cout << "j: " << j << " sum_ix: ";
1107  vector_out(std::cout,sum_ix,true);
1108  }
1109 
1110  // Remap from sum_ix to ix_old
1111  size_t cnt=0;
1112  for(size_t k=0;k<rank_old;k++) {
1113  if (spec_old[k].type==index_spec::sum) {
1114  if (cnt>=sum_ix.size()) {
1115  std::cout << "X: " << cnt << " " << sum_ix.size() << std::endl;
1116  O2SCL_ERR2("Bad sync 1 in sum_ix in ",
1117  "tensor::rearrange_and_copy()",o2scl::exc_esanity);
1118  }
1119  ix_old[k]=sum_ix[cnt];
1120  cnt++;
1121  } else if (spec_old[k].type==index_spec::trace &&
1122  spec_old[k].ix1<spec_old[k].ix2) {
1123  if (cnt>=sum_ix.size()) {
1124  std::cout << "X: " << cnt << " " << sum_ix.size() << std::endl;
1125  O2SCL_ERR2("Bad sync 2 in sum_ix in ",
1126  "tensor::rearrange_and_copy()",o2scl::exc_esanity);
1127  }
1128  ix_old[spec_old[k].ix1]=sum_ix[cnt];
1129  ix_old[spec_old[k].ix2]=sum_ix[cnt];
1130  cnt++;
1131  }
1132  }
1133 
1134  if (verbose>2) {
1135  std::cout << "Here old: ";
1136  vector_out(std::cout,ix_old,true);
1137  std::cout << "Here new: ";
1138  vector_out(std::cout,ix_new,true);
1139  }
1140  val+=get(ix_old);
1141 
1142  }
1143 
1144  // Set the new point by performing the linear interpolation
1145  t_new.set(ix_new,val);
1146  }
1147 
1148  return t_new;
1149  }
1150 
1151  };
1152 
1153  /** \brief Rank 1 tensor
1154  */
1155  template<class data_t=double, class vec_t=std::vector<data_t>,
1156  class vec_size_t=std::vector<size_t> > class tensor1 :
1157  public tensor<data_t,vec_t,vec_size_t> {
1158 
1159  public:
1160 
1161  /// Create an empty tensor
1162  tensor1() : tensor<data_t,vec_t,vec_size_t>() {}
1163 
1164  /// Create a rank 1 tensory of size \c sz
1165  tensor1(size_t sz) : tensor<data_t,vec_t,vec_size_t>() {
1166  vec_size_t sizex(1);
1167  sizex[0]=sz;
1168  this->resize(1,sizex);
1169  }
1170 
1171  /// \name Method to check for valid object
1172  //@{
1173  /** \brief Check that the \ref o2scl::tensor1 object is valid
1174  */
1175  void is_valid() const {
1177  if (this->rk>1) {
1178  O2SCL_ERR2("Rank is neither 0 nor 1 in ",
1179  "tensor1::is_valid().",
1181 
1182  }
1183  return;
1184  }
1185  //@}
1186 
1187  /// \name Specialized get and set functions
1188  //@{
1189  /// Get the element indexed by \c ix
1190  data_t &get(size_t ix) {
1192  }
1193 
1194  /// Get the element indexed by \c ix
1195  const data_t &get(size_t ix) const {
1197  }
1198 
1199  /// Set the element indexed by \c index to value \c val
1200  void set(size_t index, data_t val)
1201  { tensor<data_t,vec_t,vec_size_t>::set(&index,val); }
1202 
1203  /** \brief Set the element indexed by \c index to value \c val
1204 
1205  (We have to explicitly provide this version since the set()
1206  function is overloaded in this child of \ref tensor.)
1207  */
1208  template<class size_vec_t>
1209  void set(const size_vec_t &index, data_t val) {
1211  }
1212  //@}
1213 
1214  /// \name Specialized operator functions
1215  //@{
1216  /// Get an element using array-like indexing
1217  data_t &operator[](size_t ix) { return this->data[ix]; }
1218 
1219  /// Get an element using array-like indexing (const version)
1220  const data_t &operator[](size_t ix) const { return this->data[ix]; }
1221 
1222  /// Get an element using operator()
1223  data_t &operator()(size_t ix) { return this->data[ix]; }
1224 
1225  /// Get an element using operator() (const version)
1226  const data_t &operator()(size_t ix) const { return this->data[ix]; }
1227  //@}
1228 
1229  };
1230 
1231  /** \brief Rank 2 tensor
1232  */
1233  template<class data_t=double, class vec_t=std::vector<data_t>,
1234  class vec_size_t=std::vector<size_t> > class tensor2 :
1235  public tensor<data_t,vec_t,vec_size_t> {
1236 
1237  public:
1238 
1239  /// Create an empty tensor
1240  tensor2() : tensor<data_t,vec_t,vec_size_t>() {}
1241 
1242  /// Create a rank 2 tensor of size \c (sz,sz2)
1243  tensor2(size_t sz, size_t sz2) : tensor<data_t,vec_t,vec_size_t>() {
1244  this->rk=2;
1245  this->size.resize(2);
1246  this->size[0]=sz;
1247  this->size[1]=sz2;
1248  size_t tot=sz*sz2;
1249  this->data.resize(tot);
1250  }
1251 
1252  /// \name Method to check for valid object
1253  //@{
1254  /** \brief Check that the \ref o2scl::tensor2 object is valid
1255  */
1256  void is_valid() const {
1258  if (this->rk!=0 && this->rk!=2) {
1259  O2SCL_ERR2("Rank is neither 0 nor 2 in ",
1260  "tensor2::is_valid().",
1262 
1263  }
1264  return;
1265  }
1266  //@}
1267 
1268  /// \name Specialized get and set functions
1269  //@{
1270  /// Get the element indexed by \c (ix1,ix2)
1271  data_t &get(size_t ix1, size_t ix2) {
1272  size_t sz[2]={ix1,ix2};
1274  }
1275 
1276  /// Get the element indexed by \c (ix1,ix2)
1277  const data_t &get(size_t ix1, size_t ix2) const {
1278  size_t sz[2]={ix1,ix2};
1280  }
1281 
1282  /// Set the element indexed by \c (ix1,ix2) to value \c val
1283  void set(size_t ix1, size_t ix2, data_t val) {
1284  size_t sz[2]={ix1,ix2};
1286  return;
1287  }
1288 
1289  /** \brief Set the element indexed by \c index to value \c val
1290 
1291  (We have to explicitly provide this version since the set()
1292  function is overloaded in this child of \ref tensor.)
1293  */
1294  template<class size_vec_t>
1295  void set(const size_vec_t &index, data_t val) {
1297  return;
1298  }
1299 
1300  /// Get the element indexed by \c (ix1,ix2)
1301  data_t &operator()(size_t ix, size_t iy)
1302  { return this->data[ix*this->size[1]+iy]; }
1303 
1304  /// Get the element indexed by \c (ix1,ix2) (const version)
1305  const data_t &operator()(size_t ix, size_t iy) const
1306  { return this->data[ix*this->size[1]+iy]; }
1307  //@}
1308  };
1309 
1310  /** \brief Rank 3 tensor
1311  */
1312  template<class data_t=double, class vec_t=std::vector<data_t>,
1313  class vec_size_t=std::vector<size_t> > class tensor3 :
1314  public tensor<data_t,vec_t,vec_size_t> {
1315 
1316  public:
1317 
1318  /// Create an empty tensor
1319  tensor3() : tensor<data_t,vec_t,vec_size_t>() {}
1320 
1321  /// Create a rank 3 tensor of size \c (sz,sz2,sz3)
1322  tensor3(size_t sz, size_t sz2, size_t sz3) :
1323  tensor<data_t,vec_t,vec_size_t>() {
1324  this->rk=3;
1325  this->size.resize(3);
1326  this->size[0]=sz;
1327  this->size[1]=sz2;
1328  this->size[2]=sz3;
1329  size_t tot=sz*sz2*sz3;
1330  this->data.resize(tot);
1331  }
1332 
1333  /// \name Method to check for valid object
1334  //@{
1335  /** \brief Check that the \ref o2scl::tensor3 object is valid
1336  */
1337  void is_valid() const {
1339  if (this->rk!=0 && this->rk!=3) {
1340  O2SCL_ERR2("Rank is neither 0 nor 3 in ",
1341  "tensor3::is_valid().",
1343 
1344  }
1345  return;
1346  }
1347  //@}
1348 
1349  /// \name Specialized get and set functions
1350  //@{
1351  /// Get the element indexed by \c (ix1,ix2,ix3)
1352  data_t &get(size_t ix1, size_t ix2, size_t ix3) {
1353  size_t sz[3]={ix1,ix2,ix3};
1355  }
1356 
1357  /// Get the element indexed by \c (ix1,ix2,ix3)
1358  const data_t &get(size_t ix1, size_t ix2, size_t ix3) const {
1359  size_t sz[3]={ix1,ix2,ix3};
1361  }
1362 
1363  /// Set the element indexed by \c (ix1,ix2,ix3) to value \c val
1364  void set(size_t ix1, size_t ix2, size_t ix3, data_t val) {
1365  size_t sz[3]={ix1,ix2,ix3};
1367  return;
1368  }
1369 
1370  /** \brief Set the element indexed by \c index to value \c val
1371 
1372  (We have to explicitly provide this version since the set()
1373  function is overloaded in this child of \ref tensor.)
1374  */
1375  template<class size_vec_t>
1376  void set(const size_vec_t &index, data_t val) {
1378  return;
1379  }
1380  //@}
1381 
1382  };
1383 
1384  /** \brief Rank 4 tensor
1385  */
1386  template<class data_t=double, class vec_t=std::vector<data_t>,
1387  class vec_size_t=std::vector<size_t> > class tensor4 :
1388  public tensor<data_t,vec_t,vec_size_t> {
1389 
1390  public:
1391 
1392  /// Create an empty tensor
1393  tensor4() : tensor<data_t,vec_t,vec_size_t>() {}
1394 
1395  /// Create a rank 4 tensor of size \c (sz,sz2,sz3,sz4)
1396  tensor4(size_t sz, size_t sz2, size_t sz3, size_t sz4) :
1397  tensor<data_t,vec_t,vec_size_t>() {
1398  this->rk=4;
1399  this->size.resize(4);
1400  this->size[0]=sz;
1401  this->size[1]=sz2;
1402  this->size[2]=sz3;
1403  this->size[3]=sz4;
1404  size_t tot=sz*sz2*sz3*sz4;
1405  this->data.resize(tot);
1406  }
1407 
1408  /// \name Method to check for valid object
1409  //@{
1410  /** \brief Check that the \ref o2scl::tensor4 object is valid
1411  */
1412  void is_valid() const {
1414  if (this->rk!=0 && this->rk!=4) {
1415  O2SCL_ERR2("Rank is neither 0 nor 4 in ",
1416  "tensor4::is_valid().",
1418 
1419  }
1420  return;
1421  }
1422  //@}
1423 
1424  /// \name Specialized get and set functions
1425  //@{
1426  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1427  data_t &get(size_t ix1, size_t ix2, size_t ix3, size_t ix4) {
1428  size_t sz[4]={ix1,ix2,ix3,ix4};
1430  }
1431 
1432  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1433  const data_t &get(size_t ix1, size_t ix2, size_t ix3,
1434  size_t ix4) const {
1435  size_t sz[4]={ix1,ix2,ix3,ix4};
1437  }
1438 
1439  /// Set the element indexed by \c (ix1,ix2,ix3,ix4) to value \c val
1440  void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4,
1441  data_t val) {
1442  size_t sz[4]={ix1,ix2,ix3,ix4};
1444  return;
1445  }
1446 
1447  /** \brief Set the element indexed by \c index to value \c val
1448 
1449  (We have to explicitly provide this version since the set()
1450  function is overloaded in this child of \ref tensor.)
1451  */
1452  template<class size_vec_t>
1453  void set(const size_vec_t &index, data_t val) {
1455  return;
1456  }
1457  //@}
1458  };
1459 
1460  /// \name Tensor functions in src/base/tensor.h
1461  //@{
1462  /** \brief Output a tensor to a stream
1463  */
1464  template<class tensor_t>
1465  void tensor_out(std::ostream &os, tensor_t &t, bool pretty=true) {
1466 
1467  if (pretty) {
1468 
1469  size_t rk=t.get_rank();
1470  os << "rank: " << rk << " sizes: ";
1471  for(size_t i=0;i<rk;i++) {
1472  os << t.get_size(i) << " ";
1473  }
1474  os << std::endl;
1475  auto &data=t.get_data();
1476  std::vector<size_t> ix(rk);
1477  std::vector<std::string> sv, sv_out;
1478  for(size_t i=0;i<t.total_size();i++) {
1479  t.unpack_index(i,ix);
1480  std::string tmp="(";
1481  for(size_t j=0;j<rk;j++) {
1482  if (j!=rk-1) {
1483  tmp+=o2scl::szttos(ix[j])+",";
1484  } else {
1485  tmp+=o2scl::szttos(ix[j]);
1486  }
1487  }
1488  tmp+="): "+o2scl::dtos(data[i]);
1489  sv.push_back(tmp);
1490  }
1491  screenify(sv.size(),sv,sv_out);
1492  for(size_t k=0;k<sv_out.size();k++) {
1493  os << sv_out[k] << std::endl;
1494  }
1495 
1496  } else {
1497 
1498  size_t rk=t.get_rank();
1499  os << rk << " ";
1500  for(size_t i=0;i<rk;i++) {
1501  os << t.get_size(i) << " ";
1502  }
1503  os << std::endl;
1504  auto &data=t.get_data();
1505  for(size_t i=0;i<t.total_size();i++) {
1506  os << data[i] << " ";
1507  if (i%10==9) os << std::endl;
1508  }
1509  os << std::endl;
1510 
1511  }
1512 
1513  return;
1514  }
1515 
1516  /** \brief Compare two tensors for equality
1517  */
1518  template<class data_t, class vec_t, class vec_size_t>
1520  const tensor<data_t,vec_t,vec_size_t> &t2) {
1521  if (t1.get_rank()!=t2.get_rank()) return false;
1522  for(size_t i=0;i<t1.get_rank();i++) {
1523  if (t1.get_size(i)!=t2.get_size(i)) return false;
1524  }
1525  const vec_t &v1=t1.get_data();
1526  const vec_t &v2=t2.get_data();
1527  for(size_t i=0;i<t1.total_size();i++) {
1528  if (v1[i]!=v2[i]) return false;
1529  }
1530  return true;
1531  }
1532  //@}
1533 
1534 #ifndef DOXYGEN_NO_O2NS
1535 }
1536 #endif
1537 
1538 #endif
1539 
1540 
1541 
o2scl::tensor
Tensor class with arbitrary dimensions.
Definition: tensor.h:226
o2scl::tensor::get_rank
size_t get_rank() const
Return the rank of the tensor.
Definition: tensor.h:576
o2scl::tensor4::get
data_t & get(size_t ix1, size_t ix2, size_t ix3, size_t ix4)
Get the element indexed by (ix1,ix2,ix3,ix4)
Definition: tensor.h:1427
o2scl::index_spec::index_spec
index_spec(size_t typ, size_t i1, size_t i2=0, size_t i3=0, double v1=0.0, double v2=0.0, double v3=0.0)
Explicit full constructor.
Definition: tensor.h:109
o2scl::tensor1::operator[]
const data_t & operator[](size_t ix) const
Get an element using array-like indexing (const version)
Definition: tensor.h:1220
o2scl::tensor::tensor
tensor(size_t rank, const size_vec_t &dim)
Create a tensor of rank rank with sizes given in dim.
Definition: tensor.h:260
o2scl::tensor4::tensor4
tensor4()
Create an empty tensor.
Definition: tensor.h:1393
o2scl::index_spec::index
static const size_t index
Retain an index.
Definition: tensor.h:78
o2scl::tensor1::tensor1
tensor1()
Create an empty tensor.
Definition: tensor.h:1162
o2scl::tensor3::set
void set(const size_vec_t &index, data_t val)
Set the element indexed by index to value val.
Definition: tensor.h:1376
o2scl::tensor1::tensor1
tensor1(size_t sz)
Create a rank 1 tensory of size sz.
Definition: tensor.h:1165
o2scl::tensor::min
void min(vec_size_t &index, data_t &val)
Compute the index of the minimum value in the tensor and return the minimum.
Definition: tensor.h:685
o2scl::table3d::set
void set(size_t ix, size_t iy, std::string name, double val)
Set element in slice name at location ix,iy to value val.
o2scl::tensor1::operator()
data_t & operator()(size_t ix)
Get an element using operator()
Definition: tensor.h:1223
o2scl::tensor::max_value
data_t max_value()
Compute the maximum value in the tensor.
Definition: tensor.h:694
o2scl::dtos
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Definition: string_conv.h:77
boost::numeric::ublas::vector
The default vector type from uBlas.
Definition: vector.h:3735
o2scl::index_spec::sum
static const size_t sum
Sum over an index.
Definition: tensor.h:82
o2scl::index_spec::reverse
static const size_t reverse
Reverse an index.
Definition: tensor.h:86
o2scl::index_spec::type
size_t type
Type of specification.
Definition: tensor.h:59
o2scl::tensor::is_valid
void is_valid() const
Check that the o2scl::tensor object is valid.
Definition: tensor.h:291
o2scl::tensor4::set
void set(const size_vec_t &index, data_t val)
Set the element indexed by index to value val.
Definition: tensor.h:1453
o2scl::ix_range
index_spec ix_range(size_t ix, size_t start, size_t end)
Index covers a range of values.
o2scl::ix_grid
index_spec ix_grid(size_t ix, double begin, double end, size_t n_bins, bool log=false)
Interpolate grid with fixed number of bins into index ix (for o2scl::tensor_grid only)
o2scl::tensor2::set
void set(size_t ix1, size_t ix2, data_t val)
Set the element indexed by (ix1,ix2) to value val.
Definition: tensor.h:1283
o2scl::tensor::set
void set(const size_vec_t &index, data_t val)
Set the element indexed by index to value val.
Definition: tensor.h:360
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::operator==
bool operator==(const tensor< data_t, vec_t, vec_size_t > &t1, const tensor< data_t, vec_t, vec_size_t > &t2)
Compare two tensors for equality.
Definition: tensor.h:1519
o2scl::tensor2::get
const data_t & get(size_t ix1, size_t ix2) const
Get the element indexed by (ix1,ix2)
Definition: tensor.h:1277
o2scl::tensor1::operator()
const data_t & operator()(size_t ix) const
Get an element using operator() (const version)
Definition: tensor.h:1226
o2scl::index_spec::fixed
static const size_t fixed
Fix the value of an index.
Definition: tensor.h:80
o2scl::index_spec::grid
static const size_t grid
Interpolate a value to set a new grid (fixed bin number)
Definition: tensor.h:92
o2scl::tensor::min_index
void min_index(vec_size_t &index)
Compute the index of the minimum value in the tensor.
Definition: tensor.h:676
o2scl::tensor1::operator[]
data_t & operator[](size_t ix)
Get an element using array-like indexing.
Definition: tensor.h:1217
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::tensor::max_index
void max_index(vec_size_t &index)
Compute the index of the maximum value in the tensor.
Definition: tensor.h:700
o2scl::tensor::get
const data_t & get(const size_vec_t &index) const
Get a const reference to the element indexed by index.
Definition: tensor.h:444
o2scl::tensor::get
data_t & get(const size_vec_t &index)
Get the element indexed by index.
Definition: tensor.h:412
o2scl::tensor3::get
const data_t & get(size_t ix1, size_t ix2, size_t ix3) const
Get the element indexed by (ix1,ix2,ix3)
Definition: tensor.h:1358
o2scl::index_spec
Index specification.
Definition: tensor.h:54
o2scl::tensor::vector_slice
ubvector_slice vector_slice(size_t ix, const size_vec_t &index)
Fix all but one index to create a vector.
Definition: tensor.h:504
o2scl::index_spec::val1
double val1
First double argument.
Definition: tensor.h:67
o2scl::tensor::size
vec_size_t size
A rank-sized vector of the sizes of each dimension.
Definition: tensor.h:238
o2scl::ix_gridw
index_spec ix_gridw(size_t ix, double begin, double end, double width, bool log=false)
Interpolate grid with fixed bin width into index ix (for o2scl::tensor_grid only)
o2scl::tensor::get_size
size_t get_size(size_t i) const
Returns the size of the ith index.
Definition: tensor.h:579
o2scl::tensor2::tensor2
tensor2()
Create an empty tensor.
Definition: tensor.h:1240
o2scl::tensor::pack_indices
size_t pack_indices(const size_vec_t &index)
Pack the indices into a single vector index.
Definition: tensor.h:614
o2scl::tensor3::tensor3
tensor3()
Create an empty tensor.
Definition: tensor.h:1319
o2scl::tensor::get_size_arr
const vec_size_t & get_size_arr() const
Return the full vector of sizes.
Definition: tensor.h:590
o2scl::tensor2::operator()
const data_t & operator()(size_t ix, size_t iy) const
Get the element indexed by (ix1,ix2) (const version)
Definition: tensor.h:1305
o2scl::index_spec::trace
static const size_t trace
Perform a trace (sum over two indices)
Definition: tensor.h:84
o2scl::ix_index
index_spec ix_index(size_t ix)
Choose an index.
o2scl::tensor4::is_valid
void is_valid() const
Check that the o2scl::tensor4 object is valid.
Definition: tensor.h:1412
o2scl::ix_interp
index_spec ix_interp(size_t ix, double v)
Interpolate value v into index ix (for o2scl::tensor_grid only)
o2scl::index_spec::ix2
size_t ix2
Second argument.
Definition: tensor.h:63
o2scl::tensor::get_data
const vec_t & get_data() const
Return the full data vector.
Definition: tensor.h:595
o2scl::screenify
void screenify(size_t nin, const string_arr_t &in_cols, std::vector< std::string > &out_cols, size_t max_size=80)
Reformat the columns for output of width size.
Definition: misc.h:200
o2scl::index_spec::ix1
size_t ix1
First argument.
Definition: tensor.h:61
o2scl::tensor::total_sum
double total_sum() const
Return the sum over every element in the tensor.
Definition: tensor.h:748
o2scl::tensor2
Rank 2 tensor.
Definition: tensor.h:1234
o2scl::tensor4
Rank 4 tensor.
Definition: tensor.h:1387
o2scl::tensor1::get
const data_t & get(size_t ix) const
Get the element indexed by ix.
Definition: tensor.h:1195
o2scl::tensor1::get
data_t & get(size_t ix)
Get the element indexed by ix.
Definition: tensor.h:1190
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::tensor::minmax_index
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:724
o2scl::tensor2::get
data_t & get(size_t ix1, size_t ix2)
Get the element indexed by (ix1,ix2)
Definition: tensor.h:1271
o2scl::tensor4::tensor4
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:1396
o2scl::tensor3::get
data_t & get(size_t ix1, size_t ix2, size_t ix3)
Get the element indexed by (ix1,ix2,ix3)
Definition: tensor.h:1352
o2scl::tensor3::set
void set(size_t ix1, size_t ix2, size_t ix3, data_t val)
Set the element indexed by (ix1,ix2,ix3) to value val.
Definition: tensor.h:1364
o2scl::tensor2::set
void set(const size_vec_t &index, data_t val)
Set the element indexed by index to value val.
Definition: tensor.h:1295
o2scl::tensor1::is_valid
void is_valid() const
Check that the o2scl::tensor1 object is valid.
Definition: tensor.h:1175
o2scl::tensor3::tensor3
tensor3(size_t sz, size_t sz2, size_t sz3)
Create a rank 3 tensor of size (sz,sz2,sz3)
Definition: tensor.h:1322
o2scl::ix_trace
index_spec ix_trace(size_t ix, size_t ix2)
Perform a trace over indices ix and ix2.
o2scl::table3d::get
double & get(size_t ix, size_t iy, std::string name)
Get element in slice name at location ix,iy
o2scl::index_spec::ix3
size_t ix3
Third argument.
Definition: tensor.h:65
o2scl::tensor::clear
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor.h:348
o2scl::tensor::rk
size_t rk
Rank.
Definition: tensor.h:241
o2scl::tensor2::tensor2
tensor2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
Definition: tensor.h:1243
o2scl::tensor::min_value
data_t min_value()
Compute the minimum value in the tensor.
Definition: tensor.h:670
o2scl::table3d::set_xy
void set_xy(std::string x_name, size_t nx, const vec_t &x, std::string y_name, size_t ny, const vec2_t &y)
Initialize the x-y grid.
Definition: table3d.h:120
o2scl::tensor::convert_table3d_sum
void convert_table3d_sum(size_t ix_x, size_t ix_y, table3d &tab, std::string x_name="x", std::string y_name="y", std::string slice_name="z")
Convert to a o2scl::table3d object by summing over all but two indices.
Definition: tensor.h:764
o2scl::tensor1
Rank 1 tensor.
Definition: tensor.h:1156
o2scl::tensor::tensor
tensor()
Create an empty tensor with zero rank.
Definition: tensor.h:248
o2scl::tensor1::set
void set(const size_vec_t &index, data_t val)
Set the element indexed by index to value val.
Definition: tensor.h:1209
o2scl::tensor3
Rank 3 tensor.
Definition: tensor.h:1313
o2scl::tensor2::operator()
data_t & operator()(size_t ix, size_t iy)
Get the element indexed by (ix1,ix2)
Definition: tensor.h:1301
o2scl::exc_eindex
@ exc_eindex
Invalid index for array or matrix.
Definition: err_hnd.h:123
o2scl::index_spec::range
static const size_t range
Choose a new range for an index.
Definition: tensor.h:88
o2scl::tensor4::set
void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4, data_t val)
Set the element indexed by (ix1,ix2,ix3,ix4) to value val.
Definition: tensor.h:1440
o2scl::index_spec::val2
double val2
Second double argument.
Definition: tensor.h:69
o2scl::index_spec::val3
double val3
Third double argument.
Definition: tensor.h:71
o2scl::exc_esanity
@ exc_esanity
sanity check failed - shouldn't happen
Definition: err_hnd.h:65
o2scl::interp
Interpolation class for general vectors.
Definition: interp.h:1654
o2scl::table3d::set_slice_all
void set_slice_all(std::string name, double val)
Set all of the values in slice name to val.
o2scl::tensor::minmax
void minmax(vec_size_t &index, size_t &index_min, data_t &min, size_t &index_max, data_t &max)
Compute the indices and values of the maximum and minimum in the tensor.
Definition: tensor.h:736
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::tensor::rearrange_and_copy
tensor< data_t > rearrange_and_copy(std::vector< index_spec > spec, int verbose=0, bool err_on_fail=true)
Rearrange, sum and copy current tensor to a new tensor.
Definition: tensor.h:816
o2scl::table3d::get_size
void get_size(size_t &nx, size_t &ny) const
Get the size of the slices.
o2scl::tensor::data
vec_t data
The data.
Definition: tensor.h:235
o2scl::ix_reverse
index_spec ix_reverse(size_t ix)
Reverse index ix.
o2scl::tensor_out
void tensor_out(std::ostream &os, tensor_t &t, bool pretty=true)
Output a tensor to a stream.
Definition: tensor.h:1465
o2scl::tensor2::is_valid
void is_valid() const
Check that the o2scl::tensor2 object is valid.
Definition: tensor.h:1256
o2scl::index_spec::gridw
static const size_t gridw
Interpolate a value to set a new grid (fixed bin width)
Definition: tensor.h:94
o2scl::tensor::resize
void resize(size_t rank, const size_vec_t &dim)
Resize the tensor to rank rank with sizes given in dim.
Definition: tensor.h:549
o2scl::vector_out
void vector_out(std::ostream &os, size_t n, const vec_t &v, bool endline=false)
Output the first n elements of a vector to a stream, os.
Definition: vector.h:3145
o2scl::tensor::swap_data
void swap_data(vec_t &dat)
Swap the data vector.
Definition: tensor.h:399
o2scl::ix_sum
index_spec ix_sum(size_t ix)
Sum over index ix.
o2scl::szttos
std::string szttos(size_t x)
Convert a size_t to a string.
o2scl::index_spec::index_spec
index_spec()
Default constructor.
Definition: tensor.h:98
o2scl::index_spec::empty
static const size_t empty
Empty specification.
Definition: tensor.h:76
o2scl::tensor::max
void max(vec_size_t &index, data_t &val)
Compute the index and value of the maximum value in the tensor and return the maximum.
Definition: tensor.h:709
o2scl::tensor::total_size
size_t total_size() const
Returns the size of the tensor (the product of the sizes over every index)
Definition: tensor.h:602
o2scl::tensor4::get
const data_t & get(size_t ix1, size_t ix2, size_t ix3, size_t ix4) const
Get the element indexed by (ix1,ix2,ix3,ix4)
Definition: tensor.h:1433
o2scl::ix_fixed
index_spec ix_fixed(size_t ix, size_t ix2)
Fix index ix to value ix2.
o2scl::tensor3::is_valid
void is_valid() const
Check that the o2scl::tensor3 object is valid.
Definition: tensor.h:1337
o2scl::table3d
A data structure containing one or more slices of two-dimensional data points defined on a grid.
Definition: table3d.h:78
o2scl::tensor::set_all
void set_all(data_t x)
Set all elements in a tensor to some fixed value.
Definition: tensor.h:392
o2scl::tensor::minmax_value
void minmax_value(data_t &min, data_t &max)
Compute the minimum and maximum values in the tensor.
Definition: tensor.h:718
o2scl::tensor::unpack_index
void unpack_index(size_t ix, size_vec_t &index)
Unpack the single vector index into indices.
Definition: tensor.h:641
o2scl::tensor1::set
void set(size_t index, data_t val)
Set the element indexed by index to value val.
Definition: tensor.h:1200

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