tensor_grid.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_TENSOR_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  typedef boost::numeric::ublas::range ub_range;
66  <boost::numeric::ublas::vector<double> > ubvector_range;
68  <boost::numeric::ublas::vector<size_t> > ubvector_size_t_range;
69 
70  /** \brief Tensor class with arbitrary dimensions with a grid
71 
72  This tensor class allows one to assign the indexes to numerical
73  scales, effectively defining a data set on an n-dimensional
74  grid. To set the grid, use \ref default_grid(), \ref set_grid()
75  or \ref set_grid_packed().
76 
77  By convention, member functions ending in the <tt>_val</tt>
78  suffix return the closest grid-point to some user-specified
79  values.
80 
81  \comment
82  I like how hist_new only allows you to set the
83  grid if it's the same size as the old grid or the tensor
84  is empty. This same practice could be applied here, to
85  force the user to clear the tensor before resetting the grid.
86  (10/24/11 - Actually, maybe this is a bad idea, because
87  this class is more analogous to ubvector, which doesn't
88  behave this way.)
89  \endcomment
90 
91  \b Slicing
92 
93  New \ref o2scl::tensor_grid objects can be obtained
94  by fixing any set of indices using \ref copy_slice_interp().
95 
96  Fixing all but two indices also results in a \ref o2scl::table3d
97  object, and five functions perform this task in different ways.
98  The function \ref copy_table3d_align() copies a two-dimensional
99  slice to a \ref o2scl::table3d object presuming that the grid in
100  the \ref o2scl::table3d object has already been set and exactly
101  matches the corresponding sizes for the selected tensor indices.
102  This function does not check that the grids between the two
103  objects match, it only ensures that they have the same size. In
104  order to copy to a \ref o2scl::table3d object and set its grid
105  to match that from the unfixed indices in the \ref
106  o2scl::tensor_grid object, the function \ref
107  copy_table3d_align_setxy() can be used. The function \ref
108  copy_table3d_interp() uses interpolation to extract values from
109  the \ref o2scl::tensor_grid object. It allows the user to select
110  indices to be fixed and then uses the values in the grid in the
111  \ref o2scl::table3d object for the indices which vary.
112  Alternatively \ref copy_table3d_interp_values() allows the user
113  to specify values on the grid for the indices to be fixed and
114  uses the grid in the \ref o2scl::table3d object for the indices
115  which vary. Finally, \ref copy_table3d_interp_values_setxy() acts
116  like \ref copy_table3d_interp_values() except that it sets the
117  \ref o2scl::table3d grid to be the same as the grid in the \ref
118  o2scl::tensor_grid object which corresponds to the indices which
119  are being varied.
120 
121  \b Notes and Todos
122 
123  \note Currently, HDF5 I/O is only allowed if the tensor is
124  allocated with std::vector-based types, and the \ref
125  interpolate() function only works with ublas-based vector types.
126 
127  \todo It is possible for the user to create a tensor_grid
128  object, upcast it to a tensor object, and then use
129  tensor::resize() to resize the tensor, failing to resize the
130  grid. This can be fixed by ensuring that resize functions
131  are virtual and have a version in tensor_grid which ensure
132  that the grid and tensor data are matched.
133 
134  \future Is it really necessary that get_data() is public and not
135  const? This is used in HDF5 I/O, but the HDF5 output function
136  already appears to be a friend and it conflicts with the parent
137  version from tensor which is non-const. Also, there are some
138  applications for which it's really helpful to have access to the
139  base data object. I'm not sure what to do here.
140 
141  \future A swap function for the data object might be helpful.
142 
143  \future Only allocate space for grid if it is set.
144 
145  \future As with \ref o2scl::tensor, generalize to other
146  base data types.
147  */
148  template<class vec_t=std::vector<double>,
149  class vec_size_t=std::vector<size_t> > class tensor_grid :
150  public tensor<double,vec_t,vec_size_t> {
151 
152  public:
153 
154 #ifndef DOXYGEN_INTERNAL
155 
156  protected:
157 
158 #ifdef O2SCL_NEVER_DEFINED
159  }{
160 #endif
161 
162  /// A rank-sized set of arrays for the grid points
163  vec_t grid;
164 
165  /// If true, the grid has been set by the user
166  bool grid_set;
167 
168  /// Interpolation type
169  size_t itype;
170 
171 #endif
172 
173  public:
174 
175  /// \name Constructors and Destructors
176  //@{
177  /// Create an empty tensor with zero rank
178  tensor_grid() : tensor<double,vec_t,vec_size_t>() {
179  grid_set=false;
180  itype=itp_linear;
181  }
182 
183  /** \brief Create a tensor of rank \c rank with sizes given in \c dim
184 
185  The parameter \c dim must be a vector of sizes with length \c
186  rank. If the user requests any of the sizes to be zero, this
187  constructor will call the error handler.
188  */
189  template<class size_vec_t>
190  tensor_grid(size_t rank, const size_vec_t &dim) :
191  tensor<double,vec_t,vec_size_t>(rank,dim) {
192  grid_set=false;
193  itype=itp_linear;
194  // Note that the parent method sets rk to be equal to rank
195  for(size_t i=0;i<this->rk;i++) {
196  if (dim[i]==0) {
197  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
198  "rank for index "+szttos(i)+" in tensor_grid::"+
199  "tensor_grid(size_t,size_vec_t)").c_str(),
200  exc_einval);
201  }
202  }
203  }
204 
205  /** \brief Create a tensor with a grid defined by a set
206  of \ref o2scl::uniform_grid objects
207  */
208  tensor_grid(std::vector<uniform_grid<double> > &ugs) :
209  tensor<double,vec_t,vec_size_t>() {
210  this->rk=ugs.size();
211  itype=itp_linear;
212  size_t tot=1;
213  for(size_t j=0;j<this->rk;j++) {
214  this->size.push_back(ugs[j].get_npoints());
215  tot*=ugs[j].get_npoints();
216  }
217  this->data.resize(tot);
218  set_grid(ugs);
219  }
220 
221  /** \brief Destructor
222  */
223  virtual ~tensor_grid() {
224  }
225  //@}
226 
227  /// \name Method to check for valid object
228  //@{
229  /** \brief Check that the \ref o2scl::tensor_grid object is valid
230  */
231  void is_valid() const {
232 
234 
235  if (this->rk>0 && grid_set) {
236  size_t tot2=0;
237  for(size_t i=0;i<this->rk;i++) {
238  tot2+=this->size[i];
239  }
240 
241  if (tot2!=grid.size()) {
242  O2SCL_ERR2("Value grid_set is true but grid vector size ",
243  "is wrong in tensor_grid::is_valid().",
245  }
246  }
247 
248  if (!grid_set && grid.size()>0) {
249  O2SCL_ERR2("Value grid_set is false but grid vector size ",
250  "is not zero in tensor_grid::is_valid().",
252  }
253 
254  return;
255  }
256  //@}
257 
258  /// \name Copy constructors
259  //@{
260  /** \brief Copy using <tt>operator()</tt>
261  */
264  this->rk=t.rk;
265  this->data=t.data;
266  this->size=t.size;
267  grid=t.grid;
268  grid_set=t.grid_set;
269  itype=t.itype;
270  }
271 
272  /** \brief Copy using <tt>operator=()</tt>
273  */
276  if (this!=&t) {
277  this->rk=t.rk;
278  this->data=t.data;
279  this->size=t.size;
280  grid=t.grid;
281  grid_set=t.grid_set;
282  itype=t.itype;
283  }
284  return *this;
285  }
286  //@}
287 
288  /// \name Set functions
289  //@{
290  /// Set the element closest to grid point \c grdp to value \c val
291  template<class vec2_t>
292  void set_val(const vec2_t &grdp, double val) {
293 
294  // Find indices
295  vec_size_t index(this->rk);
296  for(size_t i=0;i<this->rk;i++) {
297  index[i]=lookup_grid(i,grdp[i]);
298  }
299 
300  // Pack
301  size_t ix=index[0];
302  for(size_t i=1;i<this->rk;i++) {
303  ix*=this->size[i];
304  ix+=index[i];
305  }
306 
307  // Set value
308  this->data[ix]=val;
309 
310  return;
311  }
312 
313  /** \brief Set the element closest to grid point \c grdp to value
314  \c val
315 
316  The parameters \c closest and \c grdp may be identical,
317  allowing one to update a vector \c grdp with the
318  closest grid point.
319  */
320  template<class vec2_t, class vec3_t>
321  void set_val(const vec2_t &grdp, double val, vec3_t &closest) {
322 
323  // Find indices
324  vec_size_t index(this->rk);
325  for(size_t i=0;i<this->rk;i++) {
326  index[i]=lookup_grid_val(i,grdp[i],closest[i]);
327  }
328 
329  // Pack
330  size_t ix=index[0];
331  for(size_t i=1;i<this->rk;i++) {
332  ix*=this->size[i];
333  ix+=index[i];
334  }
335 
336  // Set value
337  this->data[ix]=val;
338 
339  return;
340  }
341  //@}
342 
343  /// \name Get functions
344  //@{
345  /// Get the element closest to grid point \c gridp
346  template<class vec2_t> double get_val(const vec2_t &gridp) {
347 
348  // Find indices
349  vec_size_t index(this->rk);
350  for(size_t i=0;i<this->rk;i++) {
351  index[i]=lookup_grid(i,gridp[i]);
352  }
353 
354  // Pack
355  size_t ix=index[0];
356  for(size_t i=1;i<this->rk;i++) {
357  ix*=this->size[i];
358  ix+=index[i];
359  }
360 
361  // Set value
362  return this->data[ix];
363  }
364 
365  /** \brief Get the element closest to grid point \c gridp,
366  store grid values in \c closest and return value
367 
368  The parameters \c gridp and \c closest may refer to the
369  same object.
370  */
371  template<class vec2_t, class vec3_t>
372  double get_val(const vec2_t &gridp, vec3_t &closest) {
373 
374  // Find indices
375  vec_size_t index(this->rk);
376  for(size_t i=0;i<this->rk;i++) {
377  index[i]=lookup_grid_val(i,gridp[i],closest[i]);
378  }
379 
380  // Pack
381  size_t ix=index[0];
382  for(size_t i=1;i<this->rk;i++) {
383  ix*=this->size[i];
384  ix+=index[i];
385  }
386 
387  // Set value
388  return this->data[ix];
389  }
390 
391  /// Return a reference to the data (for HDF I/O)
392  vec_t &get_data() {
393  return this->data;
394  }
395  //@}
396 
397  /// \name Resize method
398  //@{
399  /** \brief Resize the tensor to rank \c rank with sizes
400  given in \c dim
401 
402  The parameter \c dim must be a vector of sizes with a length
403  equal to \c rank. This resize method is always destructive,
404  and the grid is always reset.
405 
406  If the user requests any of the sizes to be zero, this
407  function will call the error handler.
408  */
409  template<class size_vec2_t>
410  void resize(size_t rank, const size_vec2_t &dim) {
411  // Double check that none of the sizes that the user
412  // specified are zero
413  for(size_t i=0;i<rank;i++) {
414  if (dim[i]==0) {
415  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
416  "rank for index "+szttos(i)+" in tensor_grid::"+
417  "resize().").c_str(),exc_einval);
418  }
419  }
420  // Set the new rank
421  this->rk=rank;
422  // Resize the size vector
423  this->size.resize(this->rk);
424  // Reset the grid
425  grid_set=false;
426  grid.resize(0);
427  // If the new rank is zero, reset the data, otherwise,
428  // update the size vector and reset the data vector
429  if (rank==0) {
430  this->data.resize(0);
431  return;
432  } else {
433  size_t tot=1;
434  for(size_t i=0;i<this->rk;i++) {
435  this->size[i]=dim[i];
436  tot*=this->size[i];
437  }
438  this->data.resize(tot);
439  }
440  return;
441  }
442 
443  //@}
444 
445  /// \name Grid manipulation
446  //@{
447  /// Return true if the grid has been set
448  bool is_grid_set() const { return grid_set; }
449 
450  /** \brief Set the grid
451 
452  The grid must be specified for all of the dimensions at
453  once. Denote \f$ (\mathrm{size})_0 \f$ as the size
454  of the first dimension, \f$ (\mathrm{size})_1 \f$ as the
455  size of the second dimesion, and so on. Then the
456  first \f$ (\mathrm{size})_0 \f$ entries in \c grid must
457  be the grid for the first dimension, the next
458  \f$ (\mathrm{size})_1 \f$ entries must be the grid
459  for the second dimension, and so on. Thus \c grid
460  must be a vector of size
461  \f[
462  \sum_{i=0}^{\mathrm{rank}} (\mathrm{size})_i
463  \f]
464 
465  Note that the grid is copied so the function argument may
466  be destroyed by the user after calling set_grid_packed() without
467  affecting the tensor grid.
468 
469  \future Define a more generic interface for matrix types
470  */
471  template<class vec2_t>
472  void set_grid_packed(const vec2_t &grid_vec) {
473  if (this->rk==0) {
474  O2SCL_ERR2("Tried to set grid for empty tensor in ",
475  "tensor_grid::set_grid_packed().",exc_einval);
476  }
477  size_t ngrid=0;
478  for(size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
479  grid.resize(ngrid);
480  for(size_t i=0;i<ngrid;i++) {
481  grid[i]=grid_vec[i];
482  }
483  grid_set=true;
484  return;
485  }
486 
487  /** \brief Set grid from a vector of vectors of grid points
488  */
489  template<class vec_vec_t>
490  void set_grid(const vec_vec_t &grid_vecs) {
491  if (this->rk==0) {
492  O2SCL_ERR2("Tried to set grid for empty tensor in ",
493  "tensor_grid::set_grid().",exc_einval);
494  }
495  size_t ngrid=0;
496  for(size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
497  grid.resize(ngrid);
498  size_t k=0;
499  for(size_t i=0;i<this->rk;i++) {
500  for(size_t j=0;j<this->size[i];j++) {
501  grid[k]=grid_vecs[i][j];
502  k++;
503  }
504  }
505  grid_set=true;
506  return;
507  }
508 
509  /** \brief Use a default grid which just uses the index
510  */
511  void default_grid() {
512  size_t ngrid=0;
513  for(size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
514  grid.resize(ngrid);
515  size_t k=0;
516  for(size_t i=0;i<this->rk;i++) {
517  for(size_t j=0;j<this->size[i];j++) {
518  grid[k]=((double)j);
519  k++;
520  }
521  }
522  grid_set=true;
523  return;
524  }
525 
526  /** \brief Set grid for one index from a vector
527  */
528  template<class vec2_t>
529  void set_grid_i_vec(size_t ix, const vec2_t &grid_vec) {
530  if (grid_set==false) {
531  O2SCL_ERR2("Grid not already set in ",
532  "tensor_grid::set_grid_i_vec().",exc_einval);
533  }
534  if (this->rk==0) {
535  O2SCL_ERR2("Tried to set grid for empty tensor in ",
536  "tensor_grid::set_grid_i_vec().",exc_einval);
537  }
538  size_t k=0;
539  for(size_t i=0;i<this->rk;i++) {
540  for(size_t j=0;j<this->size[i];j++) {
541  if (j==ix) {
542  grid[k]=grid_vec[j];
543  }
544  k++;
545  }
546  }
547  return;
548  }
549 
550  /** \brief Set grid for one index from a function
551  */
552  template<class vec2_t>
553  void set_grid_i_func(size_t ix, std::string func) {
554  if (grid_set==false) {
555  O2SCL_ERR2("Grid not already set in ",
556  "tensor_grid::set_grid_i_func().",exc_einval);
557  }
558  if (this->rk==0) {
559  O2SCL_ERR2("Tried to set grid for empty tensor in ",
560  "tensor_grid::set_grid_i_func().",exc_einval);
561  }
562 
563  calculator calc;
564  std::map<std::string,double> vars;
565  calc.compile(func.c_str(),&vars);
566 
567  size_t k=0;
568  for(size_t i=0;i<this->rk;i++) {
569  for(size_t j=0;j<this->size[i];j++) {
570  if (j==ix) {
571  vars["i"]=((double)j);
572  grid[k]=calc.eval(&vars);
573  }
574  k++;
575  }
576  }
577 
578  return;
579  }
580 
581  /** \brief Set grid from a vector of uniform grid objects
582 
583  \note This is called by one of the constructors.
584  */
585  void set_grid(std::vector<uniform_grid<double> > &ugs) {
586  if (this->rk==0) {
587  O2SCL_ERR2("Tried to set grid for empty tensor in ",
588  "tensor_grid::set_grid().",exc_einval);
589  }
590  size_t ngrid=0;
591  for(size_t i=0;i<this->rk;i++) ngrid+=ugs[i].get_npoints();
592  grid.resize(ngrid);
593  size_t k=0;
594  for(size_t i=0;i<this->rk;i++) {
595  for(size_t j=0;j<ugs[i].get_npoints();j++) {
596  grid[k]=ugs[i][j];
597  k++;
598  }
599  }
600  grid_set=true;
601  return;
602  }
603 
604  /** \brief Copy grid for index \c i to vector \c v
605 
606  The type \c rvec_t must be a vector with a resize
607  method.
608  */
609  template<class rvec_t> void copy_grid(size_t i, rvec_t &v) {
610  v.resize(this->size[i]);
611  size_t istart=0;
612  for(size_t k=0;k<i;k++) istart+=this->size[k];
613  for(size_t k=0;k<this->size[i];k++) {
614  v[k]=grid[istart+k];
615  }
616  return;
617  }
618 
619  /// Lookup jth value on the ith grid
620  double get_grid(size_t i, size_t j) const {
621  if (!grid_set) {
622  O2SCL_ERR("Grid not set in tensor_grid::get_grid().",
623  exc_einval);
624  }
625  if (i>=this->rk) {
626  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
627  " greater than or equal to rank, "+szttos(this->rk)+
628  ", in tensor_grid::get_grid().").c_str(),
629  exc_einval);
630  }
631  size_t istart=0;
632  for(size_t k=0;k<i;k++) istart+=this->size[k];
633  return grid[istart+j];
634  }
635 
636  /// Set the jth value on the ith grid
637  void set_grid(size_t i, size_t j, double val) {
638  if (!grid_set) {
639  O2SCL_ERR("Grid not set in tensor_grid::get_grid().",
640  exc_einval);
641  }
642  if (i>=this->rk) {
643  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
644  " greater than or equal to rank, "+szttos(this->rk)+
645  ", in tensor_grid::get_grid().").c_str(),
646  exc_einval);
647  }
648  size_t istart=0;
649  for(size_t k=0;k<i;k++) istart+=this->size[k];
650  grid[istart+j]=val;
651  return;
652  }
653 
654  /** \brief Lookup index for grid closest to \c val, returning the
655  grid point
656 
657  The parameters \c val and \c val2 may refer to the
658  same object.
659  */
660  size_t lookup_grid_val(size_t i, const double &val, double &val2) {
661  if (i>=this->rk) {
662  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
663  " greater than or equal to rank, "+szttos(this->rk)+
664  ", in tensor_grid::lookup_grid_val().").c_str(),
665  exc_einval);
666  }
667  if (grid_set==false) {
668  O2SCL_ERR("Grid not set in tensor_grid::lookup_grid_val().",
669  exc_einval);
670  }
671 
672  // We have to store the grid point in a temporary in case
673  // the parameters gridp and closest refer to the same object
674  double temp=val;
675 
676  size_t istart=0;
677  for(size_t j=0;j<i;j++) istart+=this->size[j];
678  size_t best=istart;
679  double min=fabs(grid[istart]-temp);
680  val2=grid[istart];
681  for(size_t j=istart;j<istart+this->size[i];j++) {
682  if (fabs(grid[j]-temp)<min) {
683  best=j;
684  min=fabs(grid[j]-temp);
685  val2=grid[j];
686  }
687  }
688  return best-istart;
689  }
690 
691  /// Lookup index for grid closest to \c val
692  size_t lookup_grid(size_t i, double val) {
693  double val2;
694  return lookup_grid_val(i,val,val2);
695  }
696 
697  /** \brief Lookup indices for grid closest point to \c vals
698 
699  The values in \c vals are not modified by this function.
700 
701  \comment
702  This function must have a different name than
703  lookup_grid() because the template types cause
704  confusion between the two functions.
705  \endcomment
706  */
707  template<class vec2_t, class size_vec2_t>
708  void lookup_grid_vec(const vec2_t &vals, size_vec2_t &indices) const {
709  for(size_t k=0;k<this->rk;k++) {
710  indices[k]=lookup_grid(k,vals[k]);
711  }
712  return;
713  }
714 
715  /** \brief Lookup internal packed grid index for point closest
716  to \c val and store closest value in \c val2
717 
718  This version, rather than \ref
719  o2scl::tensor_grid::lookup_grid_val() can be useful because it
720  gives the index of the grid point in the internal grid vector
721  object.
722  */
723  size_t lookup_grid_packed_val(size_t i, double val, double &val2) {
724  if (!grid_set) {
725  O2SCL_ERR("Grid not set in tensor_grid::lookup_grid_packed_val().",
726  exc_einval);
727  }
728  if (i>=this->rk) {
729  O2SCL_ERR((((std::string)"Index ")+szttos(i)+" greater than rank, "+
730  szttos(this->rk)+
731  ", in tensor_grid::lookup_grid_packed_val().").c_str(),
732  exc_einval);
733  }
734  size_t istart=0;
735  for(size_t j=0;j<i;j++) istart+=this->size[j];
736  size_t best=istart;
737  double min=fabs(grid[istart]-val);
738  val2=grid[istart];
739  for(size_t j=istart;j<istart+this->size[i];j++) {
740  if (fabs(grid[j]-val)<min) {
741  best=j;
742  min=fabs(grid[j]-val);
743  val2=grid[j];
744  }
745  }
746  return best;
747  }
748 
749  /** \brief Lookup internal packed grid index for point closest
750  to \c val
751  */
752  size_t lookup_grid_packed(size_t i, double val) {
753  double val2;
754  return lookup_grid_packed_val(i,val,val2);
755  }
756  //@}
757 
758  /// \name Slicing to tensor_grid objects
759  //@{
760  /** \brief Copy an abitrary slice by fixing 1 or more indices and
761  use interpolation to return a new \ref tensor_grid object
762  */
763  template<class size_vec2_t, class vec2_t>
764  tensor_grid<> copy_slice_interp(size_vec2_t &ifix, vec2_t &vals) {
765 
766  if (this->rk<1+ifix.size()) {
767  O2SCL_ERR2("Fixed too many indices in ",
768  "tensor_grid::copy_slice_interp().",
770  }
771  if (ifix.size()!=vals.size()) {
772  O2SCL_ERR2("Mismatch between indices and values in ",
773  "tensor_grid::copy_slice_interp().",
775  }
776 
777  // Determine new rank
778  size_t rank_new=this->rk-ifix.size();
779 
780  // Determine the new sizes and new grid
781  std::vector<size_t> sz_new;
782  std::vector<std::vector<double> > grid_new;
783  for(size_t i=0;i<this->rk;i++) {
784  bool found=false;
785  for(size_t j=0;j<ifix.size();j++) {
786  if (ifix[j]==i) found=true;
787  }
788  if (found==false) {
789  sz_new.push_back(this->get_size(i));
790  std::vector<double> grid_temp;
791  for(size_t j=0;j<this->get_size(i);j++) {
792  grid_temp.push_back(this->get_grid(i,j));
793  }
794  grid_new.push_back(grid_temp);
795  }
796  }
797 
798  // Create the new tensor_grid object and set the new grid
799  tensor_grid<> tg_new(rank_new,sz_new);
800  tg_new.set_grid(grid_new);
801 
802  // Interpolate the data into the new tensor_grid object
803  std::vector<size_t> ix_new(rank_new);
804  std::vector<double> point_old(this->rk);
805 
806  // Loop over the new tensor_grid object
807  for(size_t i=0;i<tg_new.total_size();i++) {
808 
809  // Find the location in the new tensor_grid object
810  tg_new.unpack_index(i,ix_new);
811 
812  // Find the point in the old tensor object to interpolate
813  for(size_t j=0;j<this->rk;j++) {
814  int ix_found=-1;
815  for(size_t k=0;k<ifix.size();k++) {
816  if (ifix[k]==j) ix_found=k;
817  }
818  if (ix_found==-1) {
819  point_old[j]=this->get_grid(j,ix_new[j]);
820  } else {
821  point_old[j]=vals[ix_found];
822  }
823  }
824 
825  // Set the new point by performing the linear interpolation
826  tg_new.set(ix_new,this->interp_linear(point_old));
827  }
828 
829  return tg_new;
830  }
831  //@}
832 
833  /// \name Slicing and converting to table3d objects
834  //@{
835 
836  /** \brief Convert to a \ref o2scl::table3d object by
837  summing over all but two indices
838  */
839  void convert_table3d_sum
840  (size_t ix_x, size_t ix_y, table3d &tab, std::string x_name="x",
841  std::string y_name="y", std::string slice_name="z") {
842 
843  // Get current table3d grid
844  size_t nx, ny;
845  tab.get_size(nx,ny);
846 
847  if (nx==0 && ny==0) {
848 
849  if (x_name.length()==0) x_name="x";
850  if (y_name.length()==0) y_name="y";
851 
852  // If there's no grid, then create a grid in the table3d
853  // object that is the same as that in the tensor_grid object
854  std::vector<double> grid_x, grid_y;
855  copy_grid(ix_x,grid_x);
856  copy_grid(ix_y,grid_y);
857  tab.set_xy("x",grid_x.size(),grid_x,
858  "y",grid_y.size(),grid_y);
859  // Now that the grid is set, get nx and ny
860  tab.get_size(nx,ny);
861  }
862 
863  // Check that the grids are commensurate
864  if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
865  O2SCL_ERR2("Grids not commensurate in ",
866  "tensor_grid::convert_table3d_sum().",exc_einval);
867  }
868 
869  tab.set_slice_all(slice_name,0.0);
870 
871  std::vector<size_t> ix;
872  for(size_t i=0;i<this->total_size();i++) {
873  this->unpack_index(i,ix);
874  tab.set(ix[ix_x],ix[ix_y],slice_name,
875  tab.get(ix[ix_x],ix[ix_y],slice_name)+
876  this->data[i]);
877  }
878 
879  return;
880  }
881 
882  /** \brief Create a slice in a \ref o2scl::table3d object with an
883  aligned grid
884 
885  This function uses the grid associated with indices \c ix_x
886  and \c ix_y, to copy data to a slice named \c slice_name in
887  the table3d object \c tab . All other indices are fixed
888  to values specified by the user in \c index and the values
889  of <tt>index[ix_x]</tt> and <tt>index[ix_y]</tt> are
890  used for temporary storage.
891 
892  If the table3d object does not currently have a grid set, then
893  the grid is automatically set to be the same as that stored in
894  the tensor_grid object associated with ranks \c ix_x and \c
895  iy_y. If the \ref o2scl::table3d object does have a grid set,
896  then the values returned by \ref o2scl::table3d::get_nx() and
897  \ref o2scl::table3d::get_ny() must be equal to the size of the
898  tensor in indices \c ix_x and ix_y, respectively. If a
899  slice named \c slice_name is not already present in
900  \c tab, then a new slice with that name is created.
901 
902  The error handler is called if \c ix_x is the same as
903  \c ix_y, or if either of these two values is greater
904  than or equal to the tensor rank.
905  */
906  template<class size_vec2_t>
907  void copy_table3d_align(size_t ix_x, size_t ix_y, size_vec2_t &index,
908  table3d &tab, std::string slice_name="z") {
909 
910  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
911  O2SCL_ERR2("Either indices greater than rank or x and y ind",
912  "ices equal in tensor_grid::copy_table3d_align().",
913  exc_efailed);
914  }
915 
916  // Get current table3d grid
917  size_t nx, ny;
918  tab.get_size(nx,ny);
919 
920  // Check that the grids are commensurate
921  if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
922  O2SCL_ERR2("Grids not commensurate in ",
923  "tensor_grid::copy_table3d_align().",exc_einval);
924  }
925 
926  // Create slice if not already present
927  size_t is;
928  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
929 
930  // Copy over data
931  for(size_t i=0;i<nx;i++) {
932  for(size_t j=0;j<ny;j++) {
933  index[ix_x]=i;
934  index[ix_y]=j;
935  double val=this->get(index);
936  tab.set(i,j,slice_name,val);
937  }
938  }
939 
940  return;
941  }
942 
943  /** \brief Create a slice in a table3d object with a new aligned
944  grid
945  */
946  template<class size_vec2_t>
947  void copy_table3d_align_setxy
948  (size_t ix_x, size_t ix_y, size_vec2_t &index,
949  table3d &tab, std::string x_name="x", std::string y_name="y",
950  std::string slice_name="z") {
951 
952  // Get current table3d grid
953  size_t nx, ny;
954  tab.get_size(nx,ny);
955 
956  if (nx==0 && ny==0) {
957 
958  if (x_name.length()==0) x_name="x";
959  if (y_name.length()==0) y_name="y";
960 
961  // If there's no grid, then create a grid in the table3d
962  // object that is the same as that in the tensor_grid object
963  std::vector<double> grid_x, grid_y;
964  copy_grid(ix_x,grid_x);
965  copy_grid(ix_y,grid_y);
966  tab.set_xy("x",grid_x.size(),grid_x,
967  "y",grid_y.size(),grid_y);
968  // Now that the grid is set, get nx and ny
969  tab.get_size(nx,ny);
970  }
971 
972  copy_table3d_align(ix_x,ix_y,index,tab,slice_name);
973 
974  return;
975  }
976 
977  /** \brief Copy to a slice in a table3d object using interpolation
978 
979  This function uses the grid associated with indices \c ix_x
980  and \c ix_y, and the tensor interpolation function to copy the
981  tensor information to the slice named \c slice_name in the
982  table3d object \c tab . All other indices are fixed
983  to values specified by the user in \c index and the values
984  of <tt>index[ix_x]</tt> and <tt>index[ix_y]</tt> are
985  used for temporary storage.
986 
987  If a slice named \c slice_name is not already present in \c
988  tab, then a new slice with that name is created.
989 
990  The error handler is called if \c ix_x is the same as
991  \c ix_y, or if either of these two values is greater
992  than or equal to the tensor rank.
993 
994  \note This function uses the \ref tensor_grid::interp_linear()
995  for the interpolation.
996  */
997  template<class size_vec2_t>
998  void copy_table3d_interp(size_t ix_x, size_t ix_y, size_vec2_t &index,
999  table3d &tab, std::string slice_name="z") {
1000 
1001  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
1002  O2SCL_ERR2("Either indices greater than rank or x and y ",
1003  "indices equal in tensor_grid::copy_table3d_interp().",
1004  exc_efailed);
1005  }
1006 
1007  // Get current table3d grid
1008  size_t nx, ny;
1009  tab.get_size(nx,ny);
1010 
1011  if (nx==0 && ny==0) {
1012  // If there's no grid, then just use the aligned version
1013  return copy_table3d_align(ix_x,ix_y,index,tab,slice_name);
1014  }
1015 
1016  // Create vector of values to interpolate with
1017  std::vector<double> vals(this->rk);
1018  for(size_t i=0;i<this->rk;i++) {
1019  if (i!=ix_x && i!=ix_y) vals[i]=this->get_grid(i,index[i]);
1020  }
1021 
1022  // Create slice if not already present
1023  size_t is;
1024  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
1025 
1026  // Loop through the table grid to perform the interpolation
1027  for(size_t i=0;i<nx;i++) {
1028  for(size_t j=0;j<ny;j++) {
1029  vals[ix_x]=tab.get_grid_x(i);
1030  vals[ix_y]=tab.get_grid_y(j);
1031  tab.set(i,j,slice_name,this->interp_linear(vals));
1032  }
1033  }
1034 
1035  return;
1036  }
1037 
1038  /** \brief Copy to a slice in a table3d object using interpolation
1039  */
1040  template<class vec2_t>
1041  void copy_table3d_interp_values(size_t ix_x, size_t ix_y,
1042  vec2_t &values, table3d &tab,
1043  std::string slice_name="z",
1044  int verbose=0) {
1045 
1046  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
1047  O2SCL_ERR2("Either indices greater than rank or x and y ",
1048  "indices equal in tensor_grid::copy_table3d_interp().",
1049  exc_efailed);
1050  }
1051  if (values.size()!=this->rk) {
1052  O2SCL_ERR2("Values array not equal to rank ",
1053  "in tensor_grid::copy_table3d_interp_values().",
1054  exc_efailed);
1055  }
1056 
1057  if (tab.is_size_set()==false || tab.is_xy_set()==false) {
1058  O2SCL_ERR2("Grid not set in tensor_grid::",
1059  "copy_table3d_interp_value().",o2scl::exc_einval);
1060  }
1061 
1062  // Get current table3d grid
1063  size_t nx, ny;
1064  tab.get_size(nx,ny);
1065 
1066  // Create slice if not already present
1067  size_t is;
1068  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
1069 
1070  // Loop through the table grid to perform the interpolation
1071  for(size_t i=0;i<nx;i++) {
1072  for(size_t j=0;j<ny;j++) {
1073  values[ix_x]=tab.get_grid_x(i);
1074  values[ix_y]=tab.get_grid_y(j);
1075  tab.set(i,j,slice_name,this->interp_linear(values));
1076  if (verbose>0) {
1077  std::cout << "At location values: ";
1078  for(size_t k=0;k<values.size();k++) {
1079  std::cout << values[k] << " ";
1080  }
1081  std::cout << "Interpolated to get: "
1082  << i << " " << j << " " << slice_name << " "
1083  << this->interp_linear(values) << std::endl;
1084  if (verbose>1) {
1085  char ch;
1086  std::cin >> ch;
1087  }
1088  }
1089  }
1090  }
1091 
1092  return;
1093  }
1094 
1095  /** \brief Copy to a slice in a table3d object using interpolation
1096  creating a new table3d grid
1097  */
1098  template<class vec2_t>
1099  void copy_table3d_interp_values_setxy
1100  (size_t ix_x, size_t ix_y, vec2_t &values, table3d &tab,
1101  std::string x_name="x", std::string y_name="y",
1102  std::string slice_name="z") {
1103 
1104  // Get current table3d grid
1105  size_t nx, ny;
1106  tab.get_size(nx,ny);
1107 
1108  if (nx==0 && ny==0) {
1109 
1110  if (x_name.length()==0) x_name="x";
1111  if (y_name.length()==0) y_name="y";
1112 
1113  // If there's no grid, then create a grid in the table3d
1114  // object that is the same as that in the tensor_grid object
1115  std::vector<double> grid_x, grid_y;
1116  copy_grid(ix_x,grid_x);
1117  copy_grid(ix_y,grid_y);
1118  tab.set_xy("x",grid_x.size(),grid_x,
1119  "y",grid_y.size(),grid_y);
1120  // Now that the grid is set, get nx and ny
1121  tab.get_size(nx,ny);
1122  }
1123 
1124  copy_table3d_interp_values(ix_x,ix_y,values,tab,slice_name);
1125  return;
1126  }
1127  //@}
1128 
1129  /// \name Clear method
1130  //@{
1131  /// Clear the tensor of all data and free allocated memory
1132  void clear() {
1133  grid.resize(0);
1134  grid_set=false;
1136  return;
1137  }
1138  //@}
1139 
1140  /// \name Interpolation
1141  //@{
1142  /// Set interpolation type for \ref interpolate()
1143  void set_interp_type(size_t interp_type) {
1144  itype=interp_type;
1145  return;
1146  }
1147 
1148  /** \brief Interpolate values \c vals into the tensor,
1149  returning the result
1150 
1151  This is a quick and dirty implementation of n-dimensional
1152  interpolation by recursive application of the 1-dimensional
1153  routine from \ref interp_vec, using the base interpolation
1154  object specified in the template parameter \c base_interp_t.
1155  This will be very slow for sufficiently large data sets.
1156 
1157  \note This function requires a range objects to obtain ranges
1158  of vector objects. In ublas, this is done with
1159  <tt>ublas::vector_range</tt> objects, so this function will
1160  certainly work for \ref tensor_grid objects built on ublas
1161  vector types. There is no corresponding <tt>std::range</tt>,
1162  but you may be able to use either <tt>ublas::vector_range</tt>
1163  or <tt>Boost.Range</tt> in order to call this function with
1164  \ref tensor_grid objects built on <tt>std::vector</tt>.
1165  However, this is not fully tested at the moment.
1166 
1167  \future It should be straightforward to improve the scaling of
1168  this algorithm significantly by creating a "window" of local
1169  points around the point of interest. This could be done easily
1170  by constructing an initial subtensor. However, this should
1171  probably be superceded by a more generic alternative which
1172  avoids explicit use of the 1-d interpolation types.
1173  */
1174  template<class range_t=ub_range,
1175  class data_range_t=ubvector_range,
1176  class index_range_t=ubvector_size_t_range>
1177  double interpolate(double *vals) {
1178 
1179  typedef interp_vec<vec_t> interp_t;
1180 
1181  if (this->rk==1) {
1182 
1183  interp_t si(this->size[0],grid,this->data,itype);
1184  return si.eval(vals[0]);
1185 
1186  } else {
1187 
1188  // Get total number of interpolations at this level
1189  size_t ss=1;
1190  for(size_t i=1;i<this->rk;i++) ss*=this->size[i];
1191 
1192  // Create space for y vectors and interpolators
1193  std::vector<vec_t> yvec(ss);
1194  std::vector<interp_t *> si(ss);
1195  for(size_t i=0;i<ss;i++) {
1196  yvec[i].resize(this->size[0]);
1197  }
1198 
1199  // Create space for interpolation results
1200  tensor_grid tdat;
1201  index_range_t size_new(this->size,ub_range(1,this->rk));
1202  tdat.resize(this->rk-1,size_new);
1203 
1204  // Set grid for temporary tensor
1205  data_range_t grid_new(grid,ub_range(this->size[0],grid.size()));
1206  tdat.set_grid_packed(grid_new);
1207 
1208  // Create starting coordinate and counter
1209  vec_size_t co(this->rk);
1210  for(size_t i=0;i<this->rk;i++) co[i]=0;
1211  size_t cnt=0;
1212 
1213  // Loop over every interpolation
1214  bool done=false;
1215  while(done==false) {
1216 
1217  // Fill yvector with the appropriate data
1218  for(size_t i=0;i<this->size[0];i++) {
1219  co[0]=i;
1220  yvec[cnt][i]=this->get(co);
1221  }
1222 
1223  si[cnt]=new interp_t(this->size[0],grid,yvec[cnt],itype);
1224 
1225  index_range_t co2(co,ub_range(1,this->rk));
1226  tdat.set(co2,si[cnt]->eval(vals[0]));
1227 
1228  // Go to next interpolation
1229  cnt++;
1230  co[this->rk-1]++;
1231  // carry if necessary
1232  for(int j=((int)this->rk)-1;j>0;j--) {
1233  if (co[j]>=this->size[j]) {
1234  co[j]=0;
1235  co[j-1]++;
1236  }
1237  }
1238 
1239  // Test if done
1240  if (cnt==ss) done=true;
1241 
1242  // End of while loop
1243  }
1244 
1245  // Now call the next level of interpolation
1246  double res=tdat.interpolate(vals+1);
1247 
1248  for(size_t i=0;i<ss;i++) {
1249  delete si[i];
1250  }
1251 
1252  return res;
1253  }
1254  }
1255 
1256  /** \brief Perform a linear interpolation of \c v into the
1257  function implied by the tensor and grid
1258 
1259  This performs multi-dimensional linear interpolation (or
1260  extrapolation) It works by first using \ref o2scl::search_vec
1261  to find the interval containing (or closest to) the specified
1262  point in each direction and constructing the corresponding
1263  hypercube of size \f$ 2^{\mathrm{rank}} \f$ containing \c v.
1264  It then calls \ref interp_linear_power_two() to perform the
1265  interpolation in that hypercube.
1266 
1267  \future This starts with a small copy, which can be eliminated
1268  by creating a new version of interp_linear_power_two
1269  which accepts an offset vector parameter so that the
1270  first interpolation is offset. Remaining interpolations
1271  don't need to be offset because the tensor has to be
1272  created from the previous interpolation round.
1273  */
1274  template<class vec2_size_t> double interp_linear(vec2_size_t &v) {
1275 
1276  // Find the the corner of the hypercube containing v
1277  size_t rgs=0;
1278  std::vector<size_t> loc(this->rk);
1279  std::vector<double> gnew;
1280  for(size_t i=0;i<this->rk;i++) {
1281  std::vector<double> grid_unpacked(this->size[i]);
1282  for(size_t j=0;j<this->size[i];j++) {
1283  grid_unpacked[j]=grid[j+rgs];
1284  }
1285  search_vec<std::vector<double> > sv(this->size[i],grid_unpacked);
1286  loc[i]=sv.find(v[i]);
1287  gnew.push_back(grid_unpacked[loc[i]]);
1288  gnew.push_back(grid_unpacked[loc[i]+1]);
1289  rgs+=this->size[i];
1290  }
1291 
1292  // Now construct a 2^{rk}-sized tensor containing only that
1293  // hypercube
1294  std::vector<size_t> snew(this->rk);
1295  for(size_t i=0;i<this->rk;i++) {
1296  snew[i]=2;
1297  }
1298  tensor_grid tnew(this->rk,snew);
1299  tnew.set_grid_packed(gnew);
1300 
1301  // Copy over the relevant data
1302  for(size_t i=0;i<tnew.total_size();i++) {
1303  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1304  tnew.unpack_index(i,index_new);
1305  for(size_t j=0;j<this->rk;j++) index_old[j]=index_new[j]+loc[j];
1306  tnew.set(index_new,this->get(index_old));
1307  }
1308 
1309  // Now use interp_power_two()
1310  return tnew.interp_linear_power_two(v);
1311  }
1312 
1313  /** \brief Perform linear interpolation assuming that all
1314  indices can take only two values
1315 
1316  This function works by recursively slicing the hypercube of
1317  size \f$ 2^{\mathrm{rank}} \f$ into a hypercube of size \f$
1318  2^{\mathrm{rank-1}} \f$ performing linear interpolation for
1319  each pair of points.
1320 
1321  \note This is principally a function for internal use
1322  by \ref interp_linear().
1323  */
1324  template<class vec2_size_t>
1325  double interp_linear_power_two(vec2_size_t &v) {
1326 
1327  if (this->rk==1) {
1328  return this->data[0]+(this->data[1]-this->data[0])/
1329  (grid[1]-grid[0])*(v[0]-grid[0]);
1330  }
1331 
1332  size_t last=this->rk-1;
1333  double frac=(v[last]-get_grid(last,0))/
1334  (get_grid(last,1)-get_grid(last,0));
1335 
1336  // Create new size vector and grid
1337  tensor_grid tnew(this->rk-1,this->size);
1338  tnew.set_grid_packed(grid);
1339 
1340  // Create data in new tensor, removing the last index through
1341  // linear interpolation
1342  for(size_t i=0;i<tnew.total_size();i++) {
1343  std::vector<size_t> index(this->rk);
1344  tnew.unpack_index(i,index);
1345  index[this->rk-1]=0;
1346  double val_lo=this->get(index);
1347  index[this->rk-1]=1;
1348  double val_hi=this->get(index);
1349  tnew.set(index,val_lo+frac*(val_hi-val_lo));
1350  }
1351 
1352  // Recursively interpolate the smaller tensor
1353  return tnew.interp_linear_power_two(v);
1354  }
1355 
1356  /** \brief Perform a linear interpolation of <tt>v[1]</tt>
1357  to <tt>v[n-1]</tt> resulting in a vector
1358 
1359  \note The type <tt>vec2_t</tt> for the vector <tt>res</tt>
1360  must have a <tt>resize()</tt> method.
1361 
1362  This performs multi-dimensional linear interpolation (or
1363  extrapolation) in the last <tt>n-1</tt> indices of the
1364  rank-<tt>n</tt> tensor leaving the first index free and places
1365  the results in the vector \c res.
1366  */
1367  template<class vec2_size_t, class vec2_t>
1368  void interp_linear_vec0(vec2_size_t &v, vec2_t &res) {
1369 
1370  // Find the the corner of the hypercube containing v
1371  size_t rgs=0;
1372  std::vector<size_t> loc(this->rk);
1373  std::vector<double> gnew;
1374  for(size_t i=0;i<this->size[0];i++) {
1375  gnew.push_back(grid[i]);
1376  }
1377  rgs=this->size[0];
1378  loc[0]=0;
1379  for(size_t i=1;i<this->rk;i++) {
1380  std::vector<double> grid_unpacked(this->size[i]);
1381  for(size_t j=0;j<this->size[i];j++) {
1382  grid_unpacked[j]=grid[j+rgs];
1383  }
1384  search_vec<std::vector<double> > sv(this->size[i],grid_unpacked);
1385  loc[i]=sv.find(v[i]);
1386  gnew.push_back(grid_unpacked[loc[i]]);
1387  gnew.push_back(grid_unpacked[loc[i]+1]);
1388  rgs+=this->size[i];
1389  }
1390 
1391  // Now construct a n*2^{rk-1}-sized tensor containing only that
1392  // hypercube
1393  std::vector<size_t> snew(this->rk);
1394  snew[0]=this->size[0];
1395  for(size_t i=1;i<this->rk;i++) {
1396  snew[i]=2;
1397  }
1398  tensor_grid tnew(this->rk,snew);
1399  tnew.set_grid_packed(gnew);
1400 
1401  // Copy over the relevant data
1402  for(size_t i=0;i<tnew.total_size();i++) {
1403  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1404  tnew.unpack_index(i,index_new);
1405  for(size_t j=0;j<this->rk;j++) {
1406  index_old[j]=index_new[j]+loc[j];
1407  }
1408  tnew.set(index_new,this->get(index_old));
1409  }
1410 
1411  // Now use interp_power_two_vec0()
1412  tnew.interp_linear_power_two_vec0(v,res);
1413 
1414  return;
1415  }
1416 
1417  /** \brief Perform linear interpolation assuming that the last
1418  <tt>n-1</tt> indices can take only two values
1419 
1420  This function performs linear interpolation assuming that the
1421  last <tt>n-1</tt> indices can take only two values and placing
1422  the result into <tt>res</tt>.
1423 
1424  \note The type <tt>vec2_t</tt> for the vector <tt>res</tt>
1425  must have a <tt>resize()</tt> method. This is principally a
1426  function for internal use by \ref interp_linear_vec0().
1427  */
1428  template<class vec2_size_t, class vec2_t>
1429  void interp_linear_power_two_vec0(vec2_size_t &v, vec2_t &res) {
1430 
1431  if (this->rk==2) {
1432  size_t n=this->size[0];
1433  res.resize(n);
1434  vec_size_t ix0(2), ix1(2);
1435  ix0[1]=0;
1436  ix1[1]=1;
1437  for(size_t i=0;i<n;i++) {
1438  ix0[0]=i;
1439  ix1[0]=i;
1440  res[i]=this->get(ix0)+(this->get(ix1)-this->get(ix0))/
1441  (grid[n+1]-grid[n])*(v[1]-grid[n]);
1442  }
1443  return;
1444  }
1445 
1446  size_t last=this->rk-1;
1447  double frac=(v[last]-get_grid(last,0))/
1448  (get_grid(last,1)-get_grid(last,0));
1449 
1450  // Create new size vector and grid
1451  tensor_grid tnew(this->rk-1,this->size);
1452  tnew.set_grid_packed(grid);
1453 
1454  // Create data in new tensor, removing the last index through
1455  // linear interpolation
1456  for(size_t i=0;i<tnew.total_size();i++) {
1457  std::vector<size_t> index(this->rk);
1458  tnew.unpack_index(i,index);
1459  index[this->rk-1]=0;
1460  double val_lo=this->get(index);
1461  index[this->rk-1]=1;
1462  double val_hi=this->get(index);
1463  tnew.set(index,val_lo+frac*(val_hi-val_lo));
1464  }
1465 
1466  // Recursively interpolate the smaller tensor
1467  tnew.interp_linear_power_two_vec0(v,res);
1468 
1469  return;
1470  }
1471 
1472  /** \brief Perform a linear interpolation of <tt>v</tt>
1473  into the tensor leaving one index free resulting in a vector
1474 
1475  This performs multi-dimensional linear interpolation (or
1476  extrapolation) in the last <tt>n-1</tt> indices of the
1477  rank-<tt>n</tt> tensor leaving the first index free and places
1478  the results in the vector \c res.
1479 
1480  \future This function could be more efficient.
1481  */
1482  template<class vec2_size_t, class vec2_t>
1483  void interp_linear_vec(vec2_size_t &v, size_t ifree, vec2_t &res) {
1484 
1485  size_t n=this->size[ifree];
1486 
1487  // This function uses interp_linear_power_two_vec0(), so it
1488  // works by remapping the indices. This defines the remapping.
1489  std::vector<size_t> map;
1490  map.push_back(ifree);
1491  for(size_t i=0;i<this->rk;i++) {
1492  if (i!=ifree) {
1493  map.push_back(i);
1494  }
1495  }
1496 
1497  // Find the the corner of the hypercube containing v
1498  size_t rgs=0;
1499  vec_size_t loc(this->rk);
1500  loc[ifree]=0;
1501  for(size_t i=0;i<this->rk;i++) {
1502  vec_t grid_unpacked(this->size[i]);
1503  for(size_t j=0;j<this->size[i];j++) {
1504  grid_unpacked[j]=grid[j+rgs];
1505  }
1506  search_vec<vec_t> sv(this->size[i],grid_unpacked);
1507  if (i!=ifree) {
1508  loc[i]=sv.find(v[i]);
1509  }
1510  rgs+=this->size[i];
1511  }
1512 
1513  // Compute the remapped grid and interpolating vector
1514  std::vector<double> gnew, vnew;
1515  for(size_t new_ix=0;new_ix<this->rk;new_ix++) {
1516  for(size_t old_ix=0;old_ix<this->rk;old_ix++) {
1517  if (map[new_ix]==old_ix) {
1518  vnew.push_back(v[old_ix]);
1519  if (old_ix==ifree) {
1520  for(size_t j=0;j<this->size[old_ix];j++) {
1521  gnew.push_back(this->get_grid(old_ix,j));
1522  }
1523  } else {
1524  gnew.push_back(this->get_grid(old_ix,loc[old_ix]));
1525  gnew.push_back(this->get_grid(old_ix,loc[old_ix]+1));
1526  }
1527  }
1528  }
1529  }
1530 
1531  // Now construct a n*2^{rk-1}-sized tensor containing only the
1532  // hypercube needed to do the interpolation
1533 
1534  // Specify the size of each rank
1535  std::vector<size_t> snew;
1536  snew.push_back(n);
1537  for(size_t i=0;i<this->rk;i++) {
1538  if (i!=ifree) {
1539  snew.push_back(2);
1540  }
1541  }
1542 
1543  // Create the tensor and set the grid
1544  tensor_grid tnew(this->rk,snew);
1545  tnew.set_grid_packed(gnew);
1546 
1547  // Copy over the relevant data
1548  for(size_t i=0;i<tnew.total_size();i++) {
1549  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1550  tnew.unpack_index(i,index_new);
1551  for(size_t j=0;j<this->rk;j++) {
1552  index_old[map[j]]=index_new[j]+loc[map[j]];
1553  }
1554  tnew.set(index_new,this->get(index_old));
1555  }
1556 
1557  // Now use interp_power_two_vec()
1558  tnew.interp_linear_power_two_vec0(vnew,res);
1559 
1560  return;
1561  }
1562  //@}
1563 
1564  template<class vecf_t, class vecf_size_t> friend void o2scl_hdf::hdf_output
1566  std::string name);
1567 
1568  template<class vecf_t, class vecf_size_t> friend void o2scl_hdf::hdf_input
1570  std::string name);
1571 
1572  };
1573 
1574  /** \brief Rank 1 tensor with a grid
1575 
1576  \future Make rank-specific get_val and set_val functions?
1577  */
1578  template<class vec_t=std::vector<double>,
1579  class vec_size_t=std::vector<size_t> > class tensor_grid1 :
1580  public tensor_grid<vec_t,vec_size_t> {
1581 
1582  public:
1583 
1584  /// Create an empty tensor
1585  tensor_grid1() : tensor_grid<vec_t,vec_size_t>() {}
1586 
1587  /// Create a rank 2 tensor of size \c (sz,sz2,sz3)
1588  tensor_grid1(size_t sz) : tensor_grid<vec_t,vec_size_t>() {
1589  this->rk=1;
1590  this->size.resize(1);
1591  this->size[0]=sz;
1592  this->data.resize(sz);
1593  this->grid_set=false;
1594  }
1595 
1596 #ifdef O2SCL_NEVER_DEFINED
1597  }{
1598 #endif
1599 
1600  virtual ~tensor_grid1() {
1601  }
1602 
1603  /// Get the element indexed by \c (ix1)
1604  double &get(size_t ix1) {
1605  size_t sz[1]={ix1};
1607  }
1608 
1609  /// Get the element indexed by \c (ix1)
1610  const double &get(size_t ix1) const {
1611  size_t sz[1]={ix1};
1613  }
1614 
1615  /// Set the element indexed by \c (ix1) to value \c val
1616  void set(size_t ix1, double val) {
1617  size_t sz[1]={ix1};
1619  }
1620 
1621  /// Interpolate \c x and return the results
1622  template<class range_t=ub_range, class data_range_t=ubvector_range,
1623  class index_range_t=ubvector_size_t_range>
1624  double interp(double x) {
1625  return tensor_grid<vec_t,vec_size_t>::template interpolate
1626  <range_t,data_range_t,index_range_t>(&x);
1627  }
1628 
1629  /// Interpolate \c x and return the results
1630  double interp_linear(double x) {
1631  double arr[1]={x};
1633  }
1634 
1635  };
1636 
1637  /** \brief Rank 2 tensor with a grid
1638  */
1639  template<class vec_t=std::vector<double>,
1640  class vec_size_t=std::vector<size_t> > class tensor_grid2 :
1641  public tensor_grid<vec_t,vec_size_t> {
1642 
1643  public:
1644 
1645  /// Create an empty tensor
1646  tensor_grid2() : tensor_grid<vec_t,vec_size_t>() {}
1647 
1648  /// Create a rank 2 tensor of size \c (sz,sz2)
1649  tensor_grid2(size_t sz, size_t sz2) : tensor_grid<vec_t,vec_size_t>() {
1650  this->rk=2;
1651  this->size.resize(2);
1652  this->size[0]=sz;
1653  this->size[1]=sz2;
1654  size_t tot=sz*sz2;
1655  this->data.resize(tot);
1656  this->grid_set=false;
1657  }
1658 
1659 #ifdef O2SCL_NEVER_DEFINED
1660  }{
1661 #endif
1662 
1663  virtual ~tensor_grid2() {
1664  }
1665 
1666  /// Get the element indexed by \c (ix1,ix2)
1667  double &get(size_t ix1, size_t ix2) {
1668  size_t sz[2]={ix1,ix2};
1670  }
1671 
1672  /// Get the element indexed by \c (ix1,ix2)
1673  const double &get(size_t ix1, size_t ix2) const {
1674  size_t sz[2]={ix1,ix2};
1676  }
1677 
1678  /// Set the element indexed by \c (ix1,ix2) to value \c val
1679  void set(size_t ix1, size_t ix2, double val) {
1680  size_t sz[2]={ix1,ix2};
1682  return;
1683  }
1684 
1685  /// Interpolate \c (x,y) and return the results
1686  template<class range_t=ub_range, class data_range_t=ubvector_range,
1687  class index_range_t=ubvector_size_t_range>
1688  double interp(double x, double y) {
1689  double arr[2]={x,y};
1690  return tensor_grid<vec_t,vec_size_t>::template interpolate
1691  <range_t,data_range_t,index_range_t>(arr);
1692  }
1693 
1694  /// Interpolate \c (x,y) and return the results
1695  double interp_linear(double x, double y) {
1696  double arr[2]={x,y};
1698  }
1699 
1700  };
1701 
1702  /** \brief Rank 3 tensor with a grid
1703  */
1704  template<class vec_t=std::vector<double>,
1705  class vec_size_t=std::vector<size_t> > class tensor_grid3 :
1706  public tensor_grid<vec_t,vec_size_t> {
1707 
1708  public:
1709 
1710  /// Create an empty tensor
1711  tensor_grid3() : tensor_grid<vec_t,vec_size_t>() {}
1712 
1713  /// Create a rank 3 tensor of size \c (sz,sz2,sz3)
1714  tensor_grid3(size_t sz, size_t sz2, size_t sz3) :
1715  tensor_grid<vec_t,vec_size_t>() {
1716  this->rk=3;
1717  this->size.resize(3);
1718  this->size[0]=sz;
1719  this->size[1]=sz2;
1720  this->size[2]=sz3;
1721  size_t tot=sz*sz2*sz3;
1722  this->data.resize(tot);
1723  this->grid_set=false;
1724  }
1725 
1726 #ifdef O2SCL_NEVER_DEFINED
1727  }{
1728 #endif
1729 
1730  virtual ~tensor_grid3() {
1731  }
1732 
1733  /// Get the element indexed by \c (ix1,ix2,ix3)
1734  double &get(size_t ix1, size_t ix2, size_t ix3) {
1735  size_t sz[3]={ix1,ix2,ix3};
1737  }
1738 
1739  /// Get the element indexed by \c (ix1,ix2,ix3)
1740  const double &get(size_t ix1, size_t ix2, size_t ix3) const {
1741  size_t sz[3]={ix1,ix2,ix3};
1743  }
1744 
1745  /// Set the element indexed by \c (ix1,ix2,ix3) to value \c val
1746  void set(size_t ix1, size_t ix2, size_t ix3, double val) {
1747  size_t sz[3]={ix1,ix2, ix3};
1749  return;
1750  }
1751 
1752  /// Interpolate \c (x,y,z) and return the results
1753  template<class range_t=ub_range, class data_range_t=ubvector_range,
1754  class index_range_t=ubvector_size_t_range>
1755  double interp(double x, double y, double z) {
1756  double arr[3]={x,y,z};
1757  return tensor_grid<vec_t,vec_size_t>::template interpolate
1758  <range_t,data_range_t,index_range_t>(arr);
1759  }
1760 
1761  /// Interpolate \c (x,y,z) and return the results
1762  double interp_linear(double x, double y, double z) {
1763  double arr[3]={x,y,z};
1765  }
1766 
1767  };
1768 
1769  /** \brief Rank 4 tensor with a grid
1770  */
1771  template<class vec_t=std::vector<double>,
1772  class vec_size_t=std::vector<size_t> > class tensor_grid4 :
1773  public tensor_grid<vec_t,vec_size_t> {
1774 
1775  public:
1776 
1777  /// Create an empty tensor
1778  tensor_grid4() : tensor_grid<vec_t,vec_size_t>() {}
1779 
1780  /// Create a rank 4 tensor of size \c (sz,sz2,sz3,sz4)
1781  tensor_grid4(size_t sz, size_t sz2, size_t sz3,
1782  size_t sz4) : tensor_grid<vec_t,vec_size_t>() {
1783  this->rk=4;
1784  this->size.resize(4);
1785  this->size[0]=sz;
1786  this->size[1]=sz2;
1787  this->size[2]=sz3;
1788  this->size[3]=sz4;
1789  size_t tot=sz*sz2*sz3*sz4;
1790  this->data.resize(tot);
1791  this->grid_set=false;
1792  }
1793 
1794 #ifdef O2SCL_NEVER_DEFINED
1795  }{
1796 #endif
1797 
1798  virtual ~tensor_grid4() {
1799  }
1800 
1801  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1802  double &get(size_t ix1, size_t ix2, size_t ix3, size_t ix4) {
1803  size_t sz[4]={ix1,ix2,ix3,ix4};
1805  }
1806 
1807  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1808  const double &get(size_t ix1, size_t ix2, size_t ix3,
1809  size_t ix4) const {
1810  size_t sz[4]={ix1,ix2,ix3,ix4};
1812  }
1813 
1814  /// Set the element indexed by \c (ix1,ix2,ix3,ix4) to value \c val
1815  void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4,
1816  double val) {
1817  size_t sz[4]={ix1,ix2,ix3,ix4};
1819  return;
1820  }
1821 
1822  /// Interpolate \c (x,y,z,a) and return the results
1823  template<class range_t=ub_range, class data_range_t=ubvector_range,
1824  class index_range_t=ubvector_size_t_range>
1825  double interp(double x, double y, double z, double a) {
1826  double arr[4]={x,y,z,a};
1827  return tensor_grid<vec_t,vec_size_t>::template interpolate
1828  <range_t,data_range_t,index_range_t>(arr);
1829  }
1830 
1831  /// Interpolate \c (x,y,z,a) and return the results
1832  double interp_linear(double x, double y, double z, double a) {
1833  double arr[4]={x,y,z,a};
1835  }
1836 
1837  };
1838 
1839 #ifndef DOXYGEN_NO_O2NS
1840 }
1841 #endif
1842 
1843 #endif
1844 
1845 
1846 
bool is_grid_set() const
Return true if the grid has been set.
Definition: tensor_grid.h:448
double get_grid_x(size_t ix) const
Get x grid point at index ix.
void unpack_index(size_t ix, size_vec_t &index)
Unpack the single vector index into indices.
Definition: tensor.h:527
Tensor class with arbitrary dimensions with a grid.
Definition: tensor_grid.h:46
tensor_grid()
Create an empty tensor with zero rank.
Definition: tensor_grid.h:178
void interp_linear_power_two_vec0(vec2_size_t &v, vec2_t &res)
Perform linear interpolation assuming that the last n-1 indices can take only two values...
Definition: tensor_grid.h:1429
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
Rank 1 tensor with a grid.
Definition: tensor_grid.h:1579
size_t lookup_grid_packed(size_t i, double val)
Lookup internal packed grid index for point closest to val.
Definition: tensor_grid.h:752
double interpolate(double *vals)
Interpolate values vals into the tensor, returning the result.
Definition: tensor_grid.h:1177
virtual ~tensor_grid()
Destructor.
Definition: tensor_grid.h:223
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.
size_t lookup_grid_val(size_t i, const double &val, double &val2)
Lookup index for grid closest to val, returning the grid point.
Definition: tensor_grid.h:660
size_t lookup_grid_packed_val(size_t i, double val, double &val2)
Lookup internal packed grid index for point closest to val and store closest value in val2...
Definition: tensor_grid.h:723
size_t find(const double x0)
Search an increasing or decreasing vector for the interval containing x0
Definition: search_vec.h:123
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:190
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
double & get(size_t ix, size_t iy, std::string name)
Get element in slice name at location ix,iy
invalid argument supplied by user
Definition: err_hnd.h:59
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor.h:234
void set_grid(std::vector< uniform_grid< double > > &ugs)
Set grid from a vector of uniform grid objects.
Definition: tensor_grid.h:585
Rank 2 tensor with a grid.
Definition: tensor_grid.h:1640
double interp_linear(vec2_size_t &v)
Perform a linear interpolation of v into the function implied by the tensor and grid.
Definition: tensor_grid.h:1274
void is_valid() const
Check that the o2scl::tensor_grid object is valid.
Definition: tensor_grid.h:231
void copy_table3d_align(size_t ix_x, size_t ix_y, size_vec2_t &index, table3d &tab, std::string slice_name="z")
Create a slice in a o2scl::table3d object with an aligned grid.
Definition: tensor_grid.h:907
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:410
void set_interp_type(size_t interp_type)
Set interpolation type for interpolate()
Definition: tensor_grid.h:1143
bool is_slice(std::string name, size_t &ix) const
Return true if slice is already present.
bool is_xy_set() const
True if the grid has been set.
void set_slice_all(std::string name, double val)
Set all of the values in slice name to val.
generic failure
Definition: err_hnd.h:61
Rank 3 tensor with a grid.
Definition: tensor_grid.h:1705
void default_grid()
Use a default grid which just uses the index.
Definition: tensor_grid.h:511
void copy_table3d_interp(size_t ix_x, size_t ix_y, size_vec2_t &index, table3d &tab, std::string slice_name="z")
Copy to a slice in a table3d object using interpolation.
Definition: tensor_grid.h:998
size_t lookup_grid(size_t i, double val)
Lookup index for grid closest to val.
Definition: tensor_grid.h:692
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:1714
size_t total_size() const
Returns the size of the tensor (the product of the sizes over every index)
Definition: tensor.h:488
double get_grid_y(size_t iy) const
Get y grid point at index iy.
double interp(double x)
Interpolate x and return the results.
Definition: tensor_grid.h:1624
vec_t & get_data()
Return a reference to the data (for HDF I/O)
Definition: tensor_grid.h:392
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:208
void get_size(size_t &nx, size_t &ny) const
Get the size of the slices.
double interp(double x, double y, double z)
Interpolate (x,y,z) and return the results.
Definition: tensor_grid.h:1755
double get_grid(size_t i, size_t j) const
Lookup jth value on the ith grid.
Definition: tensor_grid.h:620
double interp(double x, double y, double z, double a)
Interpolate (x,y,z,a) and return the results.
Definition: tensor_grid.h:1825
double & get(const size_vec_t &index)
Get the element indexed by index.
Definition: tensor.h:298
double interp_linear(double x, double y)
Interpolate (x,y) and return the results.
Definition: tensor_grid.h:1695
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:2746
tensor_grid2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
Definition: tensor_grid.h:1649
size_t rk
Rank.
Definition: tensor.h:127
double interp(double x, double y)
Interpolate (x,y) and return the results.
Definition: tensor_grid.h:1688
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
The O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl namespace ...
Definition: table.h:53
tensor_grid4()
Create an empty tensor.
Definition: tensor_grid.h:1778
double interp_linear_power_two(vec2_size_t &v)
Perform linear interpolation assuming that all indices can take only two values.
Definition: tensor_grid.h:1325
Rank 4 tensor with a grid.
Definition: tensor_grid.h:1772
void is_valid() const
Check that the o2scl::tensor object is valid.
Definition: tensor.h:177
double eval(std::map< std::string, double > *vars=0)
Evalate the previously compiled expression using variables specified in vars.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void interp_linear_vec(vec2_size_t &v, size_t ifree, vec2_t &res)
Perform a linear interpolation of v into the tensor leaving one index free resulting in a vector...
Definition: tensor_grid.h:1483
void set_grid_packed(const vec2_t &grid_vec)
Set the grid.
Definition: tensor_grid.h:472
void interp_linear_vec0(vec2_size_t &v, vec2_t &res)
Perform a linear interpolation of v[1] to v[n-1] resulting in a vector.
Definition: tensor_grid.h:1368
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:372
void set_grid(const vec_vec_t &grid_vecs)
Set grid from a vector of vectors of grid points.
Definition: tensor_grid.h:490
void copy_grid(size_t i, rvec_t &v)
Copy grid for index i to vector v.
Definition: tensor_grid.h:609
void set_grid(size_t i, size_t j, double val)
Set the jth value on the ith grid.
Definition: tensor_grid.h:637
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)
Copy to a slice in a table3d object using interpolation.
Definition: tensor_grid.h:1041
Evaluate a mathematical expression in a string.
double interp_linear(double x, double y, double z, double a)
Interpolate (x,y,z,a) and return the results.
Definition: tensor_grid.h:1832
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor_grid.h:1132
double get_val(const vec2_t &gridp)
Get the element closest to grid point gridp.
Definition: tensor_grid.h:346
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:321
tensor_grid copy_slice_interp(size_vec2_t &ifix, vec2_t &vals)
Copy an abitrary slice by fixing 1 or more indices and use interpolation to return a new tensor_grid ...
Definition: tensor_grid.h:764
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.
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:529
Searching class for monotonic data with caching.
Definition: search_vec.h:74
size_t itype
Interpolation type.
Definition: tensor_grid.h:169
A data structure containing one or more slices of two-dimensional data points defined on a grid...
Definition: table3d.h:77
Store data in an O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$sc...
Definition: hdf_file.h:101
void new_slice(std::string name)
Add a new slice.
bool grid_set
If true, the grid has been set by the user.
Definition: tensor_grid.h:166
tensor_grid2()
Create an empty tensor.
Definition: tensor_grid.h:1646
void set(const size_vec_t &index, double val)
Set the element indexed by index to value val.
Definition: tensor.h:246
vec_t grid
A rank-sized set of arrays for the grid points.
Definition: tensor_grid.h:163
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:1781
void set_grid_i_func(size_t ix, std::string func)
Set grid for one index from a function.
Definition: tensor_grid.h:553
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:708
vec_size_t size
A rank-sized vector of the sizes of each dimension.
Definition: tensor.h:124
void set_val(const vec2_t &grdp, double val)
Set the element closest to grid point grdp to value val.
Definition: tensor_grid.h:292
void set_xy(std::string x_name, size_t nx, const vec_t &x, std::string y_name, size_t ny, const vec2_t &y)
Initialize the x-y grid.
Definition: table3d.h:119
tensor_grid3()
Create an empty tensor.
Definition: tensor_grid.h:1711
Linear interpolation (GSL)
Definition: interp.h:210
std::string szttos(size_t x)
Convert a size_t to a string.
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:146
tensor_grid1()
Create an empty tensor.
Definition: tensor_grid.h:1585
double interp_linear(double x)
Interpolate x and return the results.
Definition: tensor_grid.h:1630
Linear.
Definition: interp.h:71
bool is_size_set() const
True if the size of the table has been set.
tensor_grid1(size_t sz)
Create a rank 2 tensor of size (sz,sz2,sz3)
Definition: tensor_grid.h:1588
Tensor class with arbitrary dimensions.
Definition: tensor.h:112
double interp_linear(double x, double y, double z)
Interpolate (x,y,z) and return the results.
Definition: tensor_grid.h:1762

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