tensor_grid.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_GRID_H
24 #define O2SCL_TENSOR_GRID_H
25 
26 /** \file tensor_grid.h
27  \brief File defining \ref o2scl::tensor_grid and rank-specific children
28 */
29 
30 #include <iostream>
31 #include <cstdlib>
32 #include <string>
33 #include <fstream>
34 #include <sstream>
35 
36 #include <gsl/gsl_matrix.h>
37 #include <gsl/gsl_ieee_utils.h>
38 
39 #include <o2scl/err_hnd.h>
40 #include <o2scl/interp.h>
41 #include <o2scl/tensor.h>
42 #include <o2scl/table3d.h>
43 
44 // Forward definition of the tensor_grid class for HDF I/O
45 namespace o2scl {
46  template<class vec_t, class vec_size_t> class tensor_grid;
47 }
48 
49 // Forward definition of HDF I/O to extend friendship
50 namespace o2scl_hdf {
51  class hdf_file;
52  template<class vec_t, class vec_size_t>
53  void hdf_input(hdf_file &hf, o2scl::tensor_grid<vec_t,vec_size_t> &t,
54  std::string name);
55  template<class vec_t, class vec_size_t>
56  void hdf_output(hdf_file &hf, o2scl::tensor_grid<vec_t,vec_size_t> &t,
57  std::string name);
58 }
59 
60 #ifndef DOXYGEN_NO_O2NS
61 namespace o2scl {
62 #endif
63 
64  /** \brief A <tt>ublas::range</tt> typedef for \ref
65  o2scl::tensor_grid and related classes in src/base/tensor_grid.h
66  */
67  typedef boost::numeric::ublas::range ub_range;
68 
69  /** \brief A <tt>ublas::vector_range</tt> typedef for \ref
70  o2scl::tensor_grid and related classes in src/base/tensor_grid.h
71  */
74 
75  /** \brief A <tt>ublas::vector_range</tt> typedef (size_t version)
76  for \ref o2scl::tensor_grid and related classes
77  in src/base/tensor_grid.h
78  */
81 
82  /** \brief Tensor class with arbitrary dimensions with a grid
83 
84  This tensor class allows one to assign the indexes to numerical
85  scales, effectively defining a data set on an n-dimensional
86  grid. To set the grid, use \ref default_grid(), \ref set_grid()
87  or \ref set_grid_packed().
88 
89  By convention, member functions ending in the <tt>_val</tt>
90  suffix return the closest grid-point to some user-specified
91  values.
92 
93  \comment
94  I like how hist_new only allows you to set the
95  grid if it's the same size as the old grid or the tensor
96  is empty. This same practice could be applied here, to
97  force the user to clear the tensor before resetting the grid.
98  (10/24/11 - Actually, maybe this is a bad idea, because
99  this class is more analogous to ubvector, which doesn't
100  behave this way.)
101  \endcomment
102 
103  \b Slicing
104 
105  New \ref o2scl::tensor_grid objects can be obtained
106  by fixing any set of indices using \ref copy_slice_interp().
107 
108  Fixing all but two indices also results in a \ref o2scl::table3d
109  object, and five functions perform this task in different ways.
110  The function \ref copy_table3d_align() copies a two-dimensional
111  slice to a \ref o2scl::table3d object presuming that the grid in
112  the \ref o2scl::table3d object has already been set and exactly
113  matches the corresponding sizes for the selected tensor indices.
114  This function does not check that the grids between the two
115  objects match, it only ensures that they have the same size. In
116  order to copy to a \ref o2scl::table3d object and set its grid
117  to match that from the unfixed indices in the \ref
118  o2scl::tensor_grid object, the function \ref
119  copy_table3d_align_setxy() can be used. The function \ref
120  copy_table3d_interp() uses interpolation to extract values from
121  the \ref o2scl::tensor_grid object. It allows the user to select
122  indices to be fixed and then uses the values in the grid in the
123  \ref o2scl::table3d object for the indices which vary.
124  Alternatively \ref copy_table3d_interp_values() allows the user
125  to specify values on the grid for the indices to be fixed and
126  uses the grid in the \ref o2scl::table3d object for the indices
127  which vary. Finally, \ref copy_table3d_interp_values_setxy() acts
128  like \ref copy_table3d_interp_values() except that it sets the
129  \ref o2scl::table3d grid to be the same as the grid in the \ref
130  o2scl::tensor_grid object which corresponds to the indices which
131  are being varied.
132 
133  \b Notes and Todos
134 
135  \note Currently, HDF5 I/O is only allowed if the tensor is
136  allocated with std::vector-based types, and the \ref
137  interpolate() function only works with ublas-based vector types.
138 
139  \todo It is possible for the user to create a tensor_grid
140  object, upcast it to a tensor object, and then use
141  tensor::resize() to resize the tensor, failing to resize the
142  grid. Following this, grid access functions will access random
143  parts of memory or segfault. This can be fixed by ensuring that
144  resize functions are virtual and have a version in tensor_grid
145  which ensure that the grid and tensor data are matched. The
146  problem is that the resize functions are templates, so they
147  cannot be virtual.
148 
149  \future A swap function for the grid similar to the data swap
150  function in the parent \ref o2scl::tensor class?
151 
152  \future Only allocate space for grid if it is set.
153 
154  \future As with \ref o2scl::tensor, generalize to other
155  base data types.
156 
157  \future The function \ref interp_linear_partial() appears
158  to be a generalization of \ref copy_table3d_interp_values_setxy(),
159  so there may be some code duplication between the two that
160  can be avoided.
161  */
162  template<class vec_t=std::vector<double>,
163  class vec_size_t=std::vector<size_t> > class tensor_grid :
164  public tensor<double,vec_t,vec_size_t> {
165 
166  public:
167 
168 #ifndef DOXYGEN_INTERNAL
169 
170  protected:
171 
172 #ifdef O2SCL_NEVER_DEFINED
173  }{
174 #endif
175 
176  /// A rank-sized set of arrays for the grid points
177  vec_t grid;
178 
179  /// If true, the grid has been set by the user
180  bool grid_set;
181 
182  /// Interpolation type
183  size_t itype;
184 
185 #endif
186 
187  public:
188 
189  /// \name Constructors and Destructors
190  //@{
191  /// Create an empty tensor with zero rank
192  tensor_grid() : tensor<double,vec_t,vec_size_t>() {
193  grid_set=false;
194  itype=itp_linear;
195  }
196 
197  /** \brief Create a tensor of rank \c rank with sizes given in \c dim
198 
199  The parameter \c dim must be a vector of sizes with length \c
200  rank. If the user requests any of the sizes to be zero, this
201  constructor will call the error handler.
202  */
203  template<class size_vec_t>
204  tensor_grid(size_t rank, const size_vec_t &dim) :
205  tensor<double,vec_t,vec_size_t>(rank,dim) {
206  grid_set=false;
207  itype=itp_linear;
208  // Note that the parent method sets rk to be equal to rank
209  for(size_t i=0;i<this->rk;i++) {
210  if (dim[i]==0) {
211  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
212  "rank for index "+szttos(i)+" in tensor_grid::"+
213  "tensor_grid(size_t,size_vec_t)").c_str(),
214  exc_einval);
215  }
216  }
217  }
218 
219  /** \brief Create a tensor with a grid defined by a set
220  of \ref o2scl::uniform_grid objects
221  */
222  tensor_grid(std::vector<uniform_grid<double> > &ugs) :
223  tensor<double,vec_t,vec_size_t>() {
224  this->rk=ugs.size();
225  itype=itp_linear;
226  size_t tot=1;
227  for(size_t j=0;j<this->rk;j++) {
228  this->size.push_back(ugs[j].get_npoints());
229  tot*=ugs[j].get_npoints();
230  }
231  this->data.resize(tot);
232  set_grid(ugs);
233  }
234 
235  /** \brief Destructor
236  */
237  virtual ~tensor_grid() {
238  }
239  //@}
240 
241  /// \name Method to check for valid object
242  //@{
243  /** \brief Check that the \ref o2scl::tensor_grid object is valid
244  */
245  void is_valid() const {
246 
248 
249  if (this->rk>0 && grid_set) {
250  size_t tot2=0;
251  for(size_t i=0;i<this->rk;i++) {
252  tot2+=this->size[i];
253  }
254 
255  if (tot2!=grid.size()) {
256  O2SCL_ERR2("Value grid_set is true but grid vector size ",
257  "is wrong in tensor_grid::is_valid().",
259  }
260  }
261 
262  if (!grid_set && grid.size()>0) {
263  O2SCL_ERR2("Value grid_set is false but grid vector size ",
264  "is not zero in tensor_grid::is_valid().",
266  }
267 
268  return;
269  }
270  //@}
271 
272  /// \name Copy constructors
273  //@{
274  /** \brief Copy using <tt>operator()</tt>
275  */
278  this->rk=t.rk;
279  this->data=t.data;
280  this->size=t.size;
281  grid=t.grid;
282  grid_set=t.grid_set;
283  itype=t.itype;
284  }
285 
286  /** \brief Copy using <tt>operator=()</tt>
287  */
290  if (this!=&t) {
291  this->rk=t.rk;
292  this->data=t.data;
293  this->size=t.size;
294  grid=t.grid;
295  grid_set=t.grid_set;
296  itype=t.itype;
297  }
298  return *this;
299  }
300  //@}
301 
302  /// \name Set functions
303  //@{
304  /// Set the element closest to grid point \c grdp to value \c val
305  template<class vec2_t>
306  void set_val(const vec2_t &grdp, double val) {
307 
308  // Find indices
309  vec_size_t index(this->rk);
310  for(size_t i=0;i<this->rk;i++) {
311  index[i]=lookup_grid(i,grdp[i]);
312  }
313 
314  // Pack
315  size_t ix=index[0];
316  for(size_t i=1;i<this->rk;i++) {
317  ix*=this->size[i];
318  ix+=index[i];
319  }
320 
321  // Set value
322  this->data[ix]=val;
323 
324  return;
325  }
326 
327  /** \brief Set the element closest to grid point \c grdp to value
328  \c val
329 
330  The parameters \c closest and \c grdp may be identical,
331  allowing one to update a vector \c grdp with the
332  closest grid point.
333  */
334  template<class vec2_t, class vec3_t>
335  void set_val(const vec2_t &grdp, double val, vec3_t &closest) {
336 
337  // Find indices
338  vec_size_t index(this->rk);
339  for(size_t i=0;i<this->rk;i++) {
340  index[i]=lookup_grid_val(i,grdp[i],closest[i]);
341  }
342 
343  // Pack
344  size_t ix=index[0];
345  for(size_t i=1;i<this->rk;i++) {
346  ix*=this->size[i];
347  ix+=index[i];
348  }
349 
350  // Set value
351  this->data[ix]=val;
352 
353  return;
354  }
355  //@}
356 
357  /// \name Get functions
358  //@{
359  /// Get the element closest to grid point \c gridp
360  template<class vec2_t> double get_val(const vec2_t &gridp) {
361 
362  // Find indices
363  vec_size_t index(this->rk);
364  for(size_t i=0;i<this->rk;i++) {
365  index[i]=lookup_grid(i,gridp[i]);
366  }
367 
368  // Pack
369  size_t ix=index[0];
370  for(size_t i=1;i<this->rk;i++) {
371  ix*=this->size[i];
372  ix+=index[i];
373  }
374 
375  // Set value
376  return this->data[ix];
377  }
378 
379  /** \brief Get the element closest to grid point \c gridp,
380  store grid values in \c closest and return value
381 
382  The parameters \c gridp and \c closest may refer to the
383  same object.
384  */
385  template<class vec2_t, class vec3_t>
386  double get_val(const vec2_t &gridp, vec3_t &closest) {
387 
388  // Find indices
389  vec_size_t index(this->rk);
390  for(size_t i=0;i<this->rk;i++) {
391  index[i]=lookup_grid_val(i,gridp[i],closest[i]);
392  }
393 
394  // Pack
395  size_t ix=index[0];
396  for(size_t i=1;i<this->rk;i++) {
397  ix*=this->size[i];
398  ix+=index[i];
399  }
400 
401  // Set value
402  return this->data[ix];
403  }
404 
405  /// Get grid
406  const vec_t &get_grid() const {
407  return grid;
408  }
409  //@}
410 
411  /// \name Resize method
412  //@{
413  /** \brief Resize the tensor to rank \c rank with sizes
414  given in \c dim
415 
416  The parameter \c dim must be a vector of sizes with a length
417  equal to \c rank. This resize method is always destructive,
418  and the grid is always reset.
419 
420  If the user requests any of the sizes to be zero, this
421  function will call the error handler.
422  */
423  template<class size_vec2_t>
424  void resize(size_t rank, const size_vec2_t &dim) {
425  // Double check that none of the sizes that the user
426  // specified are zero
427  for(size_t i=0;i<rank;i++) {
428  if (dim[i]==0) {
429  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
430  "rank for index "+szttos(i)+" in tensor_grid::"+
431  "resize().").c_str(),exc_einval);
432  }
433  }
434  // Set the new rank
435  this->rk=rank;
436  // Resize the size vector
437  this->size.resize(this->rk);
438  // Reset the grid
439  grid_set=false;
440  grid.resize(0);
441  // If the new rank is zero, reset the data, otherwise,
442  // update the size vector and reset the data vector
443  if (rank==0) {
444  this->data.resize(0);
445  return;
446  } else {
447  size_t tot=1;
448  for(size_t i=0;i<this->rk;i++) {
449  this->size[i]=dim[i];
450  tot*=this->size[i];
451  }
452  this->data.resize(tot);
453  }
454  return;
455  }
456 
457  //@}
458 
459  /// \name Grid manipulation
460  //@{
461  /// Return true if the grid has been set
462  bool is_grid_set() const { return grid_set; }
463 
464  /** \brief Set the grid
465 
466  The grid must be specified for all of the dimensions at
467  once. Denote \f$ (\mathrm{size})_0 \f$ as the size
468  of the first dimension, \f$ (\mathrm{size})_1 \f$ as the
469  size of the second dimesion, and so on. Then the
470  first \f$ (\mathrm{size})_0 \f$ entries in \c grid must
471  be the grid for the first dimension, the next
472  \f$ (\mathrm{size})_1 \f$ entries must be the grid
473  for the second dimension, and so on. Thus \c grid
474  must be a vector of size
475  \f[
476  \sum_{i=0}^{\mathrm{rank}} (\mathrm{size})_i
477  \f]
478 
479  Note that the grid is copied so the function argument may
480  be destroyed by the user after calling set_grid_packed() without
481  affecting the tensor grid.
482 
483  \future Define a more generic interface for matrix types
484  */
485  template<class vec2_t>
486  void set_grid_packed(const vec2_t &grid_vec) {
487  if (this->rk==0) {
488  O2SCL_ERR2("Tried to set grid for empty tensor in ",
489  "tensor_grid::set_grid_packed().",exc_einval);
490  }
491  size_t ngrid=0;
492  for(size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
493  grid.resize(ngrid);
494  for(size_t i=0;i<ngrid;i++) {
495  grid[i]=grid_vec[i];
496  }
497  grid_set=true;
498  return;
499  }
500 
501  /** \brief Set grid from a vector of vectors of grid points
502  */
503  template<class vec_vec_t>
504  void set_grid(const vec_vec_t &grid_vecs) {
505  if (this->rk==0) {
506  O2SCL_ERR2("Tried to set grid for empty tensor in ",
507  "tensor_grid::set_grid().",exc_einval);
508  }
509  size_t ngrid=0;
510  for(size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
511  grid.resize(ngrid);
512  size_t k=0;
513  for(size_t i=0;i<this->rk;i++) {
514  for(size_t j=0;j<this->size[i];j++) {
515  grid[k]=grid_vecs[i][j];
516  k++;
517  }
518  }
519  grid_set=true;
520  return;
521  }
522 
523  /** \brief Use a default grid which just uses the index
524  */
525  void default_grid() {
526  size_t ngrid=0;
527  for(size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
528  grid.resize(ngrid);
529  size_t k=0;
530  for(size_t i=0;i<this->rk;i++) {
531  for(size_t j=0;j<this->size[i];j++) {
532  grid[k]=((double)j);
533  k++;
534  }
535  }
536  grid_set=true;
537  return;
538  }
539 
540  /** \brief Set grid for one index from a vector
541  */
542  template<class vec2_t>
543  void set_grid_i_vec(size_t ix, const vec2_t &grid_vec) {
544  if (grid_set==false) {
545  O2SCL_ERR2("Grid not already set in ",
546  "tensor_grid::set_grid_i_vec().",exc_einval);
547  }
548  if (this->rk==0) {
549  O2SCL_ERR2("Tried to set grid for empty tensor in ",
550  "tensor_grid::set_grid_i_vec().",exc_einval);
551  }
552  size_t k=0;
553  for(size_t i=0;i<this->rk;i++) {
554  for(size_t j=0;j<this->size[i];j++) {
555  if (j==ix) {
556  grid[k]=grid_vec[j];
557  }
558  k++;
559  }
560  }
561  return;
562  }
563 
564  /** \brief Set grid for one index from a function
565  */
566  template<class vec2_t>
567  void set_grid_i_func(size_t ix, std::string func) {
568  if (grid_set==false) {
569  O2SCL_ERR2("Grid not already set in ",
570  "tensor_grid::set_grid_i_func().",exc_einval);
571  }
572  if (this->rk==0) {
573  O2SCL_ERR2("Tried to set grid for empty tensor in ",
574  "tensor_grid::set_grid_i_func().",exc_einval);
575  }
576 
577  calculator calc;
578  std::map<std::string,double> vars;
579  calc.compile(func.c_str(),&vars);
580 
581  size_t k=0;
582  for(size_t i=0;i<this->rk;i++) {
583  for(size_t j=0;j<this->size[i];j++) {
584  if (j==ix) {
585  vars["i"]=((double)j);
586  grid[k]=calc.eval(&vars);
587  }
588  k++;
589  }
590  }
591 
592  return;
593  }
594 
595  /** \brief Set grid from a vector of uniform grid objects
596 
597  \note This is called by one of the constructors.
598  */
599  void set_grid(std::vector<uniform_grid<double> > &ugs) {
600  if (this->rk==0) {
601  O2SCL_ERR2("Tried to set grid for empty tensor in ",
602  "tensor_grid::set_grid().",exc_einval);
603  }
604  size_t ngrid=0;
605  for(size_t i=0;i<this->rk;i++) ngrid+=ugs[i].get_npoints();
606  grid.resize(ngrid);
607  size_t k=0;
608  for(size_t i=0;i<this->rk;i++) {
609  for(size_t j=0;j<ugs[i].get_npoints();j++) {
610  grid[k]=ugs[i][j];
611  k++;
612  }
613  }
614  grid_set=true;
615  return;
616  }
617 
618  /** \brief Copy grid for index \c i to vector \c v
619 
620  The type \c rvec_t must be a vector with a resize
621  method.
622  */
623  template<class rvec_t> void copy_grid(size_t i, rvec_t &v) const {
624  v.resize(this->size[i]);
625  size_t istart=0;
626  for(size_t k=0;k<i;k++) istart+=this->size[k];
627  for(size_t k=0;k<this->size[i];k++) {
628  v[k]=grid[istart+k];
629  }
630  return;
631  }
632 
633  /// Lookup jth value on the ith grid
634  double get_grid(size_t i, size_t j) const {
635  if (!grid_set) {
636  O2SCL_ERR("Grid not set in tensor_grid::get_grid().",
637  exc_einval);
638  }
639  if (i>=this->rk) {
640  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
641  " greater than or equal to rank, "+szttos(this->rk)+
642  ", in tensor_grid::get_grid().").c_str(),
643  exc_einval);
644  }
645  size_t istart=0;
646  for(size_t k=0;k<i;k++) istart+=this->size[k];
647  return grid[istart+j];
648  }
649 
650  /// Set the jth value on the ith grid
651  void set_grid(size_t i, size_t j, double val) {
652  if (!grid_set) {
653  O2SCL_ERR("Grid not set in tensor_grid::get_grid().",
654  exc_einval);
655  }
656  if (i>=this->rk) {
657  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
658  " greater than or equal to rank, "+szttos(this->rk)+
659  ", in tensor_grid::get_grid().").c_str(),
660  exc_einval);
661  }
662  size_t istart=0;
663  for(size_t k=0;k<i;k++) istart+=this->size[k];
664  grid[istart+j]=val;
665  return;
666  }
667 
668  /** \brief Lookup index for grid closest to \c val, returning the
669  grid point
670 
671  The parameters \c val and \c val2 may refer to the
672  same object.
673  */
674  size_t lookup_grid_val(size_t i, const double &val, double &val2) const {
675  if (i>=this->rk) {
676  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
677  " greater than or equal to rank, "+szttos(this->rk)+
678  ", in tensor_grid::lookup_grid_val().").c_str(),
679  exc_einval);
680  }
681  if (grid_set==false) {
682  O2SCL_ERR("Grid not set in tensor_grid::lookup_grid_val().",
683  exc_einval);
684  }
685 
686  // We have to store the grid point in a temporary in case
687  // the parameters gridp and closest refer to the same object
688  double temp=val;
689 
690  size_t istart=0;
691  for(size_t j=0;j<i;j++) istart+=this->size[j];
692  size_t best=istart;
693  double min=fabs(grid[istart]-temp);
694  val2=grid[istart];
695  for(size_t j=istart;j<istart+this->size[i];j++) {
696  if (fabs(grid[j]-temp)<min) {
697  best=j;
698  min=fabs(grid[j]-temp);
699  val2=grid[j];
700  }
701  }
702  return best-istart;
703  }
704 
705  /// Lookup index for grid closest to \c val
706  size_t lookup_grid(size_t i, double val) const {
707  double val2;
708  return lookup_grid_val(i,val,val2);
709  }
710 
711  /** \brief Lookup indices for grid closest point to \c vals
712 
713  The values in \c vals are not modified by this function.
714 
715  \comment
716  This function must have a different name than
717  lookup_grid() because the template types cause
718  confusion between the two functions.
719  \endcomment
720  */
721  template<class vec2_t, class size_vec2_t>
722  void lookup_grid_vec(const vec2_t &vals, size_vec2_t &indices) const {
723  for(size_t k=0;k<this->rk;k++) {
724  indices[k]=lookup_grid(k,vals[k]);
725  }
726  return;
727  }
728 
729  /** \brief Lookup internal packed grid index for point closest
730  to \c val and store closest value in \c val2
731 
732  This version, rather than \ref
733  o2scl::tensor_grid::lookup_grid_val() can be useful because it
734  gives the index of the grid point in the internal grid vector
735  object.
736  */
737  size_t lookup_grid_packed_val(size_t i, double val, double &val2) const {
738  if (!grid_set) {
739  O2SCL_ERR("Grid not set in tensor_grid::lookup_grid_packed_val().",
740  exc_einval);
741  }
742  if (i>=this->rk) {
743  O2SCL_ERR((((std::string)"Index ")+szttos(i)+" greater than rank, "+
744  szttos(this->rk)+
745  ", in tensor_grid::lookup_grid_packed_val().").c_str(),
746  exc_einval);
747  }
748  size_t istart=0;
749  for(size_t j=0;j<i;j++) istart+=this->size[j];
750  size_t best=istart;
751  double min=fabs(grid[istart]-val);
752  val2=grid[istart];
753  for(size_t j=istart;j<istart+this->size[i];j++) {
754  if (fabs(grid[j]-val)<min) {
755  best=j;
756  min=fabs(grid[j]-val);
757  val2=grid[j];
758  }
759  }
760  return best;
761  }
762 
763  /** \brief Lookup internal packed grid index for point closest
764  to \c val
765  */
766  size_t lookup_grid_packed(size_t i, double val) const {
767  double val2;
768  return lookup_grid_packed_val(i,val,val2);
769  }
770  //@}
771 
772  /// \name Slicing to tensor_grid objects
773  //@{
774  /** \brief Copy an abitrary slice by fixing 1 or more indices and
775  use interpolation to return a new \ref tensor_grid object
776  */
777  template<class size_vec2_t, class vec2_t>
778  tensor_grid<> copy_slice_interp(size_vec2_t &ifix, vec2_t &vals) const {
779 
780  if (this->rk<1+ifix.size()) {
781  O2SCL_ERR2("Fixed too many indices in ",
782  "tensor_grid::copy_slice_interp().",
784  }
785  if (ifix.size()!=vals.size()) {
786  O2SCL_ERR2("Mismatch between indices and values in ",
787  "tensor_grid::copy_slice_interp().",
789  }
790 
791  // Determine new rank
792  size_t rank_new=this->rk-ifix.size();
793 
794  // Determine the new sizes and new grid
795  std::vector<size_t> sz_new;
796  std::vector<std::vector<double> > grid_new;
797  for(size_t i=0;i<this->rk;i++) {
798  bool found=false;
799  for(size_t j=0;j<ifix.size();j++) {
800  if (ifix[j]==i) found=true;
801  }
802  if (found==false) {
803  sz_new.push_back(this->get_size(i));
804  std::vector<double> grid_temp;
805  for(size_t j=0;j<this->get_size(i);j++) {
806  grid_temp.push_back(this->get_grid(i,j));
807  }
808  grid_new.push_back(grid_temp);
809  }
810  }
811 
812  // Create the new tensor_grid object and set the new grid
813  tensor_grid<> tg_new(rank_new,sz_new);
814  tg_new.set_grid(grid_new);
815 
816  // Interpolate the data into the new tensor_grid object
817  std::vector<size_t> ix_new(rank_new);
818  std::vector<double> point_old(this->rk);
819 
820  // Loop over the new tensor_grid object
821  for(size_t i=0;i<tg_new.total_size();i++) {
822 
823  // Find the location in the new tensor_grid object
824  tg_new.unpack_index(i,ix_new);
825 
826  // Find the point in the old tensor object to interpolate
827  for(size_t j=0;j<this->rk;j++) {
828  int ix_found=-1;
829  for(size_t k=0;k<ifix.size();k++) {
830  if (ifix[k]==j) ix_found=k;
831  }
832  if (ix_found==-1) {
833  point_old[j]=this->get_grid(j,ix_new[j]);
834  } else {
835  point_old[j]=vals[ix_found];
836  }
837  }
838 
839  // Set the new point by performing the linear interpolation
840  tg_new.set(ix_new,this->interp_linear(point_old));
841  }
842 
843  return tg_new;
844  }
845  //@}
846 
847  /// \name Slicing and converting to table3d objects
848  //@{
849 
850  /** \brief Convert to a \ref o2scl::table3d object by
851  summing over all but two indices
852  */
853  void convert_table3d_sum
854  (size_t ix_x, size_t ix_y, table3d &tab, std::string x_name="x",
855  std::string y_name="y", std::string slice_name="z") const {
856 
857  // Get current table3d grid
858  size_t nx, ny;
859  tab.get_size(nx,ny);
860 
861  if (nx==0 && ny==0) {
862 
863  if (x_name.length()==0) x_name="x";
864  if (y_name.length()==0) y_name="y";
865 
866  // If there's no grid, then create a grid in the table3d
867  // object that is the same as that in the tensor_grid object
868  std::vector<double> grid_x, grid_y;
869  copy_grid(ix_x,grid_x);
870  copy_grid(ix_y,grid_y);
871  tab.set_xy("x",grid_x.size(),grid_x,
872  "y",grid_y.size(),grid_y);
873  // Now that the grid is set, get nx and ny
874  tab.get_size(nx,ny);
875  }
876 
877  // Check that the grids are commensurate
878  if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
879  O2SCL_ERR2("Grids not commensurate in ",
880  "tensor_grid::convert_table3d_sum().",exc_einval);
881  }
882 
883  tab.set_slice_all(slice_name,0.0);
884 
885  std::vector<size_t> ix;
886  for(size_t i=0;i<this->total_size();i++) {
887  this->unpack_index(i,ix);
888  tab.set(ix[ix_x],ix[ix_y],slice_name,
889  tab.get(ix[ix_x],ix[ix_y],slice_name)+
890  this->data[i]);
891  }
892 
893  return;
894  }
895 
896  /** \brief Create a slice in a \ref o2scl::table3d object with an
897  aligned grid
898 
899  This function uses the grid associated with indices \c ix_x
900  and \c ix_y, to copy data to a slice named \c slice_name in
901  the table3d object \c tab . All other indices are fixed
902  to values specified by the user in \c index and the values
903  of <tt>index[ix_x]</tt> and <tt>index[ix_y]</tt> are
904  used for temporary storage.
905 
906  If the table3d object does not currently have a grid set, then
907  the grid is automatically set to be the same as that stored in
908  the tensor_grid object associated with ranks \c ix_x and \c
909  iy_y. If the \ref o2scl::table3d object does have a grid set,
910  then the values returned by \ref o2scl::table3d::get_nx() and
911  \ref o2scl::table3d::get_ny() must be equal to the size of the
912  tensor in indices \c ix_x and ix_y, respectively. If a
913  slice named \c slice_name is not already present in
914  \c tab, then a new slice with that name is created.
915 
916  The error handler is called if \c ix_x is the same as
917  \c ix_y, or if either of these two values is greater
918  than or equal to the tensor rank.
919  */
920  template<class size_vec2_t>
921  void copy_table3d_align(size_t ix_x, size_t ix_y, size_vec2_t &index,
922  table3d &tab, std::string slice_name="z") const {
923 
924  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
925  O2SCL_ERR2("Either indices greater than rank or x and y ind",
926  "ices equal in tensor_grid::copy_table3d_align().",
927  exc_efailed);
928  }
929 
930  // Get current table3d grid
931  size_t nx, ny;
932  tab.get_size(nx,ny);
933 
934  // Check that the grids are commensurate
935  if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
936  O2SCL_ERR2("Grids not commensurate in ",
937  "tensor_grid::copy_table3d_align().",exc_einval);
938  }
939 
940  // Create slice if not already present
941  size_t is;
942  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
943 
944  // Copy over data
945  for(size_t i=0;i<nx;i++) {
946  for(size_t j=0;j<ny;j++) {
947  index[ix_x]=i;
948  index[ix_y]=j;
949  double val=this->get(index);
950  tab.set(i,j,slice_name,val);
951  }
952  }
953 
954  return;
955  }
956 
957  /** \brief Create a slice in a table3d object with a new aligned
958  grid
959  */
960  template<class size_vec2_t>
961  void copy_table3d_align_setxy
962  (size_t ix_x, size_t ix_y, size_vec2_t &index,
963  table3d &tab, std::string x_name="x", std::string y_name="y",
964  std::string slice_name="z") const {
965 
966  // Get current table3d grid
967  size_t nx, ny;
968  tab.get_size(nx,ny);
969 
970  if (nx==0 && ny==0) {
971 
972  if (x_name.length()==0) x_name="x";
973  if (y_name.length()==0) y_name="y";
974 
975  // If there's no grid, then create a grid in the table3d
976  // object that is the same as that in the tensor_grid object
977  std::vector<double> grid_x, grid_y;
978  copy_grid(ix_x,grid_x);
979  copy_grid(ix_y,grid_y);
980  tab.set_xy("x",grid_x.size(),grid_x,
981  "y",grid_y.size(),grid_y);
982  // Now that the grid is set, get nx and ny
983  tab.get_size(nx,ny);
984  }
985 
986  copy_table3d_align(ix_x,ix_y,index,tab,slice_name);
987 
988  return;
989  }
990 
991  /** \brief Copy to a slice in a table3d object using interpolation
992 
993  This function uses the grid associated with indices \c ix_x
994  and \c ix_y, and the tensor interpolation function to copy the
995  tensor information to the slice named \c slice_name in the
996  table3d object \c tab . All other indices are fixed
997  to values specified by the user in \c index and the values
998  of <tt>index[ix_x]</tt> and <tt>index[ix_y]</tt> are
999  used for temporary storage.
1000 
1001  If a slice named \c slice_name is not already present in \c
1002  tab, then a new slice with that name is created.
1003 
1004  The error handler is called if \c ix_x is the same as
1005  \c ix_y, or if either of these two values is greater
1006  than or equal to the tensor rank.
1007 
1008  \note This function uses the \ref tensor_grid::interp_linear()
1009  for the interpolation.
1010  */
1011  template<class size_vec2_t>
1012  void copy_table3d_interp(size_t ix_x, size_t ix_y, size_vec2_t &index,
1013  table3d &tab, std::string slice_name="z") const {
1014 
1015  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
1016  O2SCL_ERR2("Either indices greater than rank or x and y ",
1017  "indices equal in tensor_grid::copy_table3d_interp().",
1018  exc_efailed);
1019  }
1020 
1021  // Get current table3d grid
1022  size_t nx, ny;
1023  tab.get_size(nx,ny);
1024 
1025  if (nx==0 && ny==0) {
1026  // If there's no grid, then just use the aligned version
1027  return copy_table3d_align(ix_x,ix_y,index,tab,slice_name);
1028  }
1029 
1030  // Create vector of values to interpolate with
1031  std::vector<double> vals(this->rk);
1032  for(size_t i=0;i<this->rk;i++) {
1033  if (i!=ix_x && i!=ix_y) vals[i]=this->get_grid(i,index[i]);
1034  }
1035 
1036  // Create slice if not already present
1037  size_t is;
1038  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
1039 
1040  // Loop through the table grid to perform the interpolation
1041  for(size_t i=0;i<nx;i++) {
1042  for(size_t j=0;j<ny;j++) {
1043  vals[ix_x]=tab.get_grid_x(i);
1044  vals[ix_y]=tab.get_grid_y(j);
1045  tab.set(i,j,slice_name,this->interp_linear(vals));
1046  }
1047  }
1048 
1049  return;
1050  }
1051 
1052  /** \brief Copy to a slice in a table3d object using interpolation
1053  */
1054  template<class vec2_t>
1055  void copy_table3d_interp_values(size_t ix_x, size_t ix_y,
1056  vec2_t &values, table3d &tab,
1057  std::string slice_name="z",
1058  int verbose=0) const {
1059 
1060  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
1061  O2SCL_ERR2("Either indices greater than rank or x and y ",
1062  "indices equal in tensor_grid::copy_table3d_interp().",
1063  exc_efailed);
1064  }
1065  if (values.size()!=this->rk) {
1066  O2SCL_ERR2("Values array not equal to rank ",
1067  "in tensor_grid::copy_table3d_interp_values().",
1068  exc_efailed);
1069  }
1070 
1071  if (tab.is_size_set()==false || tab.is_xy_set()==false) {
1072  O2SCL_ERR2("Grid not set in tensor_grid::",
1073  "copy_table3d_interp_value().",o2scl::exc_einval);
1074  }
1075 
1076  // Get current table3d grid
1077  size_t nx, ny;
1078  tab.get_size(nx,ny);
1079 
1080  // Create slice if not already present
1081  size_t is;
1082  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
1083 
1084  // Loop through the table grid to perform the interpolation
1085  for(size_t i=0;i<nx;i++) {
1086  for(size_t j=0;j<ny;j++) {
1087  values[ix_x]=tab.get_grid_x(i);
1088  values[ix_y]=tab.get_grid_y(j);
1089  tab.set(i,j,slice_name,this->interp_linear(values));
1090  if (verbose>0) {
1091  std::cout << "At location values: ";
1092  for(size_t k=0;k<values.size();k++) {
1093  std::cout << values[k] << " ";
1094  }
1095  std::cout << "Interpolated to get: "
1096  << i << " " << j << " " << slice_name << " "
1097  << this->interp_linear(values) << std::endl;
1098  if (verbose>1) {
1099  char ch;
1100  std::cin >> ch;
1101  }
1102  }
1103  }
1104  }
1105 
1106  return;
1107  }
1108 
1109  /** \brief Copy to a slice in a table3d object using interpolation
1110  creating a new table3d grid
1111  */
1112  template<class vec2_t>
1113  void copy_table3d_interp_values_setxy
1114  (size_t ix_x, size_t ix_y, vec2_t &values, table3d &tab,
1115  std::string x_name="x", std::string y_name="y",
1116  std::string slice_name="z") const {
1117 
1118  // Get current table3d grid
1119  size_t nx, ny;
1120  tab.get_size(nx,ny);
1121 
1122  if (nx==0 && ny==0) {
1123 
1124  if (x_name.length()==0) x_name="x";
1125  if (y_name.length()==0) y_name="y";
1126 
1127  // If there's no grid, then create a grid in the table3d
1128  // object that is the same as that in the tensor_grid object
1129  std::vector<double> grid_x, grid_y;
1130  copy_grid(ix_x,grid_x);
1131  copy_grid(ix_y,grid_y);
1132  tab.set_xy("x",grid_x.size(),grid_x,
1133  "y",grid_y.size(),grid_y);
1134  // Now that the grid is set, get nx and ny
1135  tab.get_size(nx,ny);
1136  }
1137 
1138  copy_table3d_interp_values(ix_x,ix_y,values,tab,slice_name);
1139  return;
1140  }
1141  //@}
1142 
1143  /// \name Clear method
1144  //@{
1145  /// Clear the tensor of all data and free allocated memory
1146  void clear() {
1147  grid.resize(0);
1148  grid_set=false;
1150  return;
1151  }
1152  //@}
1153 
1154  /// \name Interpolation
1155  //@{
1156  /// Set interpolation type for \ref interpolate()
1157  void set_interp_type(size_t interp_type) {
1158  itype=interp_type;
1159  return;
1160  }
1161 
1162  /** \brief Interpolate values \c vals into the tensor,
1163  returning the result
1164 
1165  This is a quick and dirty implementation of n-dimensional
1166  interpolation by recursive application of the 1-dimensional
1167  routine from \ref interp_vec, using the base interpolation
1168  object specified in the template parameter \c base_interp_t.
1169  This will be very slow for sufficiently large data sets.
1170 
1171  \note This function requires a range objects to obtain ranges
1172  of vector objects. In ublas, this is done with
1173  <tt>ublas::vector_range</tt> objects, so this function will
1174  certainly work for \ref tensor_grid objects built on ublas
1175  vector types. There is no corresponding <tt>std::range</tt>,
1176  but you may be able to use either <tt>ublas::vector_range</tt>
1177  or <tt>Boost.Range</tt> in order to call this function with
1178  \ref tensor_grid objects built on <tt>std::vector</tt>.
1179  However, this is not fully tested at the moment.
1180 
1181  \future It should be straightforward to improve the scaling of
1182  this algorithm significantly by creating a "window" of local
1183  points around the point of interest. This could be done easily
1184  by constructing an initial subtensor. However, this should
1185  probably be superceded by a more generic alternative which
1186  avoids explicit use of the 1-d interpolation types.
1187  */
1188  template<class range_t=ub_range,
1189  class data_range_t=ubvector_range,
1190  class index_range_t=ubvector_size_t_range>
1191  double interpolate(double *vals) {
1192 
1193  typedef interp_vec<vec_t> interp_t;
1194 
1195  if (this->rk==1) {
1196 
1197  interp_t si(this->size[0],grid,this->data,itype);
1198  return si.eval(vals[0]);
1199 
1200  } else {
1201 
1202  // Get total number of interpolations at this level
1203  size_t ss=1;
1204  for(size_t i=1;i<this->rk;i++) ss*=this->size[i];
1205 
1206  // Create space for y vectors and interpolators
1207  std::vector<vec_t> yvec(ss);
1208  std::vector<interp_t *> si(ss);
1209  for(size_t i=0;i<ss;i++) {
1210  yvec[i].resize(this->size[0]);
1211  }
1212 
1213  // Create space for interpolation results
1214  tensor_grid tdat;
1215  index_range_t size_new(this->size,ub_range(1,this->rk));
1216  tdat.resize(this->rk-1,size_new);
1217 
1218  // Set grid for temporary tensor
1219  data_range_t grid_new(grid,ub_range(this->size[0],grid.size()));
1220  tdat.set_grid_packed(grid_new);
1221 
1222  // Create starting coordinate and counter
1223  vec_size_t co(this->rk);
1224  for(size_t i=0;i<this->rk;i++) co[i]=0;
1225  size_t cnt=0;
1226 
1227  // Loop over every interpolation
1228  bool done=false;
1229  while(done==false) {
1230 
1231  // Fill yvector with the appropriate data
1232  for(size_t i=0;i<this->size[0];i++) {
1233  co[0]=i;
1234  yvec[cnt][i]=this->get(co);
1235  }
1236 
1237  si[cnt]=new interp_t(this->size[0],grid,yvec[cnt],itype);
1238 
1239  index_range_t co2(co,ub_range(1,this->rk));
1240  tdat.set(co2,si[cnt]->eval(vals[0]));
1241 
1242  // Go to next interpolation
1243  cnt++;
1244  co[this->rk-1]++;
1245  // carry if necessary
1246  for(int j=((int)this->rk)-1;j>0;j--) {
1247  if (co[j]>=this->size[j]) {
1248  co[j]=0;
1249  co[j-1]++;
1250  }
1251  }
1252 
1253  // Test if done
1254  if (cnt==ss) done=true;
1255 
1256  // End of while loop
1257  }
1258 
1259  // Now call the next level of interpolation
1260  double res=tdat.interpolate(vals+1);
1261 
1262  for(size_t i=0;i<ss;i++) {
1263  delete si[i];
1264  }
1265 
1266  return res;
1267  }
1268  }
1269 
1270  /** \brief Obtain a value by looking up some indices and
1271  interpolating the others
1272 
1273  To call this function, the arguments should be
1274  of the following form
1275  - The vector \c ix_to_interp should be a list of indices to
1276  interpolate. The size of \c ix_to_interp must be at least 1
1277  or larger but smaller than or equal to the full tensor rank.
1278  All entries in \c ix_to_interp should be smaller than the
1279  full tensor rank.
1280  - The vector \c ix should have a size equal to the tensor
1281  rank, but values stored in entries corresponding to the
1282  indices in \c ix_to_interp will be ignored
1283  - The vector \c val should be a list of values to be
1284  interpolated and should have a size equal to that of
1285  \c ix_to_interp .
1286  */
1287  template<class vec2_size_t, class vec3_size_t, class vec2_t>
1288  double interp_linear_partial
1289  (const vec2_size_t &ix_to_interp,
1290  vec3_size_t &ix, const vec2_t &val) const {
1291 
1292  if (val.size()!=ix_to_interp.size()) {
1293  O2SCL_ERR2("Index and value list don't match in ",
1294  "tensor_grid::interp_linear_partial().",
1296  }
1297  if (ix_to_interp.size()>this->get_rank() ||
1298  ix_to_interp.size()==0) {
1299  O2SCL_ERR2("Index list too large or too small in ",
1300  "tensor_grid::interp_linear_partial().",
1302  }
1303 
1304  // Find the the corner of the hypercube containing val for all
1305  // the indices to be interpolated
1306  std::vector<size_t> loc(ix_to_interp.size());
1307  std::vector<double> gnew;
1308  for(size_t i=0;i<ix_to_interp.size();i++) {
1309  size_t ixi=ix_to_interp[i];
1310  if (ixi>=this->get_rank()) {
1311  O2SCL_ERR2("Index to interpolate larger than tensor rank in ",
1312  "tensor_grid::interp_linear_partial().",
1314  }
1315  std::vector<double> grid_one(this->size[ixi]);
1316  for(size_t j=0;j<this->size[ixi];j++) {
1317  grid_one[j]=this->get_grid(ixi,j);
1318  }
1319  search_vec<std::vector<double> > sv(this->size[ixi],grid_one);
1320  loc[i]=sv.find(val[i]);
1321  gnew.push_back(grid_one[loc[i]]);
1322  gnew.push_back(grid_one[loc[i]+1]);
1323  }
1324 
1325  // Now construct a 2^{rk}-sized tensor containing only that
1326  // hypercube
1327  std::vector<size_t> snew(ix_to_interp.size());
1328  for(size_t i=0;i<ix_to_interp.size();i++) {
1329  snew[i]=2;
1330  }
1331  tensor_grid tnew(ix_to_interp.size(),snew);
1332  tnew.set_grid_packed(gnew);
1333 
1334  // Copy over the relevant data
1335  for(size_t i=0;i<tnew.total_size();i++) {
1336  std::vector<size_t> index_new(ix_to_interp.size());
1337  tnew.unpack_index(i,index_new);
1338  for(size_t j=0;j<ix_to_interp.size();j++) {
1339  ix[ix_to_interp[j]]=index_new[j]+loc[j];
1340  }
1341  tnew.set(index_new,this->get(ix));
1342  }
1343 
1344  // Now use interp_power_two()
1345  return tnew.interp_linear_power_two(val);
1346  }
1347 
1348  /** \brief Perform a linear interpolation of \c v into the
1349  function implied by the tensor and grid
1350 
1351  This performs multi-dimensional linear interpolation (or
1352  extrapolation) It works by first using \ref o2scl::search_vec
1353  to find the interval containing (or closest to) the specified
1354  point in each direction and constructing the corresponding
1355  hypercube of size \f$ 2^{\mathrm{rank}} \f$ containing \c v.
1356  It then calls \ref interp_linear_power_two() to perform the
1357  interpolation in that hypercube.
1358 
1359  \future This starts with a small copy, which can be eliminated
1360  by creating a new version of interp_linear_power_two
1361  which accepts an offset vector parameter so that the
1362  first interpolation is offset. Remaining interpolations
1363  don't need to be offset because the tensor has to be
1364  created from the previous interpolation round.
1365  */
1366  template<class vec2_t> double interp_linear(vec2_t &v) const {
1367 
1368  // Find the the corner of the hypercube containing v
1369  size_t rgs=0;
1370  std::vector<size_t> loc(this->rk);
1371  std::vector<double> gnew;
1372  for(size_t i=0;i<this->rk;i++) {
1373  std::vector<double> grid_unpacked(this->size[i]);
1374  for(size_t j=0;j<this->size[i];j++) {
1375  grid_unpacked[j]=grid[j+rgs];
1376  }
1377  search_vec<std::vector<double> > sv(this->size[i],grid_unpacked);
1378  loc[i]=sv.find(v[i]);
1379  gnew.push_back(grid_unpacked[loc[i]]);
1380  gnew.push_back(grid_unpacked[loc[i]+1]);
1381  rgs+=this->size[i];
1382  }
1383 
1384  // Now construct a 2^{rk}-sized tensor containing only that
1385  // hypercube
1386  std::vector<size_t> snew(this->rk);
1387  for(size_t i=0;i<this->rk;i++) {
1388  snew[i]=2;
1389  }
1390  tensor_grid tnew(this->rk,snew);
1391  tnew.set_grid_packed(gnew);
1392 
1393  // Copy over the relevant data
1394  for(size_t i=0;i<tnew.total_size();i++) {
1395  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1396  tnew.unpack_index(i,index_new);
1397  for(size_t j=0;j<this->rk;j++) index_old[j]=index_new[j]+loc[j];
1398  tnew.set(index_new,this->get(index_old));
1399  }
1400 
1401  // Now use interp_power_two()
1402  return tnew.interp_linear_power_two(v);
1403  }
1404 
1405  /** \brief Perform linear interpolation assuming that all
1406  indices can take only two values
1407 
1408  This function works by recursively slicing the hypercube of
1409  size \f$ 2^{\mathrm{rank}} \f$ into a hypercube of size \f$
1410  2^{\mathrm{rank-1}} \f$ performing linear interpolation for
1411  each pair of points.
1412 
1413  \note This is principally a function for internal use
1414  by \ref interp_linear().
1415  */
1416  template<class vec2_t>
1417  double interp_linear_power_two(vec2_t &v) const {
1418 
1419  if (this->rk==1) {
1420  return this->data[0]+(this->data[1]-this->data[0])/
1421  (grid[1]-grid[0])*(v[0]-grid[0]);
1422  }
1423 
1424  size_t last=this->rk-1;
1425  double frac=(v[last]-get_grid(last,0))/
1426  (get_grid(last,1)-get_grid(last,0));
1427 
1428  // Create new size vector and grid
1429  tensor_grid tnew(this->rk-1,this->size);
1430  tnew.set_grid_packed(grid);
1431 
1432  // Create data in new tensor, removing the last index through
1433  // linear interpolation
1434  for(size_t i=0;i<tnew.total_size();i++) {
1435  std::vector<size_t> index(this->rk);
1436  tnew.unpack_index(i,index);
1437  index[this->rk-1]=0;
1438  double val_lo=this->get(index);
1439  index[this->rk-1]=1;
1440  double val_hi=this->get(index);
1441  tnew.set(index,val_lo+frac*(val_hi-val_lo));
1442  }
1443 
1444  // Recursively interpolate the smaller tensor
1445  return tnew.interp_linear_power_two(v);
1446  }
1447 
1448  /** \brief Perform a linear interpolation of <tt>v[1]</tt>
1449  to <tt>v[n-1]</tt> resulting in a vector
1450 
1451  \note The type <tt>vec2_t</tt> for the vector <tt>res</tt>
1452  must have a <tt>resize()</tt> method.
1453 
1454  This performs multi-dimensional linear interpolation (or
1455  extrapolation) in the last <tt>n-1</tt> indices of the
1456  rank-<tt>n</tt> tensor leaving the first index free and places
1457  the results in the vector \c res.
1458  */
1459  template<class vec2_t, class vec3_t>
1460  void interp_linear_vec0(vec2_t &v, vec3_t &res) const {
1461 
1462  // Find the the corner of the hypercube containing v
1463  size_t rgs=0;
1464  std::vector<size_t> loc(this->rk);
1465  std::vector<double> gnew;
1466  for(size_t i=0;i<this->size[0];i++) {
1467  gnew.push_back(grid[i]);
1468  }
1469  rgs=this->size[0];
1470  loc[0]=0;
1471  for(size_t i=1;i<this->rk;i++) {
1472  std::vector<double> grid_unpacked(this->size[i]);
1473  for(size_t j=0;j<this->size[i];j++) {
1474  grid_unpacked[j]=grid[j+rgs];
1475  }
1476  search_vec<std::vector<double> > sv(this->size[i],grid_unpacked);
1477  loc[i]=sv.find(v[i]);
1478  gnew.push_back(grid_unpacked[loc[i]]);
1479  gnew.push_back(grid_unpacked[loc[i]+1]);
1480  rgs+=this->size[i];
1481  }
1482 
1483  // Now construct a n*2^{rk-1}-sized tensor containing only that
1484  // hypercube
1485  std::vector<size_t> snew(this->rk);
1486  snew[0]=this->size[0];
1487  for(size_t i=1;i<this->rk;i++) {
1488  snew[i]=2;
1489  }
1490  tensor_grid tnew(this->rk,snew);
1491  tnew.set_grid_packed(gnew);
1492 
1493  // Copy over the relevant data
1494  for(size_t i=0;i<tnew.total_size();i++) {
1495  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1496  tnew.unpack_index(i,index_new);
1497  for(size_t j=0;j<this->rk;j++) {
1498  index_old[j]=index_new[j]+loc[j];
1499  }
1500  tnew.set(index_new,this->get(index_old));
1501  }
1502 
1503  // Now use interp_power_two_vec0()
1504  tnew.interp_linear_power_two_vec0(v,res);
1505 
1506  return;
1507  }
1508 
1509  /** \brief Perform linear interpolation assuming that the last
1510  <tt>n-1</tt> indices can take only two values
1511 
1512  This function performs linear interpolation assuming that the
1513  last <tt>n-1</tt> indices can take only two values and placing
1514  the result into <tt>res</tt>.
1515 
1516  \note The type <tt>vec2_t</tt> for the vector <tt>res</tt>
1517  must have a <tt>resize()</tt> method. This is principally a
1518  function for internal use by \ref interp_linear_vec0().
1519  */
1520  template<class vec2_t, class vec3_t>
1521  void interp_linear_power_two_vec0(vec2_t &v, vec3_t &res) const {
1522 
1523  if (this->rk==2) {
1524  size_t n=this->size[0];
1525  res.resize(n);
1526  vec_size_t ix0(2), ix1(2);
1527  ix0[1]=0;
1528  ix1[1]=1;
1529  for(size_t i=0;i<n;i++) {
1530  ix0[0]=i;
1531  ix1[0]=i;
1532  res[i]=this->get(ix0)+(this->get(ix1)-this->get(ix0))/
1533  (grid[n+1]-grid[n])*(v[1]-grid[n]);
1534  }
1535  return;
1536  }
1537 
1538  size_t last=this->rk-1;
1539  double frac=(v[last]-get_grid(last,0))/
1540  (get_grid(last,1)-get_grid(last,0));
1541 
1542  // Create new size vector and grid
1543  tensor_grid tnew(this->rk-1,this->size);
1544  tnew.set_grid_packed(grid);
1545 
1546  // Create data in new tensor, removing the last index through
1547  // linear interpolation
1548  for(size_t i=0;i<tnew.total_size();i++) {
1549  std::vector<size_t> index(this->rk);
1550  tnew.unpack_index(i,index);
1551  index[this->rk-1]=0;
1552  double val_lo=this->get(index);
1553  index[this->rk-1]=1;
1554  double val_hi=this->get(index);
1555  tnew.set(index,val_lo+frac*(val_hi-val_lo));
1556  }
1557 
1558  // Recursively interpolate the smaller tensor
1559  tnew.interp_linear_power_two_vec0(v,res);
1560 
1561  return;
1562  }
1563 
1564  /** \brief Perform a linear interpolation of <tt>v</tt>
1565  into the tensor leaving one index free resulting in a vector
1566 
1567  This performs multi-dimensional linear interpolation (or
1568  extrapolation) in the last <tt>n-1</tt> indices of the
1569  rank-<tt>n</tt> tensor leaving the first index free and places
1570  the results in the vector \c res.
1571 
1572  \future This function could be more efficient.
1573  */
1574  template<class vec2_t, class vec3_t>
1575  void interp_linear_vec(vec2_t &v, size_t ifree, vec3_t &res) const {
1576 
1577  size_t n=this->size[ifree];
1578 
1579  // This function uses interp_linear_power_two_vec0(), so it
1580  // works by remapping the indices. This defines the remapping.
1581  std::vector<size_t> map;
1582  map.push_back(ifree);
1583  for(size_t i=0;i<this->rk;i++) {
1584  if (i!=ifree) {
1585  map.push_back(i);
1586  }
1587  }
1588 
1589  // Find the the corner of the hypercube containing v
1590  size_t rgs=0;
1591  vec_size_t loc(this->rk);
1592  loc[ifree]=0;
1593  for(size_t i=0;i<this->rk;i++) {
1594  vec_t grid_unpacked(this->size[i]);
1595  for(size_t j=0;j<this->size[i];j++) {
1596  grid_unpacked[j]=grid[j+rgs];
1597  }
1598  search_vec<vec_t> sv(this->size[i],grid_unpacked);
1599  if (i!=ifree) {
1600  loc[i]=sv.find(v[i]);
1601  }
1602  rgs+=this->size[i];
1603  }
1604 
1605  // Compute the remapped grid and interpolating vector
1606  std::vector<double> gnew, vnew;
1607  for(size_t new_ix=0;new_ix<this->rk;new_ix++) {
1608  for(size_t old_ix=0;old_ix<this->rk;old_ix++) {
1609  if (map[new_ix]==old_ix) {
1610  vnew.push_back(v[old_ix]);
1611  if (old_ix==ifree) {
1612  for(size_t j=0;j<this->size[old_ix];j++) {
1613  gnew.push_back(this->get_grid(old_ix,j));
1614  }
1615  } else {
1616  gnew.push_back(this->get_grid(old_ix,loc[old_ix]));
1617  gnew.push_back(this->get_grid(old_ix,loc[old_ix]+1));
1618  }
1619  }
1620  }
1621  }
1622 
1623  // Now construct a n*2^{rk-1}-sized tensor containing only the
1624  // hypercube needed to do the interpolation
1625 
1626  // Specify the size of each rank
1627  std::vector<size_t> snew;
1628  snew.push_back(n);
1629  for(size_t i=0;i<this->rk;i++) {
1630  if (i!=ifree) {
1631  snew.push_back(2);
1632  }
1633  }
1634 
1635  // Create the tensor and set the grid
1636  tensor_grid tnew(this->rk,snew);
1637  tnew.set_grid_packed(gnew);
1638 
1639  // Copy over the relevant data
1640  for(size_t i=0;i<tnew.total_size();i++) {
1641  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1642  tnew.unpack_index(i,index_new);
1643  for(size_t j=0;j<this->rk;j++) {
1644  index_old[map[j]]=index_new[j]+loc[map[j]];
1645  }
1646  tnew.set(index_new,this->get(index_old));
1647  }
1648 
1649  // Now use interp_power_two_vec()
1650  tnew.interp_linear_power_two_vec0(vnew,res);
1651 
1652  return;
1653  }
1654  //@}
1655 
1656  template<class vecf_t, class vecf_size_t> friend void o2scl_hdf::hdf_output
1658  std::string name);
1659 
1660  template<class vecf_t, class vecf_size_t> friend void o2scl_hdf::hdf_input
1662  std::string name);
1663 
1664  /** \brief Rearrange, sum and copy current tensor to a new tensor
1665 
1666  \note This function doesn't work yet. Most of all, it does
1667  not yet properly construct the new grid.
1668 
1669  \future Some code duplication between this function
1670  and the one in the tensor class.
1671  */
1672  tensor_grid<> rearrange_and_copy(std::vector<index_spec> spec,
1673  int verbose=0,
1674  bool err_on_fail=true) const {
1675 
1676  // Old rank and new rank (computed later)
1677  size_t rank_old=this->rk;
1678  size_t rank_new=0;
1679 
1680  // Size of the new indices
1681  std::vector<size_t> size_new;
1682 
1683  // Reorganize the index specifications by both the old
1684  // and new indexing system
1685  std::vector<index_spec> spec_old(rank_old);
1686  std::vector<index_spec> spec_new;
1687 
1688  // Number of elements to sum over
1689  size_t n_sum_loop=1;
1690  // Number of interpolations to perform
1691  size_t n_interps=0;
1692 
1693  // Size of sums
1694  std::vector<size_t> sum_sizes;
1695 
1696  // List of indexes to interpolate
1697  std::vector<size_t> ix_to_interp;
1698 
1699  // New grid
1700  std::vector<double> new_grid;
1701 
1702  // Collect the statistics on the transformationand set the new grid
1703  for(size_t i=0;i<spec.size();i++) {
1704  if (spec[i].type==index_spec::index ||
1705  spec[i].type==index_spec::reverse) {
1706  if (spec[i].ix1>=rank_old) {
1707  if (err_on_fail) {
1708  O2SCL_ERR2("Index too large (index,reverse) in ",
1709  "tensor_grid::rearrange_and_copy().",
1711  } else {
1712  if (verbose>0) {
1713  std::cout << "Index too large (index,reverse) in "
1714  << "tensor_grid::rearrange_and_copy()."
1715  << std::endl;
1716  }
1717  return tensor_grid<>();
1718  }
1719  }
1720  size_new.push_back(this->size[spec[i].ix1]);
1721  // Use ix1 to store the destination index (which is
1722  // at this point equal to rank_new)
1723  spec_old[spec[i].ix1]=index_spec(spec[i].type,
1724  rank_new,
1725  spec[i].ix2,0,
1726  spec[i].val1);
1727  spec_new.push_back(index_spec(spec[i].type,
1728  spec[i].ix1,
1729  spec[i].ix2,0,
1730  spec[i].val1));
1731  rank_new++;
1732  // Update the new grid
1733  if (spec[i].type==index_spec::index) {
1734  for(size_t k=0;k<this->get_size(spec[i].ix1);k++) {
1735  new_grid.push_back(this->get_grid(spec[i].ix1,k));
1736  }
1737  } else {
1738  size_t nt=this->get_size(spec[i].ix1);
1739  for(size_t k=0;k<nt;k++) {
1740  new_grid.push_back(this->get_grid(spec[i].ix1,nt-1-k));
1741  }
1742  }
1743  } else if (spec[i].type==index_spec::range) {
1744  if (spec[i].ix1>=rank_old ||
1745  spec[i].ix2>=this->size[spec[i].ix1] ||
1746  spec[i].ix3>=this->size[spec[i].ix1]) {
1747  if (err_on_fail) {
1748  O2SCL_ERR2("Index too large (range) in ",
1749  "tensor::rearrange_and_copy().",o2scl::exc_einval);
1750  } else {
1751  if (verbose>0) {
1752  std::cout << "Index too large (range) in "
1753  << "tensor in tensor::rearrange_and_copy()."
1754  << std::endl;
1755  }
1756  return tensor_grid<>();
1757  }
1758  }
1759  if (spec[i].ix3>spec[i].ix2) {
1760  size_new.push_back(spec[i].ix3-spec[i].ix2+1);
1761  } else {
1762  size_new.push_back(spec[i].ix2-spec[i].ix3+1);
1763  }
1764  // Use ix1 to store the destination index (which is
1765  // at this point equal to rank_new)
1766  spec_old[spec[i].ix1]=
1767  index_spec(spec[i].type,rank_new,spec[i].ix2,
1768  spec[i].ix3,spec[i].val1);
1769  spec_new.push_back
1770  (index_spec(spec[i].type,spec[i].ix1,
1771  spec[i].ix2,spec[i].ix3,spec[i].val1));
1772  rank_new++;
1773  // Update the new grid
1774  if (spec[i].ix3>spec[i].ix2) {
1775  for(size_t k=0;k<spec[i].ix3-spec[i].ix2+1;k++) {
1776  new_grid.push_back(this->get_grid(spec[i].ix1,k+spec[i].ix2));
1777  }
1778  } else {
1779  for(size_t k=0;k<spec[i].ix2-spec[i].ix3+1;k++) {
1780  new_grid.push_back(this->get_grid(spec[i].ix1,k+spec[i].ix3));
1781  }
1782  }
1783  } else if (spec[i].type==index_spec::trace) {
1784  if (spec[i].ix1>=rank_old || spec[i].ix2>=rank_old) {
1785  if (err_on_fail) {
1786  O2SCL_ERR2("Index too large (trace) in ",
1787  "tensor_grid::rearrange_and_copy().",
1789  } else {
1790  if (verbose>0) {
1791  std::cout << "Index too large (trace) in "
1792  << "tensor_grid::rearrange_and_copy()."
1793  << std::endl;
1794  }
1795  return tensor_grid<>();
1796  }
1797  }
1798  if (this->size[spec[i].ix1]<this->size[spec[i].ix2]) {
1799  n_sum_loop*=this->size[spec[i].ix1];
1800  sum_sizes.push_back(this->size[spec[i].ix1]);
1801  } else {
1802  n_sum_loop*=this->size[spec[i].ix2];
1803  sum_sizes.push_back(this->size[spec[i].ix2]);
1804  }
1805  // We set the values of ix1 and ix2 so that ix2
1806  // always refers to the other index being traced over
1807  spec_old[spec[i].ix1]=index_spec(spec[i].type,spec[i].ix1,
1808  spec[i].ix2);
1809  spec_old[spec[i].ix2]=index_spec(spec[i].type,spec[i].ix2,
1810  spec[i].ix1);
1811  } else if (spec[i].type==index_spec::sum) {
1812  if (spec[i].ix1>=rank_old) {
1813  if (err_on_fail) {
1814  O2SCL_ERR2("Index too large (sum) in ",
1815  "tensor_grid::rearrange_and_copy().",
1817  } else {
1818  if (verbose>0) {
1819  std::cout << "Index " << spec[i].ix1
1820  << " too large (sum) in "
1821  << "tensor_grid::rearrange_and_copy()."
1822  << std::endl;
1823  }
1824  return tensor_grid<>();
1825  }
1826  }
1827  n_sum_loop*=this->size[spec[i].ix1];
1828  sum_sizes.push_back(this->size[spec[i].ix1]);
1829  spec_old[spec[i].ix1]=index_spec(spec[i].type,
1830  spec[i].ix1,
1831  spec[i].ix2,0,
1832  spec[i].val1);
1833  } else if (spec[i].type==index_spec::fixed) {
1834  if (spec[i].ix1>=rank_old ||
1835  spec[i].ix2>=this->size[spec[i].ix1]) {
1836  if (err_on_fail) {
1837  O2SCL_ERR2("Index too large (fixed) in ",
1838  "tensor_grid::rearrange_and_copy().",
1840  } else {
1841  if (verbose>0) {
1842  std::cout << "Index " << spec[i].ix1 << " or "
1843  << spec[i].ix2 << " too large (fixed) in "
1844  << "tensor_grid::rearrange_and_copy()."
1845  << std::endl;
1846  }
1847  return tensor_grid<>();
1848  }
1849  }
1850  // Use ix1 to store the destination index (which is
1851  // at this point equal to rank_new)
1852  spec_old[spec[i].ix1]=index_spec(spec[i].type,rank_new,
1853  spec[i].ix2);
1854  } else if (spec[i].type==index_spec::interp) {
1855  if (spec[i].ix1>=rank_old) {
1856  if (err_on_fail) {
1857  O2SCL_ERR2("Index too large (interp) in ",
1858  "tensor_grid::rearrange_and_copy().",
1860  } else {
1861  if (verbose>0) {
1862  std::cout << "Index " << spec[i].ix1
1863  << " too large (interp) in "
1864  << "tensor_grid::rearrange_and_copy()."
1865  << std::endl;
1866  }
1867  return tensor_grid<>();
1868  }
1869  }
1870  // Use ix1 to store the destination index (which is
1871  // at this point equal to rank_new)
1872  spec_old[spec[i].ix1]=index_spec(spec[i].type,
1873  rank_new,
1874  spec[i].ix2,0,
1875  spec[i].val1);
1876  n_interps++;
1877  ix_to_interp.push_back(spec[i].ix1);
1878  } else if (spec[i].type==index_spec::grid ||
1879  spec[i].type==index_spec::gridw) {
1880  if (spec[i].ix1>=rank_old) {
1881  if (err_on_fail) {
1882  O2SCL_ERR2("Index too large (grid) in ",
1883  "tensor_grid::rearrange_and_copy().",
1885  } else {
1886  if (verbose>0) {
1887  std::cout << "Index " << spec[i].ix1
1888  << " too large (grid) in "
1889  << "tensor_grid::rearrange_and_copy()."
1890  << std::endl;
1891  }
1892  return tensor_grid<>();
1893  }
1894  }
1895 
1896  // Setup new grid in log or linear mode and determine
1897  // spacing between grid points for later use
1898  if (spec[i].ix3==1) {
1899  // Log mode
1900  double rat;
1901  int curr_grid_size;
1902  if (spec[i].type==index_spec::grid) {
1903  rat=pow(spec[i].val2/spec[i].val1,1.0/((double)spec[i].ix2));
1904  spec[i].val3=rat;
1905  curr_grid_size=spec[i].ix2+1;
1906  } else {
1907  rat=spec[i].val3;
1908  curr_grid_size=((size_t)(log(spec[i].val2/spec[i].val1)/
1909  log(spec[i].val3)));
1910  }
1911  size_new.push_back(curr_grid_size);
1912  for(int k=0;k<curr_grid_size;k++) {
1913  if (k==0) {
1914  new_grid.push_back(spec[i].val1);
1915  } else if (curr_grid_size>1 && k==curr_grid_size-1) {
1916  new_grid.push_back(spec[i].val2);
1917  } else {
1918  new_grid.push_back(spec[i].val1*pow(rat,((double)k)));
1919  }
1920  }
1921  } else {
1922  // Linear mode
1923  double width;
1924  int curr_grid_size;
1925  if (spec[i].type==index_spec::grid) {
1926  width=(spec[i].val2-spec[i].val1)/((double)spec[i].ix2);
1927  spec[i].val3=width;
1928  curr_grid_size=spec[i].ix2+1;
1929  } else {
1930  width=spec[i].val3;
1931  curr_grid_size=((size_t)(1+(spec[i].val2-
1932  spec[i].val1)/spec[i].val3));
1933  }
1934  size_new.push_back(curr_grid_size);
1935  for(int k=0;k<curr_grid_size;k++) {
1936  if (k==0) {
1937  new_grid.push_back(spec[i].val1);
1938  } else if (curr_grid_size>1 && k==curr_grid_size-1) {
1939  new_grid.push_back(spec[i].val2);
1940  } else {
1941  new_grid.push_back(spec[i].val1+width*((double)k));
1942  }
1943  }
1944  }
1945 
1946  // Use ix1 to store the destination index (which is
1947  // at this point equal to rank_new)
1948  spec_old[spec[i].ix1]=index_spec
1949  (spec[i].type,rank_new,spec[i].ix2,spec[i].ix3,
1950  spec[i].val1,spec[i].val2,spec[i].val3);
1951  spec_new.push_back(index_spec
1952  (spec[i].type,spec[i].ix1,spec[i].ix2,
1953  spec[i].ix3,spec[i].val1,
1954  spec[i].val2,spec[i].val3));
1955 
1956  n_interps++;
1957  ix_to_interp.push_back(spec[i].ix1);
1958 
1959  rank_new++;
1960 
1961  } else {
1962  if (err_on_fail) {
1963  O2SCL_ERR2("Index specification type not allowed in ",
1964  "tensor_grid::rearrange_and_copy()",
1966  } else {
1967  if (verbose>0) {
1968  std::cout << "Index specification type not allowed in "
1969  << "tensor_grid::rearrange_and_copy()." << std::endl;
1970  }
1971  return tensor_grid<>();
1972  }
1973  }
1974  }
1975  size_t n_sums=sum_sizes.size();
1976 
1977  // Call the error handler if the input is invalid
1978  if (rank_new==0) {
1979  if (err_on_fail) {
1980  O2SCL_ERR2("Zero new indices in ",
1981  "tensor_grid::rearrange_and_copy()",
1983  } else {
1984  if (verbose>0) {
1985  std::cout << "Zero new indices in "
1986  << "tensor_grid::rearrange_and_copy()." << std::endl;
1987  }
1988  return tensor_grid<>();
1989  }
1990  }
1991  for(size_t i=0;i<rank_old;i++) {
1992  if (spec_old[i].type==index_spec::empty) {
1993  if (err_on_fail) {
1994  O2SCL_ERR2("Not all indices accounted for in ",
1995  "tensor_grid::rearrange_and_copy()",
1997  } else {
1998  if (verbose>0) {
1999  std::cout << "Index " << i << " not accounted for in "
2000  << "tensor_grid::rearrange_and_copy()." << std::endl;
2001  }
2002  return tensor_grid<>();
2003  }
2004  }
2005  }
2006 
2007  // Verbose output if necessary
2008  if (verbose>0) {
2009  std::cout << "Using a " << rank_old
2010  << " rank tensor to create a new "
2011  << rank_new << " rank tensor." << std::endl;
2012  }
2013  if (verbose>1) {
2014  for(size_t i=0;i<rank_old;i++) {
2015  std::cout << "Old index " << i;
2016  if (spec_old[i].type==index_spec::index) {
2017  std::cout << " is being remapped to new index "
2018  << spec_old[i].ix1 << "." << std::endl;
2019  } else if (spec_old[i].type==index_spec::range) {
2020  std::cout << " is being remapped to new index "
2021  << spec_old[i].ix1 << " with a range from "
2022  << spec_old[i].ix2 << " to " << spec_old[i].ix3
2023  << "." << std::endl;
2024  } else if (spec_old[i].type==index_spec::reverse) {
2025  std::cout << " is being reversed and remapped to new index "
2026  << spec_old[i].ix1 << "." << std::endl;
2027  } else if (spec_old[i].type==index_spec::trace) {
2028  std::cout << " is being traced with index "
2029  << spec_old[i].ix2 << "." << std::endl;
2030  } else if (spec_old[i].type==index_spec::sum) {
2031  std::cout << " is being summed." << std::endl;
2032  } else if (spec_old[i].type==index_spec::fixed) {
2033  std::cout << " is being fixed to " << spec_old[i].ix2
2034  << "." << std::endl;
2035  } else if (spec_old[i].type==index_spec::interp) {
2036  std::cout << " is being interpolated from value "
2037  << spec_old[i].val1 << "." << std::endl;
2038  } else if (spec_old[i].type==index_spec::grid ||
2039  spec_old[i].type==index_spec::gridw) {
2040  std::cout << " is being reinterpolated based on grid\n "
2041  << spec_old[i].val1 << " "
2042  << spec_old[i].val2 << " ";
2043  if (spec_old[i].type==index_spec::grid) {
2044  std::cout << spec_old[i].ix2;
2045  } else {
2046  std::cout << spec_old[i].val3;
2047  }
2048  if (spec_old[i].ix3==1) {
2049  std::cout << " (log)." << std::endl;
2050  } else {
2051  std::cout << "." << std::endl;
2052  }
2053  }
2054  }
2055  for(size_t i=0;i<rank_new;i++) {
2056  std::cout << "New index " << i;
2057  if (spec_new[i].type==index_spec::index) {
2058  std::cout << " was remapped from old index " << spec_new[i].ix1
2059  << "." << std::endl;
2060  } else if (spec_new[i].type==index_spec::range) {
2061  std::cout << " was remapped from old index " << spec_new[i].ix1
2062  << " using range from " << spec_new[i].ix2 << " to "
2063  << spec_new[i].ix3 << "." << std::endl;
2064  } else if (spec_new[i].type==index_spec::reverse) {
2065  std::cout << " was reversed and remapped from old index "
2066  << spec_new[i].ix1 << "." << std::endl;
2067  } else if (spec_new[i].type==index_spec::grid ||
2068  spec_new[i].type==index_spec::gridw) {
2069  std::cout << " was obtained from grid\n "
2070  << spec_new[i].val1 << " "
2071  << spec_new[i].val2 << " ";
2072  if (spec_new[i].type==index_spec::grid) {
2073  std::cout << spec_new[i].ix2;
2074  } else {
2075  std::cout << spec_new[i].val3;
2076  }
2077  if (spec_new[i].ix3==1) {
2078  std::cout << " (log)." << std::endl;
2079  } else {
2080  std::cout << "." << std::endl;
2081  }
2082  }
2083  }
2084  }
2085 
2086  // Create the new tensor object
2087  tensor_grid<> t_new(rank_new,size_new);
2088  t_new.set_grid_packed(new_grid);
2089  if (verbose>1) {
2090  std::cout << "New grid is: " << std::endl;
2091  for(size_t k=0;k<rank_new;k++) {
2092  std::cout << k << " (" << t_new.get_size(k) << "): ";
2093  if (t_new.get_size(k)>3) {
2094  std::cout << t_new.get_grid(k,0) << ", ";
2095  std::cout << t_new.get_grid(k,1) << " ... ";
2096  } else {
2097  for(int ell=0;ell<((int)t_new.get_size(k))-1;ell++) {
2098  std::cout << t_new.get_grid(k,ell) << ", ";
2099  }
2100  }
2101  std::cout << t_new.get_grid(k,t_new.get_size(k)-1) << std::endl;
2102  }
2103  std::cout << "Interpolations (" << n_interps << "): ";
2104  o2scl::vector_out(std::cout,ix_to_interp,true);
2105  }
2106 
2107  // Index arrays
2108  std::vector<size_t> ix_new(rank_new);
2109  std::vector<size_t> ix_old(rank_old);
2110  std::vector<size_t> sum_ix(n_sums);
2111 
2112  // Loop over the new tensor object
2113  for(size_t i=0;i<t_new.total_size();i++) {
2114 
2115  // Find the location in the new tensor object
2116  t_new.unpack_index(i,ix_new);
2117 
2118  // List of interpolated values
2119  std::vector<double> interp_vals;
2120 
2121  // Determine the location in the old tensor object
2122  for(size_t j=0;j<rank_old;j++) {
2123  if (spec_old[j].type==index_spec::index) {
2124  ix_old[j]=ix_new[spec_old[j].ix1];
2125  } else if (spec_old[j].type==index_spec::range) {
2126  if (spec_old[j].ix2<spec_old[j].ix3) {
2127  ix_old[j]=ix_new[spec_old[j].ix1]+spec_old[j].ix2;
2128  } else {
2129  ix_old[j]=spec_old[j].ix2-ix_new[spec_old[j].ix1];
2130  }
2131  } else if (spec_old[j].type==index_spec::reverse) {
2132  ix_old[j]=this->get_size(j)-1-ix_new[spec_old[j].ix1];
2133  } else if (spec_old[j].type==index_spec::fixed) {
2134  ix_old[j]=spec_old[j].ix2;
2135  } else if (spec_old[j].type==index_spec::interp) {
2136  interp_vals.push_back(spec_old[j].val1);
2137  } else if (spec_old[j].type==index_spec::grid) {
2138  if (spec_old[j].ix3==1) {
2139  double val=spec_old[j].val1*
2140  pow(spec_old[j].val3,ix_new[spec_old[j].ix1]);
2141  interp_vals.push_back(val);
2142  } else {
2143  double val=spec_old[j].val1+
2144  ix_new[spec_old[j].ix1]*spec_old[j].val3;
2145  interp_vals.push_back(val);
2146  }
2147  } else if (spec_old[j].type==index_spec::gridw) {
2148  if (spec_old[j].ix3==1) {
2149  double val=spec_old[j].val1*
2150  pow(spec_old[j].val3,ix_new[spec_old[j].ix1]);
2151  interp_vals.push_back(val);
2152  } else {
2153  double val=spec_old[j].val1+
2154  ix_new[spec_old[j].ix1]*spec_old[j].val3;
2155  interp_vals.push_back(val);
2156  }
2157  }
2158 
2159  }
2160 
2161  double val=0;
2162 
2163  for(size_t j=0;j<n_sum_loop;j++) {
2164 
2165  // This code is similar to tensor::unpack_index(), it unpacks
2166  // the index j to the indices which we are summing over.
2167  size_t j2=j, sub_size;
2168  for(size_t k=0;k<n_sums;k++) {
2169  if (k==n_sums-1) {
2170  sum_ix[k]=j2;
2171  } else {
2172  sub_size=1;
2173  for(size_t kk=k+1;kk<n_sums;kk++) sub_size*=sum_sizes[kk];
2174  sum_ix[k]=j2/sub_size;
2175  j2-=sub_size*(j2/sub_size);
2176  }
2177  }
2178 
2179  if (verbose>2) {
2180  std::cout << "n_sum_loop: " << n_sum_loop << " n_sums: "
2181  << n_sums << " sum_sizes: ";
2182  vector_out(std::cout,sum_sizes,true);
2183  std::cout << "j: " << j << " sum_ix: ";
2184  vector_out(std::cout,sum_ix,true);
2185  }
2186 
2187  // Remap from sum_ix to ix_old
2188  size_t cnt=0;
2189  for(size_t k=0;k<rank_old;k++) {
2190  if (spec_old[k].type==index_spec::sum) {
2191  if (cnt>=sum_ix.size()) {
2192  std::cout << "X: " << cnt << " " << sum_ix.size()
2193  << std::endl;
2194  O2SCL_ERR2("Bad sync 1 in sum_ix in ",
2195  "tensor_grid::rearrange_and_copy()",
2197  }
2198  ix_old[k]=sum_ix[cnt];
2199  cnt++;
2200  } else if (spec_old[k].type==index_spec::trace &&
2201  spec_old[k].ix1<spec_old[k].ix2) {
2202  if (cnt>=sum_ix.size()) {
2203  std::cout << "X: " << cnt << " " << sum_ix.size()
2204  << std::endl;
2205  O2SCL_ERR2("Bad sync 2 in sum_ix in ",
2206  "tensor_grid::rearrange_and_copy()",
2208  }
2209  ix_old[spec_old[k].ix1]=sum_ix[cnt];
2210  ix_old[spec_old[k].ix2]=sum_ix[cnt];
2211  cnt++;
2212  }
2213  }
2214 
2215  if (verbose>2) {
2216  std::cout << "Here old: ";
2217  vector_out(std::cout,ix_old,true);
2218  std::cout << "Here new: ";
2219  vector_out(std::cout,ix_new,true);
2220  }
2221  if (n_interps>0) {
2222  val+=this->interp_linear_partial(ix_to_interp,ix_old,interp_vals);
2223  } else {
2224  val+=this->get(ix_old);
2225  }
2226 
2227  }
2228 
2229  // Set the new point by performing the linear interpolation
2230  t_new.set(ix_new,val);
2231  }
2232 
2233  return t_new;
2234  }
2235 
2236  };
2237 
2238  /** \brief Rank 1 tensor with a grid
2239 
2240  \future Make rank-specific get_val and set_val functions?
2241  */
2242  template<class vec_t=std::vector<double>,
2243  class vec_size_t=std::vector<size_t> > class tensor_grid1 :
2244  public tensor_grid<vec_t,vec_size_t> {
2245 
2246  public:
2247 
2248  /// Create an empty tensor
2249  tensor_grid1() : tensor_grid<vec_t,vec_size_t>() {}
2250 
2251  /// Create a rank 2 tensor of size \c (sz,sz2,sz3)
2252  tensor_grid1(size_t sz) : tensor_grid<vec_t,vec_size_t>() {
2253  this->rk=1;
2254  this->size.resize(1);
2255  this->size[0]=sz;
2256  this->data.resize(sz);
2257  this->grid_set=false;
2258  }
2259 
2260 #ifdef O2SCL_NEVER_DEFINED
2261  }{
2262 #endif
2263 
2264  virtual ~tensor_grid1() {
2265  }
2266 
2267  /// Get the element indexed by \c (ix1)
2268  double &get(size_t ix1) {
2269  size_t sz[1]={ix1};
2271  }
2272 
2273  /// Get the element indexed by \c (ix1)
2274  const double &get(size_t ix1) const {
2275  size_t sz[1]={ix1};
2277  }
2278 
2279  /// Set the element indexed by \c (ix1) to value \c val
2280  void set(size_t ix1, double val) {
2281  size_t sz[1]={ix1};
2283  }
2284 
2285  /// Interpolate \c x and return the results
2286  template<class range_t=ub_range, class data_range_t=ubvector_range,
2287  class index_range_t=ubvector_size_t_range>
2288  double interp(double x) {
2289  return tensor_grid<vec_t,vec_size_t>::template interpolate
2290  <range_t,data_range_t,index_range_t>(&x);
2291  }
2292 
2293  /// Interpolate \c x and return the results
2294  double interp_linear(double x) {
2295  double arr[1]={x};
2297  }
2298 
2299  };
2300 
2301  /** \brief Rank 2 tensor with a grid
2302  */
2303  template<class vec_t=std::vector<double>,
2304  class vec_size_t=std::vector<size_t> > class tensor_grid2 :
2305  public tensor_grid<vec_t,vec_size_t> {
2306 
2307  public:
2308 
2309  /// Create an empty tensor
2310  tensor_grid2() : tensor_grid<vec_t,vec_size_t>() {}
2311 
2312  /// Create a rank 2 tensor of size \c (sz,sz2)
2313  tensor_grid2(size_t sz, size_t sz2) : tensor_grid<vec_t,vec_size_t>() {
2314  this->rk=2;
2315  this->size.resize(2);
2316  this->size[0]=sz;
2317  this->size[1]=sz2;
2318  size_t tot=sz*sz2;
2319  this->data.resize(tot);
2320  this->grid_set=false;
2321  }
2322 
2323 #ifdef O2SCL_NEVER_DEFINED
2324  }{
2325 #endif
2326 
2327  virtual ~tensor_grid2() {
2328  }
2329 
2330  /// Get the element indexed by \c (ix1,ix2)
2331  double &get(size_t ix1, size_t ix2) {
2332  size_t sz[2]={ix1,ix2};
2334  }
2335 
2336  /// Get the element indexed by \c (ix1,ix2)
2337  const double &get(size_t ix1, size_t ix2) const {
2338  size_t sz[2]={ix1,ix2};
2340  }
2341 
2342  /// Set the element indexed by \c (ix1,ix2) to value \c val
2343  void set(size_t ix1, size_t ix2, double val) {
2344  size_t sz[2]={ix1,ix2};
2346  return;
2347  }
2348 
2349  /// Interpolate \c (x,y) and return the results
2350  template<class range_t=ub_range, class data_range_t=ubvector_range,
2351  class index_range_t=ubvector_size_t_range>
2352  double interp(double x, double y) {
2353  double arr[2]={x,y};
2354  return tensor_grid<vec_t,vec_size_t>::template interpolate
2355  <range_t,data_range_t,index_range_t>(arr);
2356  }
2357 
2358  /// Interpolate \c (x,y) and return the results
2359  double interp_linear(double x, double y) {
2360  double arr[2]={x,y};
2362  }
2363 
2364  };
2365 
2366  /** \brief Rank 3 tensor with a grid
2367  */
2368  template<class vec_t=std::vector<double>,
2369  class vec_size_t=std::vector<size_t> > class tensor_grid3 :
2370  public tensor_grid<vec_t,vec_size_t> {
2371 
2372  public:
2373 
2374  /// Create an empty tensor
2375  tensor_grid3() : tensor_grid<vec_t,vec_size_t>() {}
2376 
2377  /// Create a rank 3 tensor of size \c (sz,sz2,sz3)
2378  tensor_grid3(size_t sz, size_t sz2, size_t sz3) :
2379  tensor_grid<vec_t,vec_size_t>() {
2380  this->rk=3;
2381  this->size.resize(3);
2382  this->size[0]=sz;
2383  this->size[1]=sz2;
2384  this->size[2]=sz3;
2385  size_t tot=sz*sz2*sz3;
2386  this->data.resize(tot);
2387  this->grid_set=false;
2388  }
2389 
2390 #ifdef O2SCL_NEVER_DEFINED
2391  }{
2392 #endif
2393 
2394  virtual ~tensor_grid3() {
2395  }
2396 
2397  /// Get the element indexed by \c (ix1,ix2,ix3)
2398  double &get(size_t ix1, size_t ix2, size_t ix3) {
2399  size_t sz[3]={ix1,ix2,ix3};
2401  }
2402 
2403  /// Get the element indexed by \c (ix1,ix2,ix3)
2404  const double &get(size_t ix1, size_t ix2, size_t ix3) const {
2405  size_t sz[3]={ix1,ix2,ix3};
2407  }
2408 
2409  /// Set the element indexed by \c (ix1,ix2,ix3) to value \c val
2410  void set(size_t ix1, size_t ix2, size_t ix3, double val) {
2411  size_t sz[3]={ix1,ix2, ix3};
2413  return;
2414  }
2415 
2416  /// Interpolate \c (x,y,z) and return the results
2417  template<class range_t=ub_range, class data_range_t=ubvector_range,
2418  class index_range_t=ubvector_size_t_range>
2419  double interp(double x, double y, double z) {
2420  double arr[3]={x,y,z};
2421  return tensor_grid<vec_t,vec_size_t>::template interpolate
2422  <range_t,data_range_t,index_range_t>(arr);
2423  }
2424 
2425  /// Interpolate \c (x,y,z) and return the results
2426  double interp_linear(double x, double y, double z) {
2427  double arr[3]={x,y,z};
2429  }
2430 
2431  };
2432 
2433  /** \brief Rank 4 tensor with a grid
2434  */
2435  template<class vec_t=std::vector<double>,
2436  class vec_size_t=std::vector<size_t> > class tensor_grid4 :
2437  public tensor_grid<vec_t,vec_size_t> {
2438 
2439  public:
2440 
2441  /// Create an empty tensor
2442  tensor_grid4() : tensor_grid<vec_t,vec_size_t>() {}
2443 
2444  /// Create a rank 4 tensor of size \c (sz,sz2,sz3,sz4)
2445  tensor_grid4(size_t sz, size_t sz2, size_t sz3,
2446  size_t sz4) : tensor_grid<vec_t,vec_size_t>() {
2447  this->rk=4;
2448  this->size.resize(4);
2449  this->size[0]=sz;
2450  this->size[1]=sz2;
2451  this->size[2]=sz3;
2452  this->size[3]=sz4;
2453  size_t tot=sz*sz2*sz3*sz4;
2454  this->data.resize(tot);
2455  this->grid_set=false;
2456  }
2457 
2458 #ifdef O2SCL_NEVER_DEFINED
2459  }{
2460 #endif
2461 
2462  virtual ~tensor_grid4() {
2463  }
2464 
2465  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
2466  double &get(size_t ix1, size_t ix2, size_t ix3, size_t ix4) {
2467  size_t sz[4]={ix1,ix2,ix3,ix4};
2469  }
2470 
2471  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
2472  const double &get(size_t ix1, size_t ix2, size_t ix3,
2473  size_t ix4) const {
2474  size_t sz[4]={ix1,ix2,ix3,ix4};
2476  }
2477 
2478  /// Set the element indexed by \c (ix1,ix2,ix3,ix4) to value \c val
2479  void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4,
2480  double val) {
2481  size_t sz[4]={ix1,ix2,ix3,ix4};
2483  return;
2484  }
2485 
2486  /// Interpolate \c (x,y,z,a) and return the results
2487  template<class range_t=ub_range, class data_range_t=ubvector_range,
2488  class index_range_t=ubvector_size_t_range>
2489  double interp(double x, double y, double z, double a) {
2490  double arr[4]={x,y,z,a};
2491  return tensor_grid<vec_t,vec_size_t>::template interpolate
2492  <range_t,data_range_t,index_range_t>(arr);
2493  }
2494 
2495  /// Interpolate \c (x,y,z,a) and return the results
2496  double interp_linear(double x, double y, double z, double a) {
2497  double arr[4]={x,y,z,a};
2499  }
2500 
2501  };
2502 
2503 #ifndef DOXYGEN_NO_O2NS
2504 }
2505 #endif
2506 
2507 #endif
2508 
2509 
2510 
o2scl::tensor
Tensor class with arbitrary dimensions.
Definition: tensor.h:226
o2scl::tensor_grid4::interp_linear
double interp_linear(double x, double y, double z, double a)
Interpolate (x,y,z,a) and return the results.
Definition: tensor_grid.h:2496
o2scl::index_spec::index
static const size_t index
Retain an index.
Definition: tensor.h:78
o2scl::tensor_grid::set_interp_type
void set_interp_type(size_t interp_type)
Set interpolation type for interpolate()
Definition: tensor_grid.h:1157
o2scl::tensor_grid::interp_linear_vec0
void interp_linear_vec0(vec2_t &v, vec3_t &res) const
Perform a linear interpolation of v[1] to v[n-1] resulting in a vector.
Definition: tensor_grid.h:1460
o2scl::tensor_grid2
Rank 2 tensor with a grid.
Definition: tensor_grid.h:2304
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.
boost::numeric::ublas::vector< double >
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::tensor::is_valid
void is_valid() const
Check that the o2scl::tensor object is valid.
Definition: tensor.h:291
o2scl::tensor_grid::copy_table3d_interp
void copy_table3d_interp(size_t ix_x, size_t ix_y, size_vec2_t &index, table3d &tab, std::string slice_name="z") const
Copy to a slice in a table3d object using interpolation.
Definition: tensor_grid.h:1012
o2scl::tensor_grid::lookup_grid_val
size_t lookup_grid_val(size_t i, const double &val, double &val2) const
Lookup index for grid closest to val, returning the grid point.
Definition: tensor_grid.h:674
o2scl::tensor_grid4::get
const double & 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_grid.h:2472
o2scl::tensor_grid::lookup_grid
size_t lookup_grid(size_t i, double val) const
Lookup index for grid closest to val.
Definition: tensor_grid.h:706
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::tensor_grid3
Rank 3 tensor with a grid.
Definition: tensor_grid.h:2369
o2scl::tensor_grid::get_val
double get_val(const vec2_t &gridp)
Get the element closest to grid point gridp.
Definition: tensor_grid.h:360
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::table3d::new_slice
void new_slice(std::string name)
Add a new slice.
o2scl::tensor_grid::is_valid
void is_valid() const
Check that the o2scl::tensor_grid object is valid.
Definition: tensor_grid.h:245
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_grid::tensor_grid
tensor_grid(std::vector< uniform_grid< double > > &ugs)
Create a tensor with a grid defined by a set of o2scl::uniform_grid objects.
Definition: tensor_grid.h:222
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::index_spec
Index specification.
Definition: tensor.h:54
o2scl::tensor_grid3::tensor_grid3
tensor_grid3(size_t sz, size_t sz2, size_t sz3)
Create a rank 3 tensor of size (sz,sz2,sz3)
Definition: tensor_grid.h:2378
o2scl::tensor_grid::get_grid
const vec_t & get_grid() const
Get grid.
Definition: tensor_grid.h:406
o2scl::tensor_grid3::get
double & get(size_t ix1, size_t ix2, size_t ix3)
Get the element indexed by (ix1,ix2,ix3)
Definition: tensor_grid.h:2398
o2scl::tensor_grid2::get
const double & get(size_t ix1, size_t ix2) const
Get the element indexed by (ix1,ix2)
Definition: tensor_grid.h:2337
o2scl::tensor_grid1::interp_linear
double interp_linear(double x)
Interpolate x and return the results.
Definition: tensor_grid.h:2294
o2scl::tensor_grid1::get
const double & get(size_t ix1) const
Get the element indexed by (ix1)
Definition: tensor_grid.h:2274
o2scl::tensor::size
vec_size_t size
A rank-sized vector of the sizes of each dimension.
Definition: tensor.h:238
o2scl::tensor_grid::resize
void resize(size_t rank, const size_vec2_t &dim)
Resize the tensor to rank rank with sizes given in dim.
Definition: tensor_grid.h:424
o2scl::tensor_grid::lookup_grid_packed
size_t lookup_grid_packed(size_t i, double val) const
Lookup internal packed grid index for point closest to val.
Definition: tensor_grid.h:766
o2scl::tensor_grid::tensor_grid
tensor_grid()
Create an empty tensor with zero rank.
Definition: tensor_grid.h:192
o2scl::tensor_grid2::get
double & get(size_t ix1, size_t ix2)
Get the element indexed by (ix1,ix2)
Definition: tensor_grid.h:2331
o2scl::tensor_grid::clear
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor_grid.h:1146
o2scl::itp_linear
@ itp_linear
Linear.
Definition: interp.h:71
o2scl::interp_linear
Linear interpolation (GSL)
Definition: interp.h:210
o2scl::tensor_grid2::interp
double interp(double x, double y)
Interpolate (x,y) and return the results.
Definition: tensor_grid.h:2352
o2scl_hdf::hdf_input
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
Input a o2scl::table object from a hdf_file.
Definition: hdf_io.h:148
o2scl::tensor_grid::get_val
double get_val(const vec2_t &gridp, vec3_t &closest)
Get the element closest to grid point gridp, store grid values in closest and return value.
Definition: tensor_grid.h:386
o2scl::tensor_grid4
Rank 4 tensor with a grid.
Definition: tensor_grid.h:2436
o2scl::index_spec::trace
static const size_t trace
Perform a trace (sum over two indices)
Definition: tensor.h:84
o2scl::tensor_grid::copy_table3d_interp_values
void copy_table3d_interp_values(size_t ix_x, size_t ix_y, vec2_t &values, table3d &tab, std::string slice_name="z", int verbose=0) const
Copy to a slice in a table3d object using interpolation.
Definition: tensor_grid.h:1055
o2scl::interp_vec
Interpolation class for pre-specified vector.
Definition: interp.h:1795
o2scl::tensor_grid4::get
double & get(size_t ix1, size_t ix2, size_t ix3, size_t ix4)
Get the element indexed by (ix1,ix2,ix3,ix4)
Definition: tensor_grid.h:2466
o2scl::calculator::compile
void compile(const char *expr, std::map< std::string, double > *vars=0, bool debug=false, std::map< std::string, int > opPrec=opPrecedence)
Compile expression expr using variables specified in vars and return an integer to indicate success o...
o2scl::tensor_grid2::tensor_grid2
tensor_grid2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
Definition: tensor_grid.h:2313
o2scl::tensor_grid::set_grid
void set_grid(size_t i, size_t j, double val)
Set the jth value on the ith grid.
Definition: tensor_grid.h:651
o2scl::tensor_grid1::interp
double interp(double x)
Interpolate x and return the results.
Definition: tensor_grid.h:2288
o2scl::tensor_grid::grid_set
bool grid_set
If true, the grid has been set by the user.
Definition: tensor_grid.h:180
o2scl::tensor_grid::get_grid
double get_grid(size_t i, size_t j) const
Lookup jth value on the ith grid.
Definition: tensor_grid.h:634
o2scl::tensor_grid2::tensor_grid2
tensor_grid2()
Create an empty tensor.
Definition: tensor_grid.h:2310
o2scl::tensor_grid4::interp
double interp(double x, double y, double z, double a)
Interpolate (x,y,z,a) and return the results.
Definition: tensor_grid.h:2489
o2scl::tensor_grid::lookup_grid_packed_val
size_t lookup_grid_packed_val(size_t i, double val, double &val2) const
Lookup internal packed grid index for point closest to val and store closest value in val2.
Definition: tensor_grid.h:737
o2scl::tensor_grid::set_grid_i_func
void set_grid_i_func(size_t ix, std::string func)
Set grid for one index from a function.
Definition: tensor_grid.h:567
o2scl::tensor_grid::grid
vec_t grid
A rank-sized set of arrays for the grid points.
Definition: tensor_grid.h:177
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::ub_range
boost::numeric::ublas::range ub_range
A ublas::range typedef for o2scl::tensor_grid and related classes in src/base/tensor_grid....
Definition: tensor_grid.h:67
o2scl::calculator
Evaluate a mathematical expression in a string.
Definition: shunting_yard.h:119
o2scl::tensor_grid3::interp
double interp(double x, double y, double z)
Interpolate (x,y,z) and return the results.
Definition: tensor_grid.h:2419
o2scl::tensor_grid::rearrange_and_copy
tensor_grid rearrange_and_copy(std::vector< index_spec > spec, int verbose=0, bool err_on_fail=true) const
Rearrange, sum and copy current tensor to a new tensor.
Definition: tensor_grid.h:1672
o2scl::ubvector_size_t_range
boost::numeric::ublas::vector_range< boost::numeric::ublas::vector< size_t > > ubvector_size_t_range
A ublas::vector_range typedef (size_t version) for o2scl::tensor_grid and related classes in src/base...
Definition: tensor_grid.h:80
o2scl::tensor_grid::~tensor_grid
virtual ~tensor_grid()
Destructor.
Definition: tensor_grid.h:237
o2scl::tensor_grid::interp_linear_power_two_vec0
void interp_linear_power_two_vec0(vec2_t &v, vec3_t &res) const
Perform linear interpolation assuming that the last n-1 indices can take only two values.
Definition: tensor_grid.h:1521
o2scl::index_spec::interp
static const size_t interp
Interpolate a value to fix an index.
Definition: tensor.h:90
o2scl_hdf
The O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl namespace ...
Definition: table.h:59
o2scl::tensor_grid
Tensor class with arbitrary dimensions with a grid.
Definition: tensor_grid.h:46
o2scl::table3d::get
double & get(size_t ix, size_t iy, std::string name)
Get element in slice name at location ix,iy
o2scl::table3d::get_grid_x
double get_grid_x(size_t ix) const
Get x grid point at index ix.
o2scl::tensor_grid3::tensor_grid3
tensor_grid3()
Create an empty tensor.
Definition: tensor_grid.h:2375
o2scl::tensor_grid::itype
size_t itype
Interpolation type.
Definition: tensor_grid.h:183
o2scl::tensor_grid::set_grid_packed
void set_grid_packed(const vec2_t &grid_vec)
Set the grid.
Definition: tensor_grid.h:486
o2scl::uniform_grid< double >
o2scl::tensor_grid3::interp_linear
double interp_linear(double x, double y, double z)
Interpolate (x,y,z) and return the results.
Definition: tensor_grid.h:2426
o2scl::tensor::clear
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor.h:348
o2scl::search_vec::find
size_t find(const double x0)
Search an increasing or decreasing vector for the interval containing x0
Definition: search_vec.h:123
o2scl::search_vec
Searching class for monotonic data with caching.
Definition: search_vec.h:74
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::table3d::is_size_set
bool is_size_set() const
True if the size of the table has been set.
o2scl::table3d::is_slice
bool is_slice(std::string name, size_t &ix) const
Return true if slice is already present.
o2scl::tensor_grid::interp_linear
double interp_linear(vec2_t &v) const
Perform a linear interpolation of v into the function implied by the tensor and grid.
Definition: tensor_grid.h:1366
o2scl::tensor_grid::set_val
void set_val(const vec2_t &grdp, double val, vec3_t &closest)
Set the element closest to grid point grdp to value val.
Definition: tensor_grid.h:335
o2scl::tensor_grid::interpolate
double interpolate(double *vals)
Interpolate values vals into the tensor, returning the result.
Definition: tensor_grid.h:1191
o2scl::index_spec::range
static const size_t range
Choose a new range for an index.
Definition: tensor.h:88
o2scl::table3d::is_xy_set
bool is_xy_set() const
True if the grid has been set.
o2scl::tensor_grid::copy_table3d_align
void copy_table3d_align(size_t ix_x, size_t ix_y, size_vec2_t &index, table3d &tab, std::string slice_name="z") const
Create a slice in a o2scl::table3d object with an aligned grid.
Definition: tensor_grid.h:921
o2scl::tensor_grid::copy_grid
void copy_grid(size_t i, rvec_t &v) const
Copy grid for index i to vector v.
Definition: tensor_grid.h:623
o2scl::tensor_grid1::get
double & get(size_t ix1)
Get the element indexed by (ix1)
Definition: tensor_grid.h:2268
o2scl_hdf::hdf_file
Store data in an O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$sc...
Definition: hdf_file.h:105
o2scl::vector_range
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:3221
o2scl::exc_esanity
@ exc_esanity
sanity check failed - shouldn't happen
Definition: err_hnd.h:65
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_grid::set_grid
void set_grid(std::vector< uniform_grid< double > > &ugs)
Set grid from a vector of uniform grid objects.
Definition: tensor_grid.h:599
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::tensor_grid::interp_linear_vec
void interp_linear_vec(vec2_t &v, size_t ifree, vec3_t &res) const
Perform a linear interpolation of v into the tensor leaving one index free resulting in a vector.
Definition: tensor_grid.h:1575
o2scl::table3d::get_size
void get_size(size_t &nx, size_t &ny) const
Get the size of the slices.
o2scl::tensor_grid::interp_linear_power_two
double interp_linear_power_two(vec2_t &v) const
Perform linear interpolation assuming that all indices can take only two values.
Definition: tensor_grid.h:1417
o2scl::tensor_grid::set_val
void set_val(const vec2_t &grdp, double val)
Set the element closest to grid point grdp to value val.
Definition: tensor_grid.h:306
o2scl::table3d::get_grid_y
double get_grid_y(size_t iy) const
Get y grid point at index iy.
o2scl::calculator::eval
double eval(std::map< std::string, double > *vars=0)
Evalate the previously compiled expression using variables specified in vars.
o2scl::tensor_grid::is_grid_set
bool is_grid_set() const
Return true if the grid has been set.
Definition: tensor_grid.h:462
o2scl::ubvector_range
boost::numeric::ublas::vector_range< boost::numeric::ublas::vector< double > > ubvector_range
A ublas::vector_range typedef for o2scl::tensor_grid and related classes in src/base/tensor_grid....
Definition: tensor_grid.h:73
o2scl::tensor_grid1::set
void set(size_t ix1, double val)
Set the element indexed by (ix1) to value val.
Definition: tensor_grid.h:2280
o2scl::tensor_grid1
Rank 1 tensor with a grid.
Definition: tensor_grid.h:2243
o2scl::tensor_grid::lookup_grid_vec
void lookup_grid_vec(const vec2_t &vals, size_vec2_t &indices) const
Lookup indices for grid closest point to vals.
Definition: tensor_grid.h:722
o2scl::tensor_grid::set_grid_i_vec
void set_grid_i_vec(size_t ix, const vec2_t &grid_vec)
Set grid for one index from a vector.
Definition: tensor_grid.h:543
o2scl::tensor_grid3::set
void set(size_t ix1, size_t ix2, size_t ix3, double val)
Set the element indexed by (ix1,ix2,ix3) to value val.
Definition: tensor_grid.h:2410
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_grid4::tensor_grid4
tensor_grid4()
Create an empty tensor.
Definition: tensor_grid.h:2442
o2scl::tensor_grid1::tensor_grid1
tensor_grid1(size_t sz)
Create a rank 2 tensor of size (sz,sz2,sz3)
Definition: tensor_grid.h:2252
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_grid::set_grid
void set_grid(const vec_vec_t &grid_vecs)
Set grid from a vector of vectors of grid points.
Definition: tensor_grid.h:504
o2scl::tensor_grid1::tensor_grid1
tensor_grid1()
Create an empty tensor.
Definition: tensor_grid.h:2249
o2scl::szttos
std::string szttos(size_t x)
Convert a size_t to a string.
o2scl::tensor_grid4::set
void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4, double val)
Set the element indexed by (ix1,ix2,ix3,ix4) to value val.
Definition: tensor_grid.h:2479
o2scl::index_spec::empty
static const size_t empty
Empty specification.
Definition: tensor.h:76
o2scl::tensor_grid3::get
const double & get(size_t ix1, size_t ix2, size_t ix3) const
Get the element indexed by (ix1,ix2,ix3)
Definition: tensor_grid.h:2404
o2scl::tensor_grid2::interp_linear
double interp_linear(double x, double y)
Interpolate (x,y) and return the results.
Definition: tensor_grid.h:2359
o2scl::tensor_grid2::set
void set(size_t ix1, size_t ix2, double val)
Set the element indexed by (ix1,ix2) to value val.
Definition: tensor_grid.h:2343
o2scl::tensor_grid::default_grid
void default_grid()
Use a default grid which just uses the index.
Definition: tensor_grid.h:525
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_grid::copy_slice_interp
tensor_grid copy_slice_interp(size_vec2_t &ifix, vec2_t &vals) const
Copy an abitrary slice by fixing 1 or more indices and use interpolation to return a new tensor_grid ...
Definition: tensor_grid.h:778
o2scl::tensor_grid4::tensor_grid4
tensor_grid4(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_grid.h:2445
o2scl::tensor_grid::tensor_grid
tensor_grid(size_t rank, const size_vec_t &dim)
Create a tensor of rank rank with sizes given in dim.
Definition: tensor_grid.h:204

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