tensor.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, 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 Tensor class with arbitrary dimensions
53 
54  The elements of a tensor are typically specified as a list of
55  <tt>size_t</tt> numbers with length equal to the tensor rank.
56  For a rank-4 tensor named \c t, the element
57  <tt>t[1][2][0][3]</tt> can be obtained with something similar to
58  \code
59  size_t ix[4]={1,2,0,3};
60  double x=t.get(ix);
61  \endcode
62 
63  Empty tensors have zero rank.
64 
65  The type <tt>vec_t</tt> can be any vector type with
66  <tt>operator[]</tt>, <tt>size()</tt> and <tt>resize()</tt>
67  methods. The type <tt>vec_size_t</tt> can be any integer-like
68  vector type with <tt>operator[]</tt>, <tt>size()</tt> and
69  <tt>resize()</tt> methods.
70 
71  For I/O with tensors, see \ref o2scl_hdf::hdf_file::setd_ten()
72  and \ref o2scl_hdf::hdf_file::getd_ten() . See also
73  the discussion in the sections \ref tensor_subsect and
74  \ref vec_io_cont_subsect of the user's guide.
75 
76  The storage pattern is a generalization of row-major order.
77  In the case of a 4-rank tensor, the location of a generic
78  element is
79  \f[
80  \left( \left( i_0 s_1 + i_1 \right) s_2 + i_2 \right) s_3 + i_3 \, .
81  \f]
82  In this case the distance between two elements \f$(i_0,i_1,
83  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
84  \f$, the distance between two elements \f$(i_0,i_1,i_2, i_3)\f$
85  and \f$ (i_0,i_1+1,i_2,i_3) \f$ is \f$ s_2 s_3 \f$, and the
86  elements \f$(i_0,i_1,i_2,i_3)\f$ and \f$ (i_0,i_1,i_2,i_3+1) \f$
87  are adjacent.
88 
89  \note Slices of tensors are subsets obtained from fixing the
90  index of several dimensions, while letting other dimensions
91  vary. For a slice with one dimension not fixed, see \ref
92  vector_slice(). The \ref o2scl::tensor::vector_slice() function
93  should clearly work for uBlas vectors, and seems to work with
94  std::vector objects also, but latter use has not been fully
95  tested.
96 
97  \future Create an operator[] for tensor and not just tensor1?
98 
99  \future Could implement arithmetic operators + and - and some
100  different products.
101 
102  \future Implement copies to and from vector
103  and matrices
104 
105  \future Implement tensor contractions, i.e. tensor
106  = tensor * tensor
107 
108  \future Could be interesting to write an iterator for this class.
109 
110  */
111  template<class data_t=double, class vec_t=std::vector<data_t>,
112  class vec_size_t=std::vector<size_t> > class tensor {
113 
114  public:
115 
116 #ifndef DOXYGEN_INTERNAL
117 
118  protected:
119 
120  /// The data
121  vec_t data;
122 
123  /// A rank-sized vector of the sizes of each dimension
124  vec_size_t size;
125 
126  /// Rank
127  size_t rk;
128 
129 #endif
130 
131  public:
132 
133  /// Create an empty tensor with zero rank
134  tensor() {
135  rk=0;
136  }
137 
138  /** \brief Create a tensor of rank \c rank with sizes given in \c dim
139 
140  The parameter \c dim must be a pointer to a vector of sizes with
141  length \c rank. If the user requests any of the sizes to be
142  zero, this constructor will call the error handler, create an
143  empty tensor, and will allocate no memory.
144  */
145  template<class size_vec_t>
146  tensor(size_t rank, const size_vec_t &dim) {
147  if (rank==0) {
148  rk=0;
149  } else {
150  rk=rank;
151  for(size_t i=0;i<rk;i++) {
152  if (dim[i]==0) {
153  rk=0;
154  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
155  "rank for index "+szttos(i)+
156  " in tensor::tensor(size_t,size_vec_t)").c_str(),
157  exc_einval);
158  }
159  }
160  size.resize(rk);
161  size_t tot=1;
162  for(size_t i=0;i<rk;i++) {
163  size[i]=dim[i];
164  tot*=size[i];
165  }
166  data.resize(tot);
167  }
168  }
169 
170  ~tensor() {
171  }
172 
173  /// \name Method to check for valid object
174  //@{
175  /** \brief Check that the \ref o2scl::tensor object is valid
176  */
177  void is_valid() const {
178  if (rk==0) {
179  if (data.size()!=0) {
180  O2SCL_ERR2("Rank is zero but the data vector has non-zero size ",
181  "in tensor::is_valid().",o2scl::exc_esanity);
182  }
183  }
184 
185  if (rk!=size.size()) {
186  O2SCL_ERR2("Rank does not match size vector size ",
187  "in tensor::is_valid().",o2scl::exc_esanity);
188  }
189 
190  if (rk>0) {
191  size_t tot=1;
192  for(size_t i=0;i<rk;i++) tot*=size[i];
193  if (tot==0) {
194  O2SCL_ERR2("One entry in the size vector is zero ",
195  "in tensor::is_valid().",o2scl::exc_esanity);
196  }
197  if (tot!=data.size()) {
198  O2SCL_ERR2("Product of size vector entries does not match data ",
199  "vector size in tensor::is_valid().",o2scl::exc_esanity);
200  }
201  }
202 
203  return;
204  }
205  //@}
206 
207  /// \name Copy constructors
208  //@{
209  /** \brief Copy using <tt>operator()</tt>
210  */
213  rk=t.rk;
214  data=t.data;
215  size=t.size;
216  }
217 
218  /** \brief Copy using <tt>operator=()</tt>
219  */
222  if (this!=&t) {
223  rk=t.rk;
224  data=t.data;
225  size=t.size;
226  }
227  return *this;
228  }
229  //@}
230 
231  /// \name Clear method
232  //@{
233  /// Clear the tensor of all data and free allocated memory
234  void clear() {
235  rk=0;
236  data.resize(0);
237  size.resize(0);
238  return;
239  }
240  //@}
241 
242  /// \name Set functions
243  //@{
244  /// Set the element indexed by \c index to value \c val
245  template<class size_vec_t>
246  void set(const size_vec_t &index, data_t val) {
247 #if O2SCL_NO_RANGE_CHECK
248 #else
249  if (rk==0) {
250  O2SCL_ERR("Empty tensor in tensor::set().",exc_einval);
251  }
252  if (index[0]>=size[0]) {
253  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
254  " greater than or equal to size[0]="+
255  szttos(size[0])+" in tensor::set().").c_str(),
256  exc_eindex);
257  }
258 #endif
259  size_t ix=index[0];
260  for(size_t i=1;i<rk;i++) {
261 #if O2SCL_NO_RANGE_CHECK
262 #else
263  if (index[i]>=size[i]) {
264  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
265  szttos(index[i])+" greater than or equal to size "+
266  szttos(size[i])+" in tensor::set().").c_str(),
267  exc_eindex);
268  }
269 #endif
270  ix*=size[i];
271  ix+=index[i];
272  }
273  data[ix]=val;
274  return;
275  }
276 
277  /// Set all elements in a tensor to some fixed value
278  void set_all(data_t x) {
279  for(size_t i=0;i<total_size();i++) data[i]=x;
280  return;
281  }
282 
283  /** \brief Swap the data vector
284  */
285  void swap_data(vec_t &dat) {
286  if (data.size()!=dat.size()) {
287  O2SCL_ERR2("Size of new vector does not equal tensor size in ",
288  "tensor::swap_data().",o2scl::exc_einval);
289  }
290  std::swap(dat,data);
291  return;
292  }
293  //@}
294 
295  /// \name Get functions
296  //@{
297  /// Get the element indexed by \c index
298  template<class size_vec_t> data_t &get(const size_vec_t &index) {
299 #if O2SCL_NO_RANGE_CHECK
300 #else
301  if (rk==0) {
302  O2SCL_ERR("Empty tensor in tensor::get().",exc_einval);
303  }
304  if (index[0]>=size[0]) {
305  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
306  " greater than or equal to size[0]="+
307  szttos(size[0])+" in tensor::get().").c_str(),
308  exc_eindex);
309  }
310 #endif
311  size_t ix=index[0];
312  for(size_t i=1;i<rk;i++) {
313 #if O2SCL_NO_RANGE_CHECK
314 #else
315  if (index[i]>=size[i]) {
316  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
317  szttos(index[i])+" greater than or equal to size "+
318  szttos(size[i])+" in tensor::get().").c_str(),
319  exc_eindex);
320  }
321 #endif
322  ix*=size[i];
323  ix+=index[i];
324  }
325  return data[ix];
326  }
327 
328  /// Get a const reference to the element indexed by \c index
329  template<class size_vec_t>
330  data_t const &get(const size_vec_t &index) const {
331 #if O2SCL_NO_RANGE_CHECK
332 #else
333  if (rk==0) {
334  O2SCL_ERR("Empty tensor in tensor::get() (const).",exc_einval);
335  }
336  if (index[0]>=size[0]) {
337  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
338  " greater than or equal to size[0]="+
339  szttos(size[0])+" in tensor::get() (const).").c_str(),
340  exc_eindex);
341  }
342 #endif
343  size_t ix=index[0];
344  for(size_t i=1;i<rk;i++) {
345 #if O2SCL_NO_RANGE_CHECK
346 #else
347  if (index[i]>=size[i]) {
348  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
349  szttos(index[i])+" greater than or equal to size "+
350  szttos(size[i])+" in tensor::get() (const).").c_str(),
351  exc_eindex);
352  }
353 #endif
354  ix*=size[i];
355  ix+=index[i];
356  }
357  return data[ix];
358  }
359  //@}
360 
361  typedef boost::numeric::ublas::vector_slice<
362  boost::numeric::ublas::vector<data_t> > ubvector_slice;
363  typedef boost::numeric::ublas::slice slice;
364 
365  /// \name Slice function
366  //@{
367  /** \brief Fix all but one index to create a vector
368 
369  This fixes all of the indices to the values given in \c index
370  except for the index number \c ix, and returns the
371  corresponding vector, whose length is equal to the size of the
372  tensor in that index. The value <tt>index[ix]</tt> is ignored.
373 
374  For example, for a rank 3 tensor allocated with
375  \code
376  tensor t;
377  size_t dim[3]={3,4,5};
378  t.resize(3,dim);
379  \endcode
380  the following code
381  \code
382  size_t index[3]={1,0,3};
383  ubvector_view v=t.vector_slice(1,index);
384  \endcode
385  Gives a vector \c v of length 4 which refers to the values
386  <tt>t(1,0,3)</tt>, <tt>t(1,1,3)</tt>, <tt>t(1,2,3)</tt>, and
387  <tt>t(1,3,3)</tt>.
388  */
389  template<class size_vec_t>
390  ubvector_slice vector_slice(size_t ix, const size_vec_t &index) {
391  if (ix+1>rk) {
392  O2SCL_ERR((((std::string)"Specified index ")+szttos(ix)+
393  " greater than or equal to rank "+szttos(rk)+
394  " in tensor::vector_slice()").c_str(),
395  exc_eindex);
396  }
397  size_t start;
398  if (ix==0) start=0;
399  else start=index[0];
400  for(size_t i=1;i<rk;i++) {
401  start*=size[i];
402  if (i!=ix) start+=index[i];
403  }
404  size_t stride=1;
405  for(size_t i=ix+1;i<rk;i++) stride*=size[i];
406  return ubvector_slice(data,slice(start,stride,size[ix]));
407  }
408  //@}
409 
410 #ifdef O2SCL_NEVER_DEFINED
411 
412  /** \brief Fix all but two indices to create a matrix
413 
414  One could use ublas::make_matrix_from_pointer() and
415  then ublas matrix slicing to perform this operation, but
416  only if the template types are built on ublas objects.
417  */
418  template<class size_vec_t> ubmatrix_array matrix_slice() {
419  }
420 
421 #endif
422 
423  /// \name Resize method
424  //@{
425  /** \brief Resize the tensor to rank \c rank with sizes
426  given in \c dim
427 
428  The parameter \c dim must be a vector of sizes with a length
429  equal to \c rank. This resize method is always destructive.
430 
431  If the user requests any of the sizes to be zero, this
432  function will call the error handler.
433  */
434  template<class size_vec_t>
435  void resize(size_t rank, const size_vec_t &dim) {
436  for(size_t i=0;i<rank;i++) {
437  if (dim[i]==0) {
438  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
439  "rank for index "+szttos(i)+" in tensor::"+
440  "resize().").c_str(),exc_einval);
441  }
442  }
443  rk=rank;
444  size.resize(rk);
445  if (rk==0) {
446  data.resize(0);
447  } else {
448  size_t tot=1;
449  for(size_t i=0;i<rk;i++) {
450  size[i]=dim[i];
451  tot*=size[i];
452  }
453  data.resize(tot);
454  }
455  return;
456  }
457  //@}
458 
459  /// \name Size functions
460  //@{
461  /// Return the rank of the tensor
462  size_t get_rank() const { return rk; }
463 
464  /// Returns the size of the ith index
465  size_t get_size(size_t i) const {
466  if (i>=rk) {
467  O2SCL_ERR((((std::string)"Specified index ")+szttos(i)+
468  " greater than or equal to rank "+szttos(rk)+
469  " in tensor::get_size()").c_str(),
470  exc_einval);
471  }
472  return size[i];
473  }
474 
475  /// Return the full vector of sizes
476  const vec_size_t &get_size_arr() const {
477  return size;
478  }
479 
480  /// Return the full data vector
481  const vec_t &get_data() const {
482  return data;
483  }
484 
485  /** \brief Returns the size of the tensor (the product of
486  the sizes over every index)
487  */
488  size_t total_size() const {
489  if (rk==0) return 0;
490  size_t tot=1;
491  for(size_t i=0;i<rk;i++) tot*=size[i];
492  return tot;
493  }
494  //@}
495 
496  /// \name Index manipulation
497  //@{
498  /// Pack the indices into a single vector index
499  template<class size_vec_t>
500  size_t pack_indices(const size_vec_t &index) {
501  if (rk==0) {
502  O2SCL_ERR("Empty tensor in tensor::pack_indices().",exc_einval);
503  return 0;
504  }
505  if (index[0]>=size[0]) {
506  O2SCL_ERR((((std::string)"Value of index[0]=")+szttos(index[0])+
507  " greater than or equal to size[0]="+
508  szttos(size[0])+" in tensor::pack_indices().").c_str(),
509  exc_eindex);
510  }
511  size_t ix=index[0];
512  for(size_t i=1;i<rk;i++) {
513  if (index[i]>=size[i]) {
514  O2SCL_ERR((((std::string)"Value of index[")+szttos(i)+"]="+
515  szttos(index[i])+" greater than or equal to size "+
516  szttos(size[i])+" in tensor::pack_indices().").c_str(),
517  exc_eindex);
518  }
519  ix*=size[i];
520  ix+=index[i];
521  }
522  return ix;
523  }
524 
525  /// Unpack the single vector index into indices
526  template<class size_vec_t>
527  void unpack_index(size_t ix, size_vec_t &index) {
528  if (ix>=total_size()) {
529  O2SCL_ERR((((std::string)"Value of index ")+szttos(ix)+
530  " greater than or equal to total size"+
531  szttos(total_size())+
532  " in tensor::unpack_index().").c_str(),
533  exc_eindex);
534  return;
535  }
536  size_t ix2, sub_size;
537  for(size_t i=0;i<rk;i++) {
538  if (i==rk-1) {
539  index[i]=ix;
540  } else {
541  sub_size=1;
542  for(size_t j=i+1;j<rk;j++) sub_size*=size[j];
543  index[i]=ix/sub_size;
544  // (Remember we're doing integer arithmetic here.)
545  ix-=sub_size*(ix/sub_size);
546  }
547  }
548  return;
549  }
550  //@}
551 
552  /// \name Minimum, maximum, and sum
553  //@{
554  /** \brief Compute the minimum value in the tensor
555  */
556  data_t min_value() {
557  return o2scl::vector_min_value<vec_t,data_t>(total_size(),data);
558  }
559 
560  /** \brief Compute the index of the minimum value in the tensor
561  */
562  void min_index(vec_size_t &index) {
563  size_t ix=o2scl::vector_min_index<vec_t,data_t>(total_size(),data);
564  unpack_index(ix,index);
565  return;
566  }
567 
568  /** \brief Compute the index of the minimum value in the tensor
569  and return the minimum
570  */
571  void min(vec_size_t &index, data_t &val) {
572  size_t ix;
573  o2scl::vector_min<vec_t,data_t>(total_size(),data,ix,val);
574  unpack_index(ix,index);
575  return;
576  }
577 
578  /** \brief Compute the maximum value in the tensor
579  */
580  data_t max_value() {
581  return o2scl::vector_max_value<vec_t,data_t>(total_size(),data);
582  }
583 
584  /** \brief Compute the index of the maximum value in the tensor
585  */
586  void max_index(vec_size_t &index) {
587  size_t ix=o2scl::vector_max_index<vec_t,data_t>(total_size(),data);
588  unpack_index(ix,index);
589  return;
590  }
591 
592  /** \brief Compute the index and value of the maximum value in the tensor
593  and return the maximum
594  */
595  void max(vec_size_t &index, data_t &val) {
596  size_t ix;
597  o2scl::vector_max<vec_t,data_t>(total_size(),data,ix,val);
598  unpack_index(ix,index);
599  return;
600  }
601 
602  /** \brief Compute the minimum and maximum values in the tensor
603  */
604  void minmax_value(data_t &min, data_t &max) {
605  return o2scl::vector_minmax_value<vec_t,data_t>(total_size(),data,min,max);
606  }
607 
608  /** \brief Compute the indices of the minimum and maximum values in the tensor
609  */
610  void minmax_index(vec_size_t &index_min, vec_size_t &index_max) {
611  size_t ix_min, ix_max;
612  o2scl::vector_minmax_index<vec_t,data_t>
613  (total_size(),data,ix_min,ix_max);
614  unpack_index(ix_min,index_min);
615  unpack_index(ix_max,index_max);
616  return;
617  }
618 
619  /** \brief Compute the indices and values of the maximum and minimum
620  in the tensor
621  */
622  void minmax(vec_size_t &index, size_t &index_min, data_t &min,
623  size_t &index_max, data_t &max) {
624  size_t ix_min, ix_max;
625  o2scl::vector_minmax<vec_t,data_t>(total_size(),data,ix_min,min,
626  ix_max,max);
627  unpack_index(ix_min,index_min);
628  unpack_index(ix_max,index_max);
629  return;
630  }
631 
632  /** \brief Return the sum over every element in the tensor
633  */
634  double total_sum() const {
635  if (rk==0) return 0.0;
636  double tot=0.0;
637  for(size_t i=0;i<data.size();i++) {
638  tot+=data[i];
639  }
640  return tot;
641  }
642  //@}
643 
644  /// \name Slicing and converting to table3d objects
645  //@{
646  /** \brief Convert to a \ref o2scl::table3d object by
647  summing over all but two indices
648  */
650  (size_t ix_x, size_t ix_y, table3d &tab, std::string x_name="x",
651  std::string y_name="y", std::string slice_name="z") {
652 
653  // Get current table3d grid
654  size_t nx, ny;
655  tab.get_size(nx,ny);
656 
657  if (nx==0 && ny==0) {
658 
659  if (x_name.length()==0) x_name="x";
660  if (y_name.length()==0) y_name="y";
661 
662  // If there's no grid, then create a grid in the table3d
663  // object which just enumerates the indices
664  std::vector<double> grid_x(size[ix_x]), grid_y(size[ix_y]);
665  for(size_t i=0;i<size[ix_x];i++) {
666  grid_x[i]=((double)i);
667  }
668  for(size_t i=0;i<size[ix_y];i++) {
669  grid_y[i]=((double)i);
670  }
671  tab.set_xy("x",grid_x.size(),grid_x,
672  "y",grid_y.size(),grid_y);
673  // Now that the grid is set, get nx and ny
674  tab.get_size(nx,ny);
675  }
676 
677  // Check that the grids are commensurate
678  if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
679  O2SCL_ERR2("Grids not commensurate in ",
680  "tensor_grid::convert_table3d_sum().",exc_einval);
681  }
682 
683  tab.set_slice_all(slice_name,0.0);
684 
685  std::vector<size_t> ix;
686  for(size_t i=0;i<this->total_size();i++) {
687  this->unpack_index(i,ix);
688  tab.set(ix[ix_x],ix[ix_y],slice_name,
689  tab.get(ix[ix_x],ix[ix_y],slice_name)+
690  this->data[i]);
691  }
692 
693  return;
694  }
695  //@}
696 
697  /// \name Slicing and summing to create tensor_grid objects
698  //@{
699 #ifdef O2SCL_NEVER_DEFINED
700  // commenting this function out because clang
701  // had a problem with ix_new=ix_old[mapping[j]];
702 
703 
704  /** \brief Sum over some (but not all) tensor indices to obtain
705  a smaller rank tensor
706  */
707  template<class size_vec2_t, class vec2_t>
708  tensor<> sum_slice(size_vec2_t &isum) {
709 
710  if (isum.size()==0) {
711  O2SCL_ERR2("Specified empty vector in ",
712  "tensor::copy_slice_interp().",
714  }
715  if (this->rk<1+isum.size()) {
716  O2SCL_ERR2("Summed too many indices in ",
717  "tensor::copy_slice_interp().",
719  }
720 
721  // Determine new rank
722  size_t rank_new=this->rk-isum.size();
723 
724  // Determine the mapping between indices and the new size array
725  std::vector<size_t> mapping;
726  std::vector<size_t> sz_new;
727 
728  for(size_t i=0;i<this->rk;i++) {
729  bool found=false;
730  for(size_t j=0;j<isum.size();j++) {
731  if (isum[j]==i) found=true;
732  }
733  if (found==false) {
734  sz_new.push_back(this->get_size(i));
735  mapping.push_back(i);
736  }
737  }
738 
739  // Create the new tensor object
740  tensor<> tg_new(rank_new,sz_new);
741 
742  // Set all to zero
743  tg_new.set_all(0.0);
744 
745  // Loop over the old tensor object
746  for(size_t i=0;i<total_size();i++) {
747 
748  std::vector<size_t> ix_old(this->rk), ix_new(rank_new);
749  this->unpack_index(i,ix_old);
750  for(size_t j=0;j<rank_new;j++) {
751  ix_new=ix_old[mapping[j]];
752  }
753 
754  size_t k=tg_new.pack_indices(ix_new);
755 
756  tg_new.set(ix_new,tg_new.get(ix_new)+this->data[k]);
757  }
758 
759  return tg_new;
760  }
761 #endif
762 
763 #ifdef O2SCL_NEVER_DEFINED
764 
765  class index_spec {
766 
767  public:
768 
769  size_t type;
770  size_t ix1;
771  size_t ix2;
772  double val;
773 
774  static const size_t index=1;
775  static const size_t fixed=2;
776  static const size_t sum=3;
777  static const size_t contract=4;
778  static const size_t reverse=5;
779  static const size_t range=6;
780 
781  // When a grid is available
782  static const size_t interp=7;
783  static const size_t grid=8;
784 
785  index_spec(size_t typ, size_t i1, size_t i2, double v) {
786  type=typ;
787  ix1=i1;
788  ix2=i2;
789  val=v;
790  }
791 
792  };
793 
794  index_spec &ix_index(size_t ix) {
795  return index_spec(index_spec::index,ix,0,0.0);
796  }
797 
798  index_spec &ix_fixed(size_t ix, size_t ix2) {
799  return index_spec(index_spec::fixed,ix,ix2,0.0);
800  }
801 
802  index_spec &ix_sum(size_t ix) {
803  return index_spec(index_spec::sum,ix,0,0.0);
804  }
805 
806  index_spec &ix_contract(size_t ix, size_t ix2) {
807  return index_spec(index_spec::contract,ix,ix2,0.0);
808  }
809 
810  index_spec &ix_reverse(size_t ix) {
811  return index_spec(index_spec::reverse,ix,0,0.0);
812  }
813 
814  index_spec &ix_interp(size_t ix, double v) {
815  return index_spec(index_spec::interp,ix,0,v);
816  }
817 
818  index_spec &ix_range(size_t ix, double v) {
819  return index_spec(index_spec::range,ix,0,v);
820  }
821 
822  index_spec &ix_grid(size_t ix, double v) {
823  return index_spec(index_spec::grid,ix,0,v);
824  }
825 
826  // AWS: I'm waiting on this function because I would like
827  // to generalize it by creating a new index_spec class
828  // which allows for more general contractions, sums, and
829  // index rearrangements, etc.
830 
831  /** \brief Create a copy by contracting, summing over, fixing,
832  and rearranging tensor indices
833  */
834  tensor<> copy(std::vector<index_spec> &spec) {
835 
836  size_t rank_old=this->rk;
837  size_t rank_new=0;
838  permutation p(rank_old);
839  std::vector<size_t> size_new;
840  size_t ip=0;
841  size_t n_inner_loop=1;
842 
843  for(size_t i=0;i<spec.size();i++) {
844  if (spec[i].type==index_spec::index ||
845  spec[i].type==index_spec::reverse) {
846  rank_new++;
847  p[ip]=spec[i].ix1;
848  ip++;
849  size_new.push_back(this->size[spec[i].ix1]);
850  } else if (spec_type==index_spec::contract) {
851  p[ip]=spec[i].ix1;
852  ip++;
853  p[ip]=spec[i].ix2;
854  ip++;
855  if (size[spec[i].ix1]<size[spec[i].ix2]) {
856  n_inner_loop*=size[spec[i].ix1];
857  } else {
858  n_inner_loop*=size[spec[i].ix2];
859  }
860  } else if (spec[i].type==index_spec::sum) {
861  p[ip]=spec[i].ix1;
862  ip++;
863  n_inner_loop*=size[spec[i].ix1];
864  } else if (spec[i].type==index_spec::fixed) {
865  p[ip]=spec[i].ix1;
866  ip++;
867  }
868  }
869 
870  if (rank_new==0) {
871  O2SCL_ERR("Zero new indices in tensor::copy().",
873  }
874  if (p.valid()==false) {
875  O2SCL_ERR("Not all indices accounted for in tensor::copy().",
877  }
878 
879  // Create the new tensor_grid object and set the new grid
880  tensor<> t_new(rank_new,size_new);
881 
882  // Index arrays
883  std::vector<size_t> ix_new(rank_new);
884  std::vector<size_t> ix_old(rank_old);
885 
886  // Loop over the new tensor object
887  for(size_t i=0;i<t_new.total_size();i++) {
888 
889  // Find the location in the new tensor_grid object
890  t_new.unpack_index(i,ix_new);
891 
892  double val=0.0;
893 
894  for(size_t j=0;j<n_inner_loop;j++) {
895  }
896 
897  // Set the new point by performing the linear interpolation
898  tg_new.set(ix_new,val);
899  }
900 
901  return tg_new;
902  }
903  //@}
904 #endif
905 
906  };
907 
908  /** \brief Rank 1 tensor
909  */
910  template<class data_t=double, class vec_t=std::vector<data_t>,
911  class vec_size_t=std::vector<size_t> > class tensor1 :
912  public tensor<data_t,vec_t,vec_size_t> {
913 
914  public:
915 
916  /// Create an empty tensor
917  tensor1() : tensor<data_t,vec_t,vec_size_t>() {}
918 
919  /// Create a rank 1 tensory of size \c sz
920  tensor1(size_t sz) : tensor<data_t,vec_t,vec_size_t>() {
921  vec_size_t sizex(1);
922  sizex[0]=sz;
923  this->resize(1,sizex);
924  }
925 
926  /// \name Method to check for valid object
927  //@{
928  /** \brief Check that the \ref o2scl::tensor1 object is valid
929  */
930  void is_valid() const {
932  if (this->rk>1) {
933  O2SCL_ERR2("Rank is neither 0 nor 1 in ",
934  "tensor1::is_valid().",
936 
937  }
938  return;
939  }
940  //@}
941 
942  /// \name Specialized get and set functions
943  //@{
944  /// Get the element indexed by \c ix
945  data_t &get(size_t ix) {
947  }
948 
949  /// Get the element indexed by \c ix
950  const data_t &get(size_t ix) const {
952  }
953 
954  /// Set the element indexed by \c index to value \c val
955  void set(size_t index, data_t val)
956  { tensor<data_t,vec_t,vec_size_t>::set(&index,val); }
957 
958  /** \brief Set the element indexed by \c index to value \c val
959 
960  (We have to explicitly provide this version since the set()
961  function is overloaded in this child of \ref tensor.)
962  */
963  template<class size_vec_t>
964  void set(const size_vec_t &index, data_t val) {
966  }
967  //@}
968 
969  /// \name Specialized operator functions
970  //@{
971  /// Get an element using array-like indexing
972  data_t &operator[](size_t ix) { return this->data[ix]; }
973 
974  /// Get an element using array-like indexing (const version)
975  const data_t &operator[](size_t ix) const { return this->data[ix]; }
976 
977  /// Get an element using operator()
978  data_t &operator()(size_t ix) { return this->data[ix]; }
979 
980  /// Get an element using operator() (const version)
981  const data_t &operator()(size_t ix) const { return this->data[ix]; }
982  //@}
983 
984  };
985 
986  /** \brief Rank 2 tensor
987  */
988  template<class data_t=double, class vec_t=std::vector<data_t>,
989  class vec_size_t=std::vector<size_t> > class tensor2 :
990  public tensor<data_t,vec_t,vec_size_t> {
991 
992  public:
993 
994  /// Create an empty tensor
995  tensor2() : tensor<data_t,vec_t,vec_size_t>() {}
996 
997  /// Create a rank 2 tensor of size \c (sz,sz2)
998  tensor2(size_t sz, size_t sz2) : tensor<data_t,vec_t,vec_size_t>() {
999  this->rk=2;
1000  this->size.resize(2);
1001  this->size[0]=sz;
1002  this->size[1]=sz2;
1003  size_t tot=sz*sz2;
1004  this->data.resize(tot);
1005  }
1006 
1007  /// \name Method to check for valid object
1008  //@{
1009  /** \brief Check that the \ref o2scl::tensor2 object is valid
1010  */
1011  void is_valid() const {
1013  if (this->rk!=0 && this->rk!=2) {
1014  O2SCL_ERR2("Rank is neither 0 nor 2 in ",
1015  "tensor2::is_valid().",
1017 
1018  }
1019  return;
1020  }
1021  //@}
1022 
1023  /// \name Specialized get and set functions
1024  //@{
1025  /// Get the element indexed by \c (ix1,ix2)
1026  data_t &get(size_t ix1, size_t ix2) {
1027  size_t sz[2]={ix1,ix2};
1029  }
1030 
1031  /// Get the element indexed by \c (ix1,ix2)
1032  const data_t &get(size_t ix1, size_t ix2) const {
1033  size_t sz[2]={ix1,ix2};
1035  }
1036 
1037  /// Set the element indexed by \c (ix1,ix2) to value \c val
1038  void set(size_t ix1, size_t ix2, data_t val) {
1039  size_t sz[2]={ix1,ix2};
1041  return;
1042  }
1043 
1044  /** \brief Set the element indexed by \c index to value \c val
1045 
1046  (We have to explicitly provide this version since the set()
1047  function is overloaded in this child of \ref tensor.)
1048  */
1049  template<class size_vec_t>
1050  void set(const size_vec_t &index, data_t val) {
1052  return;
1053  }
1054 
1055  /// Get the element indexed by \c (ix1,ix2)
1056  data_t &operator()(size_t ix, size_t iy)
1057  { return this->data[ix*this->size[1]+iy]; }
1058 
1059  /// Get the element indexed by \c (ix1,ix2) (const version)
1060  const data_t &operator()(size_t ix, size_t iy) const
1061  { return this->data[ix*this->size[1]+iy]; }
1062  //@}
1063  };
1064 
1065  /** \brief Rank 3 tensor
1066  */
1067  template<class data_t=double, class vec_t=std::vector<data_t>,
1068  class vec_size_t=std::vector<size_t> > class tensor3 :
1069  public tensor<data_t,vec_t,vec_size_t> {
1070 
1071  public:
1072 
1073  /// Create an empty tensor
1074  tensor3() : tensor<data_t,vec_t,vec_size_t>() {}
1075 
1076  /// Create a rank 3 tensor of size \c (sz,sz2,sz3)
1077  tensor3(size_t sz, size_t sz2, size_t sz3) :
1078  tensor<data_t,vec_t,vec_size_t>() {
1079  this->rk=3;
1080  this->size.resize(3);
1081  this->size[0]=sz;
1082  this->size[1]=sz2;
1083  this->size[2]=sz3;
1084  size_t tot=sz*sz2*sz3;
1085  this->data.resize(tot);
1086  }
1087 
1088  /// \name Method to check for valid object
1089  //@{
1090  /** \brief Check that the \ref o2scl::tensor3 object is valid
1091  */
1092  void is_valid() const {
1094  if (this->rk!=0 && this->rk!=3) {
1095  O2SCL_ERR2("Rank is neither 0 nor 3 in ",
1096  "tensor3::is_valid().",
1098 
1099  }
1100  return;
1101  }
1102  //@}
1103 
1104  /// \name Specialized get and set functions
1105  //@{
1106  /// Get the element indexed by \c (ix1,ix2,ix3)
1107  data_t &get(size_t ix1, size_t ix2, size_t ix3) {
1108  size_t sz[3]={ix1,ix2,ix3};
1110  }
1111 
1112  /// Get the element indexed by \c (ix1,ix2,ix3)
1113  const data_t &get(size_t ix1, size_t ix2, size_t ix3) const {
1114  size_t sz[3]={ix1,ix2,ix3};
1116  }
1117 
1118  /// Set the element indexed by \c (ix1,ix2,ix3) to value \c val
1119  void set(size_t ix1, size_t ix2, size_t ix3, data_t val) {
1120  size_t sz[3]={ix1,ix2,ix3};
1122  return;
1123  }
1124 
1125  /** \brief Set the element indexed by \c index to value \c val
1126 
1127  (We have to explicitly provide this version since the set()
1128  function is overloaded in this child of \ref tensor.)
1129  */
1130  template<class size_vec_t>
1131  void set(const size_vec_t &index, data_t val) {
1133  return;
1134  }
1135  //@}
1136 
1137  };
1138 
1139  /** \brief Rank 4 tensor
1140  */
1141  template<class data_t=double, class vec_t=std::vector<data_t>,
1142  class vec_size_t=std::vector<size_t> > class tensor4 :
1143  public tensor<data_t,vec_t,vec_size_t> {
1144 
1145  public:
1146 
1147  /// Create an empty tensor
1148  tensor4() : tensor<data_t,vec_t,vec_size_t>() {}
1149 
1150  /// Create a rank 4 tensor of size \c (sz,sz2,sz3,sz4)
1151  tensor4(size_t sz, size_t sz2, size_t sz3, size_t sz4) :
1152  tensor<data_t,vec_t,vec_size_t>() {
1153  this->rk=4;
1154  this->size.resize(4);
1155  this->size[0]=sz;
1156  this->size[1]=sz2;
1157  this->size[2]=sz3;
1158  this->size[3]=sz4;
1159  size_t tot=sz*sz2*sz3*sz4;
1160  this->data.resize(tot);
1161  }
1162 
1163  /// \name Method to check for valid object
1164  //@{
1165  /** \brief Check that the \ref o2scl::tensor4 object is valid
1166  */
1167  void is_valid() const {
1169  if (this->rk!=0 && this->rk!=4) {
1170  O2SCL_ERR2("Rank is neither 0 nor 4 in ",
1171  "tensor4::is_valid().",
1173 
1174  }
1175  return;
1176  }
1177  //@}
1178 
1179  /// \name Specialized get and set functions
1180  //@{
1181  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1182  data_t &get(size_t ix1, size_t ix2, size_t ix3, size_t ix4) {
1183  size_t sz[4]={ix1,ix2,ix3,ix4};
1185  }
1186 
1187  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1188  const data_t &get(size_t ix1, size_t ix2, size_t ix3,
1189  size_t ix4) const {
1190  size_t sz[4]={ix1,ix2,ix3,ix4};
1192  }
1193 
1194  /// Set the element indexed by \c (ix1,ix2,ix3,ix4) to value \c val
1195  void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4,
1196  data_t val) {
1197  size_t sz[4]={ix1,ix2,ix3,ix4};
1199  return;
1200  }
1201 
1202  /** \brief Set the element indexed by \c index to value \c val
1203 
1204  (We have to explicitly provide this version since the set()
1205  function is overloaded in this child of \ref tensor.)
1206  */
1207  template<class size_vec_t>
1208  void set(const size_vec_t &index, data_t val) {
1210  return;
1211  }
1212  //@}
1213  };
1214 
1215  /** \brief Output a tensor to a stream
1216  */
1217  template<class tensor_t>
1218  void tensor_out(std::ostream &os, tensor_t &t, bool pretty=true) {
1219 
1220  if (pretty) {
1221 
1222  size_t rk=t.get_rank();
1223  os << "rank: " << rk << " sizes: ";
1224  for(size_t i=0;i<rk;i++) {
1225  os << t.get_size(i) << " ";
1226  }
1227  os << std::endl;
1228  auto &data=t.get_data();
1229  std::vector<size_t> ix(rk);
1230  std::vector<std::string> sv, sv_out;
1231  for(size_t i=0;i<t.total_size();i++) {
1232  t.unpack_index(i,ix);
1233  std::string tmp="(";
1234  for(size_t j=0;j<rk;j++) {
1235  if (j!=rk-1) {
1236  tmp+=o2scl::szttos(ix[j])+",";
1237  } else {
1238  tmp+=o2scl::szttos(ix[j]);
1239  }
1240  }
1241  tmp+="): "+o2scl::dtos(data[i]);
1242  sv.push_back(tmp);
1243  }
1244  screenify(sv.size(),sv,sv_out);
1245  for(size_t k=0;k<sv_out.size();k++) {
1246  os << sv_out[k] << std::endl;
1247  }
1248 
1249  } else {
1250 
1251  size_t rk=t.get_rank();
1252  os << rk << " ";
1253  for(size_t i=0;i<rk;i++) {
1254  os << t.get_size(i) << " ";
1255  }
1256  os << std::endl;
1257  auto &data=t.get_data();
1258  for(size_t i=0;i<t.total_size();i++) {
1259  os << data[i] << " ";
1260  if (i%10==9) os << std::endl;
1261  }
1262  os << std::endl;
1263 
1264  }
1265 
1266  return;
1267  }
1268 
1269 #ifndef DOXYGEN_NO_O2NS
1270 }
1271 #endif
1272 
1273 #endif
1274 
1275 
1276 
void resize(size_t rank, const size_vec_t &dim)
Resize the tensor to rank rank with sizes given in dim.
Definition: tensor.h:435
void unpack_index(size_t ix, size_vec_t &index)
Unpack the single vector index into indices.
Definition: tensor.h:527
data_t min_value()
Compute the minimum value in the tensor.
Definition: tensor.h:556
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:610
Rank 1 tensor.
Definition: tensor.h:911
tensor3(size_t sz, size_t sz2, size_t sz3)
Create a rank 3 tensor of size (sz,sz2,sz3)
Definition: tensor.h:1077
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:1068
tensor1(size_t sz)
Create a rank 1 tensory of size sz.
Definition: tensor.h:920
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:622
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:1151
tensor1()
Create an empty tensor.
Definition: tensor.h:917
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
double & get(size_t ix, size_t iy, std::string name)
Get element in slice name at location ix,iy
size_t get_rank() const
Return the rank of the tensor.
Definition: tensor.h:462
invalid argument supplied by user
Definition: err_hnd.h:59
const data_t & operator[](size_t ix) const
Get an element using array-like indexing (const version)
Definition: tensor.h:975
tensor2()
Create an empty tensor.
Definition: tensor.h:995
A class for representing permutations.
Definition: permutation.h:70
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor.h:234
ubvector_slice vector_slice(size_t ix, const size_vec_t &index)
Fix all but one index to create a vector.
Definition: tensor.h:390
void minmax_value(data_t &min, data_t &max)
Compute the minimum and maximum values in the tensor.
Definition: tensor.h:604
data_t max_value()
Compute the maximum value in the tensor.
Definition: tensor.h:580
void set_slice_all(std::string name, double val)
Set all of the values in slice name to val.
tensor2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
Definition: tensor.h:998
size_t total_size() const
Returns the size of the tensor (the product of the sizes over every index)
Definition: tensor.h:488
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:650
void is_valid() const
Check that the o2scl::tensor4 object is valid.
Definition: tensor.h:1167
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:595
void min_index(vec_size_t &index)
Compute the index of the minimum value in the tensor.
Definition: tensor.h:562
void get_size(size_t &nx, size_t &ny) const
Get the size of the slices.
void is_valid() const
Check that the o2scl::tensor3 object is valid.
Definition: tensor.h:1092
vec_t data
The data.
Definition: tensor.h:121
data_t & get(const size_vec_t &index)
Get the element indexed by index.
Definition: tensor.h:298
size_t rk
Rank.
Definition: tensor.h:127
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
tensor4()
Create an empty tensor.
Definition: tensor.h:1148
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
void max_index(vec_size_t &index)
Compute the index of the maximum value in the tensor.
Definition: tensor.h:586
data_t & operator()(size_t ix, size_t iy)
Get the element indexed by (ix1,ix2)
Definition: tensor.h:1056
void is_valid() const
Check that the o2scl::tensor object is valid.
Definition: tensor.h:177
const vec_size_t & get_size_arr() const
Return the full vector of sizes.
Definition: tensor.h:476
const data_t & operator()(size_t ix) const
Get an element using operator() (const version)
Definition: tensor.h:981
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:127
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
data_t & operator[](size_t ix)
Get an element using array-like indexing.
Definition: tensor.h:972
const vec_t & get_data() const
Return the full data vector.
Definition: tensor.h:481
size_t pack_indices(const size_vec_t &index)
Pack the indices into a single vector index.
Definition: tensor.h:500
void is_valid() const
Check that the o2scl::tensor2 object is valid.
Definition: tensor.h:1011
double total_sum() const
Return the sum over every element in the tensor.
Definition: tensor.h:634
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:571
tensor3()
Create an empty tensor.
Definition: tensor.h:1074
data_t & operator()(size_t ix)
Get an element using operator()
Definition: tensor.h:978
tensor(size_t rank, const size_vec_t &dim)
Create a tensor of rank rank with sizes given in dim.
Definition: tensor.h:146
const data_t & operator()(size_t ix, size_t iy) const
Get the element indexed by (ix1,ix2) (const version)
Definition: tensor.h:1060
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.
A data structure containing one or more slices of two-dimensional data points defined on a grid...
Definition: table3d.h:77
void set(const size_vec_t &index, data_t val)
Set the element indexed by index to value val.
Definition: tensor.h:246
void tensor_out(std::ostream &os, tensor_t &t, bool pretty=true)
Output a tensor to a stream.
Definition: tensor.h:1218
Rank 2 tensor.
Definition: tensor.h:989
size_t get_size(size_t i) const
Returns the size of the ith index.
Definition: tensor.h:465
Rank 4 tensor.
Definition: tensor.h:1142
void swap_data(vec_t &dat)
Swap the data vector.
Definition: tensor.h:285
Invalid index for array or matrix.
Definition: err_hnd.h:123
tensor()
Create an empty tensor with zero rank.
Definition: tensor.h:134
vec_size_t size
A rank-sized vector of the sizes of each dimension.
Definition: tensor.h:124
bool valid() const
Check to see that a permutation is valid.
Definition: permutation.h:257
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:119
std::string szttos(size_t x)
Convert a size_t to a string.
The default vector type from uBlas.
Definition: vector.h:3260
Interpolation class for general vectors.
Definition: interp.h:1654
Tensor class with arbitrary dimensions.
Definition: tensor.h:112
void is_valid() const
Check that the o2scl::tensor1 object is valid.
Definition: tensor.h:930
void set_all(data_t x)
Set all elements in a tensor to some fixed value.
Definition: tensor.h:278

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