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 set_grid() or \ref
75  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  \note Currently, HDF5 I/O is only allowed if the tensor is
92  allocated with std::vector-based types, and the \ref
93  interpolate() function only works with ublas-based vector types.
94 
95  \todo It is possible for the user to create a tensor_grid
96  object, upcast it to a tensor object, and then use
97  tensor::resize() to resize the tensor, failing to resize the
98  grid. This probably needs fixing.
99 
100  \future Is it really necessary that get_data() is public?
101  \future Only allocate space for grid if it is set
102  \future Consider creating a new set_grid() function which
103  takes grids from an object like uniform_grid. Maybe make a
104  constructor for a tensor_grid object which just takes
105  as input a set of grids?
106  */
107  template<class vec_t=std::vector<double>,
108  class vec_size_t=std::vector<size_t> > class tensor_grid :
109  public tensor<vec_t,vec_size_t> {
110 
111  public:
112 
113 #ifndef DOXYGEN_INTERNAL
114 
115  protected:
116 
117 #ifdef O2SCL_NEVER_DEFINED
118  }{
119 #endif
120 
121  /// A rank-sized set of arrays for the grid points
122  vec_t grid;
123 
124  /// If true, the grid has been set by the user
125  bool grid_set;
126 
127  /// Interpolation type
128  size_t itype;
129 
130 #endif
131 
132  public:
133 
134  /// Create an empty tensor with zero rank
135  tensor_grid() : tensor<vec_t,vec_size_t>() {
136  grid_set=false;
137  itype=itp_linear;
138  }
139 
140  /** \brief Create a tensor of rank \c rank with sizes given in \c dim
141 
142  The parameter \c dim must be a vector of sizes with length \c
143  rank. If the user requests any of the sizes to be zero, this
144  constructor will call the error handler.
145  */
146  template<class size_vec_t>
147  tensor_grid(size_t rank, const size_vec_t &dim) :
148  tensor<vec_t,vec_size_t>(rank,dim) {
149  grid_set=false;
150  itype=itp_linear;
151  // Note that the parent method sets rk to be equal to rank
152  for(size_t i=0;i<this->rk;i++) {
153  if (dim[i]==0) {
154  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
155  "rank for index "+szttos(i)+" in tensor_grid::"+
156  "tensor_grid(size_t,size_vec_t)").c_str(),
157  exc_einval);
158  }
159  }
160  }
161 
162  virtual ~tensor_grid() {
163  }
164 
165  /// \name Set functions
166  //@{
167  /// Set the element closest to grid point \c grdp to value \c val
168  template<class vec2_t>
169  void set_val(const vec2_t &grdp, double val) {
170 
171  // Find indices
172  vec_size_t index(this->rk);
173  for(size_t i=0;i<this->rk;i++) {
174  index[i]=lookup_grid(i,grdp[i]);
175  }
176 
177  // Pack
178  size_t ix=index[0];
179  for(size_t i=1;i<this->rk;i++) {
180  ix*=this->size[i];
181  ix+=index[i];
182  }
183 
184  // Set value
185  this->data[ix]=val;
186 
187  return;
188  }
189 
190  /** \brief Set the element closest to grid point \c grdp to value
191  \c val
192 
193  The parameters \c closest and \c grdp may be identical,
194  allowing one to update a vector \c grdp with the
195  closest grid point.
196  */
197  template<class vec2_t, class vec3_t>
198  void set_val(const vec2_t &grdp, double val, vec3_t &closest) {
199 
200  // Find indices
201  vec_size_t index(this->rk);
202  for(size_t i=0;i<this->rk;i++) {
203  index[i]=lookup_grid_val(i,grdp[i],closest[i]);
204  }
205 
206  // Pack
207  size_t ix=index[0];
208  for(size_t i=1;i<this->rk;i++) {
209  ix*=this->size[i];
210  ix+=index[i];
211  }
212 
213  // Set value
214  this->data[ix]=val;
215 
216  return;
217  }
218  //@}
219 
220  /// \name Get functions
221  //@{
222  /// Get the element closest to grid point \c gridp
223  template<class vec2_t> double get_val(const vec2_t &gridp) {
224 
225  // Find indices
226  vec_size_t index(this->rk);
227  for(size_t i=0;i<this->rk;i++) {
228  index[i]=lookup_grid(i,gridp[i]);
229  }
230 
231  // Pack
232  size_t ix=index[0];
233  for(size_t i=1;i<this->rk;i++) {
234  ix*=this->size[i];
235  ix+=index[i];
236  }
237 
238  // Set value
239  return this->data[ix];
240  }
241 
242  /** \brief Get the element closest to grid point \c gridp,
243  store grid values in \c closest and return value
244 
245  The parameters \c gridp and \c closest may refer to the
246  same object.
247  */
248  template<class vec2_t, class vec3_t>
249  double get_val(const vec2_t &gridp, vec3_t &closest) {
250 
251  // Find indices
252  vec_size_t index(this->rk);
253  for(size_t i=0;i<this->rk;i++) {
254  index[i]=lookup_grid_val(i,gridp[i],closest[i]);
255  }
256 
257  // Pack
258  size_t ix=index[0];
259  for(size_t i=1;i<this->rk;i++) {
260  ix*=this->size[i];
261  ix+=index[i];
262  }
263 
264  // Set value
265  return this->data[ix];
266  }
267  //@}
268 
269  /// \name Resize method
270  //@{
271  /** \brief Resize the tensor to rank \c rank with sizes
272  given in \c dim
273 
274  The parameter \c dim must be a vector of sizes with a length
275  equal to \c rank. This resize method is always destructive,
276  and the grid is always reset.
277 
278  If the user requests any of the sizes to be zero, this
279  function will call the error handler.
280  */
281  template<class size_vec2_t>
282  void resize(size_t rank, const size_vec2_t &dim) {
283  // Double check that none of the sizes that the user
284  // specified are zero
285  for(size_t i=0;i<rank;i++) {
286  if (dim[i]==0) {
287  O2SCL_ERR((((std::string)"Requested zero size with non-zero ")+
288  "rank for index "+szttos(i)+" in tensor_grid::"+
289  "resize().").c_str(),exc_einval);
290  }
291  }
292  // Set the new rank
293  this->rk=rank;
294  // Resize the size vector
295  this->size.resize(this->rk);
296  // Reset the grid
297  grid_set=false;
298  grid.resize(0);
299  // If the new rank is zero, reset the data, otherwise,
300  // update the size vector and reset the data vector
301  if (rank==0) {
302  this->data.resize(0);
303  return;
304  } else {
305  size_t tot=1;
306  for(size_t i=0;i<this->rk;i++) {
307  this->size[i]=dim[i];
308  tot*=this->size[i];
309  }
310  this->data.resize(tot);
311  }
312  return;
313  }
314 
315  //@}
316 
317  /// \name Grid manipulation
318  //@{
319  /// Return true if the grid has been set
320  bool is_grid_set() const { return grid_set; }
321 
322  /** \brief Set the grid
323 
324  The grid must be specified for all of the dimensions at
325  once. Denote \f$ (\mathrm{size})_0 \f$ as the size
326  of the first dimension, \f$ (\mathrm{size})_1 \f$ as the
327  size of the second dimesion, and so on. Then the
328  first \f$ (\mathrm{size})_0 \f$ entries in \c grid must
329  be the grid for the first dimension, the next
330  \f$ (\mathrm{size})_1 \f$ entries must be the grid
331  for the second dimension, and so on. Thus \c grid
332  must be a vector of size
333  \f[
334  \sum_{i=0}^{\mathrm{rank}} (\mathrm{size})_i
335  \f]
336 
337  Note that the grid is copied so the function argument may
338  be destroyed by the user after calling set_grid_packed() without
339  affecting the tensor grid.
340 
341  \future Define a more generic interface for matrix types
342  */
343  template<class vec2_t>
344  void set_grid_packed(const vec2_t &grid_vec) {
345  if (this->rk==0) {
346  O2SCL_ERR2("Tried to set grid for empty tensor in ",
347  "tensor_grid::set_grid_packed().",exc_einval);
348  }
349  size_t ngrid=0;
350  for(size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
351  grid.resize(ngrid);
352  for(size_t i=0;i<ngrid;i++) {
353  grid[i]=grid_vec[i];
354  }
355  grid_set=true;
356  return;
357  }
358 
359  /** \brief Set grid from a vector of vectors of grid points
360  */
361  template<class vec_vec_t>
362  void set_grid(const vec_vec_t &grid_vecs) {
363  if (this->rk==0) {
364  O2SCL_ERR2("Tried to set grid for empty tensor in ",
365  "tensor_grid::set_grid().",exc_einval);
366  }
367  size_t ngrid=0;
368  for(size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
369  grid.resize(ngrid);
370  size_t k=0;
371  for(size_t i=0;i<this->rk;i++) {
372  for(size_t j=0;j<this->size[i];j++) {
373  grid[k]=grid_vecs[i][j];
374  k++;
375  }
376  }
377  grid_set=true;
378  return;
379  }
380 
381  /** \brief Copy grid for index \c i to vector \c v
382 
383  The type \c rvec_t must be a vector with a resize
384  method.
385  */
386  template<class rvec_t> void copy_grid(size_t i, rvec_t &v) {
387  v.resize(this->size[i]);
388  size_t istart=0;
389  for(size_t k=0;k<i;k++) istart+=this->size[k];
390  for(size_t k=0;k<this->size[i];k++) {
391  v[k]=grid[istart+k];
392  }
393  return;
394  }
395 
396  /// Lookup jth value on the ith grid
397  double get_grid(size_t i, size_t j) const {
398  if (!grid_set) {
399  O2SCL_ERR("Grid not set in tensor_grid::get_grid().",
400  exc_einval);
401  }
402  if (i>=this->rk) {
403  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
404  " greater than or equal to rank, "+szttos(this->rk)+
405  ", in tensor_grid::get_grid().").c_str(),
406  exc_einval);
407  }
408  size_t istart=0;
409  for(size_t k=0;k<i;k++) istart+=this->size[k];
410  return grid[istart+j];
411  }
412 
413  /// Set the jth value on the ith grid
414  void set_grid(size_t i, size_t j, double val) {
415  if (!grid_set) {
416  O2SCL_ERR("Grid not set in tensor_grid::get_grid().",
417  exc_einval);
418  }
419  if (i>=this->rk) {
420  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
421  " greater than or equal to rank, "+szttos(this->rk)+
422  ", in tensor_grid::get_grid().").c_str(),
423  exc_einval);
424  }
425  size_t istart=0;
426  for(size_t k=0;k<i;k++) istart+=this->size[k];
427  grid[istart+j]=val;
428  return;
429  }
430 
431  /// Lookup index for grid closest to \c val
432  size_t lookup_grid(size_t i, double val) {
433  if (!grid_set) {
434  O2SCL_ERR("Grid not set in tensor_grid::lookup_grid().",
435  exc_einval);
436  }
437  if (i>=this->rk) {
438  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
439  " greater than or equal to rank, "+szttos(this->rk)+
440  ", in tensor_grid::lookup_grid().").c_str(),
441  exc_einval);
442  }
443  size_t istart=0;
444 
445  for(size_t j=0;j<i;j++) {
446  istart+=this->size[j];
447  }
448  size_t best=istart;
449  double min=fabs(grid[istart]-val);
450  for(size_t j=istart;j<istart+this->size[i];j++) {
451  if (fabs(grid[j]-val)<min) {
452  best=j;
453  min=fabs(grid[j]-val);
454  }
455  }
456  return best-istart;
457  }
458 
459  /** \brief Lookup indices for grid closest point to \c vals
460 
461  The values in \c vals are not modified by this function.
462 
463  \comment
464  This function must have a different name than
465  lookup_grid() because the template types cause
466  confusion between the two functions.
467  \endcomment
468  */
469  template<class vec2_t, class size_vec2_t>
470  void lookup_grid_vec(const vec2_t &vals, size_vec2_t &indices) const {
471  for(size_t k=0;k<this->rk;k++) {
472  indices[k]=lookup_grid(k,vals[k]);
473  }
474  return;
475  }
476 
477  /** \brief Lookup index for grid closest to \c val, returning the
478  grid point
479 
480  The parameters \c val and \c val2 may refer to the
481  same object.
482  */
483  size_t lookup_grid_val(size_t i, const double &val, double &val2) {
484  if (i>=this->rk) {
485  O2SCL_ERR((((std::string)"Index ")+szttos(i)+
486  " greater than or equal to rank, "+szttos(this->rk)+
487  ", in tensor_grid::lookup_grid_val().").c_str(),
488  exc_einval);
489  }
490  if (grid_set==false) {
491  O2SCL_ERR("Grid not set in tensor_grid::lookup_grid_val().",
492  exc_einval);
493  }
494 
495  // We have to store the grid point in a temporary in case
496  // the parameters gridp and closest refer to the same object
497  double temp=val;
498 
499  size_t istart=0;
500  for(size_t j=0;j<i;j++) istart+=this->size[j];
501  size_t best=istart;
502  double min=fabs(grid[istart]-temp);
503  val2=grid[istart];
504  for(size_t j=istart;j<istart+this->size[i];j++) {
505  if (fabs(grid[j]-temp)<min) {
506  best=j;
507  min=fabs(grid[j]-temp);
508  val2=grid[j];
509  }
510  }
511  return best-istart;
512  }
513 
514  /// Lookup index for grid closest to \c val
515  size_t lookup_grid_packed(size_t i, double val) {
516  if (!grid_set) {
517  O2SCL_ERR("Grid not set in tensor_grid::lookup_grid_packed().",
518  exc_einval);
519  }
520  if (i>=this->rk) {
521  O2SCL_ERR((((std::string)"Index ")+szttos(i)+" greater than rank, "+
522  szttos(this->rk)+
523  ", in tensor_grid::lookup_grid_packed().").c_str(),
524  exc_einval);
525  }
526  size_t istart=0;
527  for(size_t j=0;j<i;j++) istart+=this->size[j];
528  size_t best=istart;
529  double min=fabs(grid[istart]-val);
530  for(size_t j=istart;j<istart+this->size[i];j++) {
531  if (fabs(grid[j]-val)<min) {
532  best=j;
533  min=fabs(grid[j]-val);
534  }
535  }
536  return best;
537  }
538 
539  /// Lookup index for grid closest to \c val
540  size_t lookup_grid_packed_val(size_t i, double val, double &val2) {
541  if (!grid_set) {
542  O2SCL_ERR("Grid not set in tensor_grid::lookup_grid_packed().",
543  exc_einval);
544  }
545  if (i>=this->rk) {
546  O2SCL_ERR((((std::string)"Index ")+szttos(i)+" greater than rank, "+
547  szttos(this->rk)+
548  ", in tensor_grid::lookup_grid_packed().").c_str(),
549  exc_einval);
550  }
551  size_t istart=0;
552  for(size_t j=0;j<i;j++) istart+=this->size[j];
553  size_t best=istart;
554  double min=fabs(grid[istart]-val);
555  val2=grid[istart];
556  for(size_t j=istart;j<istart+this->size[i];j++) {
557  if (fabs(grid[j]-val)<min) {
558  best=j;
559  min=fabs(grid[j]-val);
560  val2=grid[j];
561  }
562  }
563  return best;
564  }
565  //@}
566 
567  /// Return a reference to the data (for HDF I/O)
568  vec_t &get_data() {
569  return this->data;
570  }
571 
572  /// \name Slicing
573  //@{
574  /** \brief Create a slice in a table3d object with an aligned
575  grid
576 
577  This function uses the grid associated with indices \c ix_x
578  and \c ix_y, to copy data to a slice named \c slice_name in
579  the table3d object \c tab . All other indices are fixed
580  to values specified by the user in \c index and the values
581  of <tt>index[ix_x]</tt> and <tt>index[ix_y]</tt> are
582  used for temporary storage.
583 
584  If the table3d object does not currently have a grid set, then
585  the grid is automatically set to be the same as that stored in
586  the tensor_grid object associated with ranks \c ix_x and \c
587  iy_y. If the \ref o2scl::table3d object does have a grid set,
588  then the values returned by \ref o2scl::table3d::get_nx() and
589  \ref o2scl::table3d::get_ny() must be equal to the size of the
590  tensor in indices \c ix_x and ix_y, respectively. If a
591  slice named \c slice_name is not already present in
592  \c tab, then a new slice with that name is created.
593 
594  The error handler is called if \c ix_x is the same as
595  \c ix_y, or if either of these two values is greater
596  than or equal to the tensor rank.
597  */
598  template<class size_vec2_t>
599  void copy_slice_align(size_t ix_x, size_t ix_y, size_vec2_t &index,
600  table3d &tab, std::string slice_name) {
601 
602  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
603  O2SCL_ERR2("Either indices greater than rank or x and y ind",
604  "ices equal in tensor_grid::copy_slice_align().",
605  exc_efailed);
606  }
607 
608  // Get current table3d grid
609  size_t nx, ny;
610  tab.get_size(nx,ny);
611 
612  if (nx==0 && ny==0) {
613 
614  // If there's no grid, just create it
615  std::vector<double> gx, gy;
616  for(size_t i=0;i<this->size[ix_x];i++) {
617  gx.push_back(this->get_grid(ix_x,i));
618  }
619  for(size_t i=0;i<this->size[ix_y];i++) {
620  gy.push_back(this->get_grid(ix_y,i));
621  }
622  nx=gx.size();
623  ny=gy.size();
624  tab.set_xy("x",nx,gx,"y",ny,gy);
625  }
626 
627  // Check that the grids are commensurate
628  if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
629  O2SCL_ERR2("Grids not commensurate in ",
630  "tensor_grid::copy_slice_align().",exc_einval);
631  }
632 
633  // Create slice if not already present
634  size_t is;
635  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
636 
637  // Copy over data
638  for(size_t i=0;i<nx;i++) {
639  for(size_t j=0;j<ny;j++) {
640  index[ix_x]=i;
641  index[ix_y]=j;
642  double val=this->get(index);
643  tab.set(i,j,slice_name,val);
644  }
645  }
646 
647  return;
648  }
649 
650  /** \brief Copy to a slice in a table3d object using interpolation
651 
652  This function uses the grid associated with indices \c ix_x
653  and \c ix_y, and the tensor interpolation function to copy the
654  tensor information to the slice named \c slice_name in the
655  table3d object \c tab . All other indices are fixed
656  to values specified by the user in \c index and the values
657  of <tt>index[ix_x]</tt> and <tt>index[ix_y]</tt> are
658  used for temporary storage.
659 
660  If a slice named \c slice_name is not already present in \c
661  tab, then a new slice with that name is created.
662 
663  The error handler is called if \c ix_x is the same as
664  \c ix_y, or if either of these two values is greater
665  than or equal to the tensor rank.
666 
667  \note This function uses the \ref tensor_grid::interp_linear()
668  for the interpolation.
669  */
670  template<class size_vec2_t>
671  void copy_slice_interp(size_t ix_x, size_t ix_y, size_vec2_t &index,
672  table3d &tab, std::string slice_name) {
673 
674  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
675  O2SCL_ERR2("Either indices greater than rank or x and y ",
676  "indices equal in tensor_grid::copy_slice_interp().",
677  exc_efailed);
678  }
679 
680  // Get current table3d grid
681  size_t nx, ny;
682  tab.get_size(nx,ny);
683 
684  if (nx==0 && ny==0) {
685  // If there's no grid, then just use the aligned version
686  return copy_slice_align(ix_x,ix_y,index,tab,slice_name);
687  }
688 
689  // Create vector of values to interpolate with
690  std::vector<double> vals(this->rk);
691  for(size_t i=0;i<this->rk;i++) {
692  if (i!=ix_x && i!=ix_y) vals[i]=this->get_grid(i,index[i]);
693  }
694 
695  // Create slice if not already present
696  size_t is;
697  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
698 
699  // Loop through the table grid to perform the interpolation
700  for(size_t i=0;i<nx;i++) {
701  for(size_t j=0;j<ny;j++) {
702  vals[ix_x]=tab.get_grid_x(i);
703  vals[ix_y]=tab.get_grid_y(j);
704  tab.set(i,j,slice_name,this->interp_linear(vals));
705  }
706  }
707 
708  return;
709  }
710 
711  /** \brief Copy to a slice in a table3d object using interpolation
712  */
713  template<class vec2_t>
714  void copy_slice_interp_values(size_t ix_x, size_t ix_y,
715  vec2_t &values, table3d &tab,
716  std::string slice_name) {
717 
718  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
719  O2SCL_ERR2("Either indices greater than rank or x and y ",
720  "indices equal in tensor_grid::copy_slice_interp().",
721  exc_efailed);
722  }
723  if (values.size()!=this->rk) {
724  O2SCL_ERR2("Values array not equal to rank ",
725  "in tensor_grid::copy_slice_interp_values().",
726  exc_efailed);
727  }
728 
729  // Get current table3d grid
730  size_t nx, ny;
731  tab.get_size(nx,ny);
732 
733  //if (nx==0 && ny==0) {
734  // If there's no grid, then just use the aligned version
735  //return copy_slice_align(ix_x,ix_y,index,tab,slice_name);
736  //}
737 
738  // Create slice if not already present
739  size_t is;
740  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
741 
742  // Loop through the table grid to perform the interpolation
743  for(size_t i=0;i<nx;i++) {
744  for(size_t j=0;j<ny;j++) {
745  values[ix_x]=tab.get_grid_x(i);
746  values[ix_y]=tab.get_grid_y(j);
747  tab.set(i,j,slice_name,this->interp_linear(values));
748  }
749  }
750 
751  return;
752  }
753  //@}
754 
755  /// \name Clear method
756  //@{
757  /// Clear the tensor of all data and free allocated memory
758  void clear() {
759  grid.resize(0);
760  grid_set=false;
762  return;
763  }
764  //@}
765 
766  /// \name Interpolation
767  //@{
768  /// Set interpolation type for \ref interpolate()
769  void set_interp_type(size_t interp_type) {
770  itype=interp_type;
771  return;
772  }
773 
774  /** \brief Interpolate values \c vals into the tensor,
775  returning the result
776 
777  This is a quick and dirty implementation of n-dimensional
778  interpolation by recursive application of the 1-dimensional
779  routine from \ref interp_vec, using the base interpolation
780  object specified in the template parameter \c base_interp_t.
781  This will be very slow for sufficiently large data sets.
782 
783  \note This function requires a range objects to obtain ranges
784  of vector objects. In ublas, this is done with
785  <tt>ublas::vector_range</tt> objects, so this function will
786  certainly work for \ref tensor_grid objects built on ublas
787  vector types. There is no corresponding <tt>std::range</tt>,
788  but you may be able to use either <tt>ublas::vector_range</tt>
789  or <tt>Boost.Range</tt> in order to call this function with
790  \ref tensor_grid objects built on <tt>std::vector</tt>.
791  However, this is not fully tested at the moment.
792 
793  \future It should be straightforward to improve the scaling of
794  this algorithm significantly by creating a "window" of local
795  points around the point of interest. This could be done easily
796  by constructing an initial subtensor. However, this should
797  probably be superceded by a more generic alternative which
798  avoids explicit use of the 1-d interpolation types.
799  */
800  template<class range_t=ub_range,
801  class data_range_t=ubvector_range,
802  class index_range_t=ubvector_size_t_range>
803  double interpolate(double *vals) {
804 
805  typedef interp_vec<vec_t> interp_t;
806 
807  if (this->rk==1) {
808 
809  interp_t si(this->size[0],grid,this->data,itype);
810  return si.eval(vals[0]);
811 
812  } else {
813 
814  // Get total number of interpolations at this level
815  size_t ss=1;
816  for(size_t i=1;i<this->rk;i++) ss*=this->size[i];
817 
818  // Create space for y vectors and interpolators
819  std::vector<vec_t> yvec(ss);
820  std::vector<interp_t *> si(ss);
821  for(size_t i=0;i<ss;i++) {
822  yvec[i].resize(this->size[0]);
823  }
824 
825  // Create space for interpolation results
826  tensor_grid tdat;
827  index_range_t size_new(this->size,ub_range(1,this->rk));
828  tdat.resize(this->rk-1,size_new);
829 
830  // Set grid for temporary tensor
831  data_range_t grid_new(grid,ub_range(this->size[0],grid.size()));
832  tdat.set_grid_packed(grid_new);
833 
834  // Create starting coordinate and counter
835  vec_size_t co(this->rk);
836  for(size_t i=0;i<this->rk;i++) co[i]=0;
837  size_t cnt=0;
838 
839  // Loop over every interpolation
840  bool done=false;
841  while(done==false) {
842 
843  // Fill yvector with the appropriate data
844  for(size_t i=0;i<this->size[0];i++) {
845  co[0]=i;
846  yvec[cnt][i]=this->get(co);
847  }
848 
849  si[cnt]=new interp_t(this->size[0],grid,yvec[cnt],itype);
850 
851  index_range_t co2(co,ub_range(1,this->rk));
852  tdat.set(co2,si[cnt]->eval(vals[0]));
853 
854  // Go to next interpolation
855  cnt++;
856  co[this->rk-1]++;
857  // carry if necessary
858  for(int j=((int)this->rk)-1;j>0;j--) {
859  if (co[j]>=this->size[j]) {
860  co[j]=0;
861  co[j-1]++;
862  }
863  }
864 
865  // Test if done
866  if (cnt==ss) done=true;
867 
868  // End of while loop
869  }
870 
871  // Now call the next level of interpolation
872  double res=tdat.interpolate(vals+1);
873 
874  for(size_t i=0;i<ss;i++) {
875  delete si[i];
876  }
877 
878  return res;
879  }
880  }
881 
882  /** \brief Perform a linear interpolation of \c v into the
883  function implied by the tensor and grid
884 
885  This performs multi-dimensional linear interpolation (or
886  extrapolation) It works by first using \ref o2scl::search_vec
887  to find the interval containing (or closest to) the specified
888  point in each direction and constructing the corresponding
889  hypercube of size \f$ 2^{\mathrm{rank}} \f$ containing \c v.
890  It then calls \ref interp_linear_power_two() to perform the
891  interpolation in that hypercube.
892 
893  \future This starts with a small copy, which can be eliminated
894  by creating a new version of interp_linear_power_two
895  which accepts an offset vector parameter so that the
896  first interpolation is offset. Remaining interpolations
897  don't need to be offset because the tensor has to be
898  created from the previous interpolation round.
899  */
900  template<class vec2_size_t> double interp_linear(vec2_size_t &v) {
901 
902  // Find the the corner of the hypercube containing v
903  size_t rgs=0;
904  std::vector<size_t> loc(this->rk);
905  std::vector<double> gnew;
906  for(size_t i=0;i<this->rk;i++) {
907  std::vector<double> grid_unpacked(this->size[i]);
908  for(size_t j=0;j<this->size[i];j++) {
909  grid_unpacked[j]=grid[j+rgs];
910  }
911  search_vec<std::vector<double> > sv(this->size[i],grid_unpacked);
912  loc[i]=sv.find(v[i]);
913  gnew.push_back(grid_unpacked[loc[i]]);
914  gnew.push_back(grid_unpacked[loc[i]+1]);
915  rgs+=this->size[i];
916  }
917 
918  // Now construct a 2^{rk}-sized tensor containing only that
919  // hypercube
920  std::vector<size_t> snew(this->rk);
921  for(size_t i=0;i<this->rk;i++) {
922  snew[i]=2;
923  }
924  tensor_grid tnew(this->rk,snew);
925  tnew.set_grid_packed(gnew);
926 
927  // Copy over the relevant data
928  for(size_t i=0;i<tnew.total_size();i++) {
929  std::vector<size_t> index_new(this->rk), index_old(this->rk);
930  tnew.unpack_indices(i,index_new);
931  for(size_t j=0;j<this->rk;j++) index_old[j]=index_new[j]+loc[j];
932  tnew.set(index_new,this->get(index_old));
933  }
934 
935  // Now use interp_power_two()
936  return tnew.interp_linear_power_two(v);
937  }
938 
939  /** \brief Perform linear interpolation assuming that all
940  indices can take only two values
941 
942  This function works by recursively slicing the hypercube of
943  size \f$ 2^{\mathrm{rank}} \f$ into a hypercube of size \f$
944  2^{\mathrm{rank-1}} \f$ performing linear interpolation for
945  each pair of points.
946 
947  \note This is principally a function for internal use
948  by \ref interp_linear().
949  */
950  template<class vec2_size_t>
951  double interp_linear_power_two(vec2_size_t &v) {
952 
953  if (this->rk==1) {
954  return this->data[0]+(this->data[1]-this->data[0])/
955  (grid[1]-grid[0])*(v[0]-grid[0]);
956  }
957 
958  size_t last=this->rk-1;
959  double frac=(v[last]-get_grid(last,0))/
960  (get_grid(last,1)-get_grid(last,0));
961 
962  // Create new size vector and grid
963  tensor_grid tnew(this->rk-1,this->size);
964  tnew.set_grid_packed(grid);
965 
966  // Create data in new tensor, removing the last index through
967  // linear interpolation
968  for(size_t i=0;i<tnew.total_size();i++) {
969  std::vector<size_t> index(this->rk);
970  tnew.unpack_indices(i,index);
971  index[this->rk-1]=0;
972  double val_lo=this->get(index);
973  index[this->rk-1]=1;
974  double val_hi=this->get(index);
975  tnew.set(index,val_lo+frac*(val_hi-val_lo));
976  }
977 
978  // Recursively interpolate the smaller tensor
979  return tnew.interp_linear_power_two(v);
980  }
981 
982  /** \brief Perform a linear interpolation of <tt>v[1]</tt>
983  to <tt>v[n-1]</tt> resulting in a vector
984 
985  \note The type <tt>vec2_t</tt> for the vector <tt>res</tt>
986  must have a <tt>resize()</tt> method.
987 
988  This performs multi-dimensional linear interpolation (or
989  extrapolation) in the last <tt>n-1</tt> indices of the
990  rank-<tt>n</tt> tensor leaving the first index free and places
991  the results in the vector \c res.
992  */
993  template<class vec2_size_t, class vec2_t>
994  void interp_linear_vec0(vec2_size_t &v, vec2_t &res) {
995 
996  // Find the the corner of the hypercube containing v
997  size_t rgs=0;
998  std::vector<size_t> loc(this->rk);
999  std::vector<double> gnew;
1000  for(size_t i=0;i<this->size[0];i++) {
1001  gnew.push_back(grid[i]);
1002  }
1003  rgs=this->size[0];
1004  loc[0]=0;
1005  for(size_t i=1;i<this->rk;i++) {
1006  std::vector<double> grid_unpacked(this->size[i]);
1007  for(size_t j=0;j<this->size[i];j++) {
1008  grid_unpacked[j]=grid[j+rgs];
1009  }
1010  search_vec<std::vector<double> > sv(this->size[i],grid_unpacked);
1011  loc[i]=sv.find(v[i]);
1012  gnew.push_back(grid_unpacked[loc[i]]);
1013  gnew.push_back(grid_unpacked[loc[i]+1]);
1014  rgs+=this->size[i];
1015  }
1016 
1017  // Now construct a n*2^{rk-1}-sized tensor containing only that
1018  // hypercube
1019  std::vector<size_t> snew(this->rk);
1020  snew[0]=this->size[0];
1021  for(size_t i=1;i<this->rk;i++) {
1022  snew[i]=2;
1023  }
1024  tensor_grid tnew(this->rk,snew);
1025  tnew.set_grid_packed(gnew);
1026 
1027  // Copy over the relevant data
1028  for(size_t i=0;i<tnew.total_size();i++) {
1029  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1030  tnew.unpack_indices(i,index_new);
1031  for(size_t j=0;j<this->rk;j++) {
1032  index_old[j]=index_new[j]+loc[j];
1033  }
1034  tnew.set(index_new,this->get(index_old));
1035  }
1036 
1037  // Now use interp_power_two_vec0()
1038  tnew.interp_linear_power_two_vec0(v,res);
1039 
1040  return;
1041  }
1042 
1043  /** \brief Perform linear interpolation assuming that the last
1044  <tt>n-1</tt> indices can take only two values
1045 
1046  This function performs linear interpolation assuming that the
1047  last <tt>n-1</tt> indices can take only two values and placing
1048  the result into <tt>res</tt>.
1049 
1050  \note The type <tt>vec2_t</tt> for the vector <tt>res</tt>
1051  must have a <tt>resize()</tt> method. This is principally a
1052  function for internal use by \ref interp_linear_vec0().
1053  */
1054  template<class vec2_size_t, class vec2_t>
1055  void interp_linear_power_two_vec0(vec2_size_t &v, vec2_t &res) {
1056 
1057  if (this->rk==2) {
1058  size_t n=this->size[0];
1059  res.resize(n);
1060  vec_size_t ix0(2), ix1(2);
1061  ix0[1]=0;
1062  ix1[1]=1;
1063  for(size_t i=0;i<n;i++) {
1064  ix0[0]=i;
1065  ix1[0]=i;
1066  res[i]=this->get(ix0)+(this->get(ix1)-this->get(ix0))/
1067  (grid[n+1]-grid[n])*(v[1]-grid[n]);
1068  }
1069  return;
1070  }
1071 
1072  size_t last=this->rk-1;
1073  double frac=(v[last]-get_grid(last,0))/
1074  (get_grid(last,1)-get_grid(last,0));
1075 
1076  // Create new size vector and grid
1077  tensor_grid tnew(this->rk-1,this->size);
1078  tnew.set_grid_packed(grid);
1079 
1080  // Create data in new tensor, removing the last index through
1081  // linear interpolation
1082  for(size_t i=0;i<tnew.total_size();i++) {
1083  std::vector<size_t> index(this->rk);
1084  tnew.unpack_indices(i,index);
1085  index[this->rk-1]=0;
1086  double val_lo=this->get(index);
1087  index[this->rk-1]=1;
1088  double val_hi=this->get(index);
1089  tnew.set(index,val_lo+frac*(val_hi-val_lo));
1090  }
1091 
1092  // Recursively interpolate the smaller tensor
1093  tnew.interp_linear_power_two_vec0(v,res);
1094 
1095  return;
1096  }
1097 
1098  /** \brief Perform a linear interpolation of <tt>v</tt>
1099  into the tensor leaving one index free resulting in a vector
1100 
1101  This performs multi-dimensional linear interpolation (or
1102  extrapolation) in the last <tt>n-1</tt> indices of the
1103  rank-<tt>n</tt> tensor leaving the first index free and places
1104  the results in the vector \c res.
1105 
1106  \future This function could be more efficient.
1107  */
1108  template<class vec2_size_t, class vec2_t>
1109  void interp_linear_vec(vec2_size_t &v, size_t ifree, vec2_t &res) {
1110 
1111  size_t n=this->size[ifree];
1112 
1113  // This function uses interp_linear_power_two_vec0(), so it
1114  // works by remapping the indices. This defines the remapping.
1115  std::vector<size_t> map;
1116  map.push_back(ifree);
1117  for(size_t i=0;i<this->rk;i++) {
1118  if (i!=ifree) {
1119  map.push_back(i);
1120  }
1121  }
1122 
1123  // Find the the corner of the hypercube containing v
1124  size_t rgs=0;
1125  vec_size_t loc(this->rk);
1126  loc[ifree]=0;
1127  for(size_t i=0;i<this->rk;i++) {
1128  vec_t grid_unpacked(this->size[i]);
1129  for(size_t j=0;j<this->size[i];j++) {
1130  grid_unpacked[j]=grid[j+rgs];
1131  }
1132  search_vec<vec_t> sv(this->size[i],grid_unpacked);
1133  if (i!=ifree) {
1134  loc[i]=sv.find(v[i]);
1135  }
1136  rgs+=this->size[i];
1137  }
1138 
1139  // Compute the remapped grid and interpolating vector
1140  std::vector<double> gnew, vnew;
1141  for(size_t new_ix=0;new_ix<this->rk;new_ix++) {
1142  for(size_t old_ix=0;old_ix<this->rk;old_ix++) {
1143  if (map[new_ix]==old_ix) {
1144  vnew.push_back(v[old_ix]);
1145  if (old_ix==ifree) {
1146  for(size_t j=0;j<this->size[old_ix];j++) {
1147  gnew.push_back(this->get_grid(old_ix,j));
1148  }
1149  } else {
1150  gnew.push_back(this->get_grid(old_ix,loc[old_ix]));
1151  gnew.push_back(this->get_grid(old_ix,loc[old_ix]+1));
1152  }
1153  }
1154  }
1155  }
1156 
1157  // Now construct a n*2^{rk-1}-sized tensor containing only the
1158  // hypercube needed to do the interpolation
1159 
1160  // Specify the size of each rank
1161  std::vector<size_t> snew;
1162  snew.push_back(n);
1163  for(size_t i=0;i<this->rk;i++) {
1164  if (i!=ifree) {
1165  snew.push_back(2);
1166  }
1167  }
1168 
1169  // Create the tensor and set the grid
1170  tensor_grid tnew(this->rk,snew);
1171  tnew.set_grid_packed(gnew);
1172 
1173  // Copy over the relevant data
1174  for(size_t i=0;i<tnew.total_size();i++) {
1175  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1176  tnew.unpack_indices(i,index_new);
1177  for(size_t j=0;j<this->rk;j++) {
1178  index_old[map[j]]=index_new[j]+loc[map[j]];
1179  }
1180  tnew.set(index_new,this->get(index_old));
1181  }
1182 
1183  // Now use interp_power_two_vec()
1184  tnew.interp_linear_power_two_vec0(vnew,res);
1185 
1186  return;
1187  }
1188  //@}
1189 
1190  template<class vecf_t, class vecf_size_t> friend void o2scl_hdf::hdf_output
1192  std::string name);
1193 
1194  template<class vecf_t, class vecf_size_t> friend void o2scl_hdf::hdf_input
1196  std::string name);
1197 
1198  };
1199 
1200  /** \brief Rank 1 tensor with a grid
1201 
1202  \future Make rank-specific get_val and set_val functions?
1203  */
1204  template<class vec_t=std::vector<double>,
1205  class vec_size_t=std::vector<size_t> > class tensor_grid1 :
1206  public tensor_grid<vec_t,vec_size_t> {
1207 
1208  public:
1209 
1210  /// Create an empty tensor
1211  tensor_grid1() : tensor_grid<vec_t,vec_size_t>() {}
1212 
1213  /// Create a rank 2 tensor of size \c (sz,sz2,sz3)
1214  tensor_grid1(size_t sz) : tensor_grid<vec_t,vec_size_t>() {
1215  this->rk=1;
1216  this->size.resize(1);
1217  this->size[0]=sz;
1218  this->data.resize(sz);
1219  this->grid_set=false;
1220  }
1221 
1222 #ifdef O2SCL_NEVER_DEFINED
1223  }{
1224 #endif
1225 
1226  virtual ~tensor_grid1() {
1227  }
1228 
1229  /// Get the element indexed by \c (ix1)
1230  double &get(size_t ix1) {
1231  size_t sz[1]={ix1};
1233  }
1234 
1235  /// Get the element indexed by \c (ix1)
1236  const double &get(size_t ix1) const {
1237  size_t sz[1]={ix1};
1239  }
1240 
1241  /// Set the element indexed by \c (ix1) to value \c val
1242  void set(size_t ix1, double val) {
1243  size_t sz[1]={ix1};
1245  }
1246 
1247  /// Interpolate \c x and return the results
1248  template<class range_t=ub_range, class data_range_t=ubvector_range,
1249  class index_range_t=ubvector_size_t_range>
1250  double interp(double x) {
1251  return tensor_grid<vec_t,vec_size_t>::template interpolate
1252  <range_t,data_range_t,index_range_t>(&x);
1253  }
1254 
1255  /// Interpolate \c x and return the results
1256  double interp_linear(double x) {
1257  double arr[1]={x};
1259  }
1260 
1261  };
1262 
1263  /** \brief Rank 2 tensor with a grid
1264  */
1265  template<class vec_t=std::vector<double>,
1266  class vec_size_t=std::vector<size_t> > class tensor_grid2 :
1267  public tensor_grid<vec_t,vec_size_t> {
1268 
1269  public:
1270 
1271  /// Create an empty tensor
1272  tensor_grid2() : tensor_grid<vec_t,vec_size_t>() {}
1273 
1274  /// Create a rank 2 tensor of size \c (sz,sz2)
1275  tensor_grid2(size_t sz, size_t sz2) : tensor_grid<vec_t,vec_size_t>() {
1276  this->rk=2;
1277  this->size.resize(2);
1278  this->size[0]=sz;
1279  this->size[1]=sz2;
1280  size_t tot=sz*sz2;
1281  this->data.resize(tot);
1282  this->grid_set=false;
1283  }
1284 
1285 #ifdef O2SCL_NEVER_DEFINED
1286  }{
1287 #endif
1288 
1289  virtual ~tensor_grid2() {
1290  }
1291 
1292  /// Get the element indexed by \c (ix1,ix2)
1293  double &get(size_t ix1, size_t ix2) {
1294  size_t sz[2]={ix1,ix2};
1296  }
1297 
1298  /// Get the element indexed by \c (ix1,ix2)
1299  const double &get(size_t ix1, size_t ix2) const {
1300  size_t sz[2]={ix1,ix2};
1302  }
1303 
1304  /// Set the element indexed by \c (ix1,ix2) to value \c val
1305  void set(size_t ix1, size_t ix2, double val) {
1306  size_t sz[2]={ix1,ix2};
1308  return;
1309  }
1310 
1311  /// Interpolate \c (x,y) and return the results
1312  template<class range_t=ub_range, class data_range_t=ubvector_range,
1313  class index_range_t=ubvector_size_t_range>
1314  double interp(double x, double y) {
1315  double arr[2]={x,y};
1316  return tensor_grid<vec_t,vec_size_t>::template interpolate
1317  <range_t,data_range_t,index_range_t>(arr);
1318  }
1319 
1320  /// Interpolate \c (x,y) and return the results
1321  double interp_linear(double x, double y) {
1322  double arr[2]={x,y};
1324  }
1325 
1326  };
1327 
1328  /** \brief Rank 3 tensor with a grid
1329  */
1330  template<class vec_t=std::vector<double>,
1331  class vec_size_t=std::vector<size_t> > class tensor_grid3 :
1332  public tensor_grid<vec_t,vec_size_t> {
1333 
1334  public:
1335 
1336  /// Create an empty tensor
1337  tensor_grid3() : tensor_grid<vec_t,vec_size_t>() {}
1338 
1339  /// Create a rank 3 tensor of size \c (sz,sz2,sz3)
1340  tensor_grid3(size_t sz, size_t sz2, size_t sz3) :
1341  tensor_grid<vec_t,vec_size_t>() {
1342  this->rk=3;
1343  this->size.resize(3);
1344  this->size[0]=sz;
1345  this->size[1]=sz2;
1346  this->size[2]=sz3;
1347  size_t tot=sz*sz2*sz3;
1348  this->data.resize(tot);
1349  this->grid_set=false;
1350  }
1351 
1352 #ifdef O2SCL_NEVER_DEFINED
1353  }{
1354 #endif
1355 
1356  virtual ~tensor_grid3() {
1357  }
1358 
1359  /// Get the element indexed by \c (ix1,ix2,ix3)
1360  double &get(size_t ix1, size_t ix2, size_t ix3) {
1361  size_t sz[3]={ix1,ix2,ix3};
1363  }
1364 
1365  /// Get the element indexed by \c (ix1,ix2,ix3)
1366  const double &get(size_t ix1, size_t ix2, size_t ix3) const {
1367  size_t sz[3]={ix1,ix2,ix3};
1369  }
1370 
1371  /// Set the element indexed by \c (ix1,ix2,ix3) to value \c val
1372  void set(size_t ix1, size_t ix2, size_t ix3, double val) {
1373  size_t sz[3]={ix1,ix2, ix3};
1375  return;
1376  }
1377 
1378  /// Interpolate \c (x,y,z) and return the results
1379  template<class range_t=ub_range, class data_range_t=ubvector_range,
1380  class index_range_t=ubvector_size_t_range>
1381  double interp(double x, double y, double z) {
1382  double arr[3]={x,y,z};
1383  return tensor_grid<vec_t,vec_size_t>::template interpolate
1384  <range_t,data_range_t,index_range_t>(arr);
1385  }
1386 
1387  /// Interpolate \c (x,y,z) and return the results
1388  double interp_linear(double x, double y, double z) {
1389  double arr[3]={x,y,z};
1391  }
1392 
1393  };
1394 
1395  /** \brief Rank 4 tensor with a grid
1396  */
1397  template<class vec_t=std::vector<double>,
1398  class vec_size_t=std::vector<size_t> > class tensor_grid4 :
1399  public tensor_grid<vec_t,vec_size_t> {
1400 
1401  public:
1402 
1403  /// Create an empty tensor
1404  tensor_grid4() : tensor_grid<vec_t,vec_size_t>() {}
1405 
1406  /// Create a rank 4 tensor of size \c (sz,sz2,sz3,sz4)
1407  tensor_grid4(size_t sz, size_t sz2, size_t sz3,
1408  size_t sz4) : tensor_grid<vec_t,vec_size_t>() {
1409  this->rk=4;
1410  this->size.resize(4);
1411  this->size[0]=sz;
1412  this->size[1]=sz2;
1413  this->size[2]=sz3;
1414  this->size[3]=sz4;
1415  size_t tot=sz*sz2*sz3*sz4;
1416  this->data.resize(tot);
1417  this->grid_set=false;
1418  }
1419 
1420 #ifdef O2SCL_NEVER_DEFINED
1421  }{
1422 #endif
1423 
1424  virtual ~tensor_grid4() {
1425  }
1426 
1427  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1428  double &get(size_t ix1, size_t ix2, size_t ix3, size_t ix4) {
1429  size_t sz[4]={ix1,ix2,ix3,ix4};
1431  }
1432 
1433  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1434  const double &get(size_t ix1, size_t ix2, size_t ix3,
1435  size_t ix4) const {
1436  size_t sz[4]={ix1,ix2,ix3,ix4};
1438  }
1439 
1440  /// Set the element indexed by \c (ix1,ix2,ix3,ix4) to value \c val
1441  void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4,
1442  double val) {
1443  size_t sz[4]={ix1,ix2,ix3,ix4};
1445  return;
1446  }
1447 
1448  /// Interpolate \c (x,y,z,a) and return the results
1449  template<class range_t=ub_range, class data_range_t=ubvector_range,
1450  class index_range_t=ubvector_size_t_range>
1451  double interp(double x, double y, double z, double a) {
1452  double arr[4]={x,y,z,a};
1453  return tensor_grid<vec_t,vec_size_t>::template interpolate
1454  <range_t,data_range_t,index_range_t>(arr);
1455  }
1456 
1457  /// Interpolate \c (x,y,z,a) and return the results
1458  double interp_linear(double x, double y, double z, double a) {
1459  double arr[4]={x,y,z,a};
1461  }
1462 
1463  };
1464 
1465 #ifndef DOXYGEN_NO_O2NS
1466 }
1467 #endif
1468 
1469 #endif
1470 
1471 
1472 
bool is_grid_set() const
Return true if the grid has been set.
Definition: tensor_grid.h:320
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:135
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:1055
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:1205
size_t lookup_grid_packed(size_t i, double val)
Lookup index for grid closest to val.
Definition: tensor_grid.h:515
double interpolate(double *vals)
Interpolate values vals into the tensor, returning the result.
Definition: tensor_grid.h:803
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:483
size_t lookup_grid_packed_val(size_t i, double val, double &val2)
Lookup index for grid closest to val.
Definition: tensor_grid.h:540
double get_grid_x(size_t ix)
Get x grid point at index ix.
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:147
void set(const size_vec_t &index, double val)
Set the element indexed by index to value val.
Definition: tensor.h:196
invalid argument supplied by user
Definition: err_hnd.h:59
void unpack_indices(size_t ix, size_vec_t &index)
Unpack the single vector index into indices.
Definition: tensor.h:466
Rank 2 tensor with a grid.
Definition: tensor_grid.h:1266
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:900
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:282
void set_interp_type(size_t interp_type)
Set interpolation type for interpolate()
Definition: tensor_grid.h:769
bool is_slice(std::string name, size_t &ix) const
Return true if slice is already present.
generic failure
Definition: err_hnd.h:61
Rank 3 tensor with a grid.
Definition: tensor_grid.h:1331
size_t lookup_grid(size_t i, double val)
Lookup index for grid closest to val.
Definition: tensor_grid.h:432
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:1340
double interp(double x)
Interpolate x and return the results.
Definition: tensor_grid.h:1250
vec_t & get_data()
Return a reference to the data (for HDF I/O)
Definition: tensor_grid.h:568
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:1381
double get_grid(size_t i, size_t j) const
Lookup jth value on the ith grid.
Definition: tensor_grid.h:397
double interp(double x, double y, double z, double a)
Interpolate (x,y,z,a) and return the results.
Definition: tensor_grid.h:1451
double interp_linear(double x, double y)
Interpolate (x,y) and return the results.
Definition: tensor_grid.h:1321
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:2593
tensor_grid2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
Definition: tensor_grid.h:1275
double interp(double x, double y)
Interpolate (x,y) and return the results.
Definition: tensor_grid.h:1314
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
void copy_slice_align(size_t ix_x, size_t ix_y, size_vec2_t &index, table3d &tab, std::string slice_name)
Create a slice in a table3d object with an aligned grid.
Definition: tensor_grid.h:599
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:1404
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:951
Rank 4 tensor with a grid.
Definition: tensor_grid.h:1398
#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:1109
void set_grid_packed(const vec2_t &grid_vec)
Set the grid.
Definition: tensor_grid.h:344
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:994
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:249
void set_grid(const vec_vec_t &grid_vecs)
Set grid from a vector of vectors of grid points.
Definition: tensor_grid.h:362
void copy_grid(size_t i, rvec_t &v)
Copy grid for index i to vector v.
Definition: tensor_grid.h:386
void set_grid(size_t i, size_t j, double val)
Set the jth value on the ith grid.
Definition: tensor_grid.h:414
double interp_linear(double x, double y, double z, double a)
Interpolate (x,y,z,a) and return the results.
Definition: tensor_grid.h:1458
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor_grid.h:758
double get_val(const vec2_t &gridp)
Get the element closest to grid point gridp.
Definition: tensor_grid.h:223
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:198
double get_grid_y(size_t iy)
Get y grid point at index iy.
void copy_slice_interp_values(size_t ix_x, size_t ix_y, vec2_t &values, table3d &tab, std::string slice_name)
Copy to a slice in a table3d object using interpolation.
Definition: tensor_grid.h:714
size_t total_size() const
Returns the size of the tensor (the product of the sizes over every index)
Definition: tensor.h:427
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.
Searching class for monotonic data with caching.
Definition: search_vec.h:74
size_t itype
Interpolation type.
Definition: tensor_grid.h:128
A data structure containing many slices of two-dimensional data points defined on a grid...
Definition: table3d.h:78
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:100
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:125
double & get(const size_vec_t &index)
Get the element indexed by index.
Definition: tensor.h:237
tensor_grid2()
Create an empty tensor.
Definition: tensor_grid.h:1272
void copy_slice_interp(size_t ix_x, size_t ix_y, size_vec2_t &index, table3d &tab, std::string slice_name)
Copy to a slice in a table3d object using interpolation.
Definition: tensor_grid.h:671
vec_t grid
A rank-sized set of arrays for the grid points.
Definition: tensor_grid.h:122
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:1407
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor.h:184
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:470
void set_val(const vec2_t &grdp, double val)
Set the element closest to grid point grdp to value val.
Definition: tensor_grid.h:169
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:129
tensor_grid3()
Create an empty tensor.
Definition: tensor_grid.h:1337
Linear interpolation (GSL)
Definition: interp.h:207
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:59
tensor_grid1()
Create an empty tensor.
Definition: tensor_grid.h:1211
double interp_linear(double x)
Interpolate x and return the results.
Definition: tensor_grid.h:1256
Linear.
Definition: interp.h:70
tensor_grid1(size_t sz)
Create a rank 2 tensor of size (sz,sz2,sz3)
Definition: tensor_grid.h:1214
Tensor class with arbitrary dimensions.
Definition: tensor.h:120
double interp_linear(double x, double y, double z)
Interpolate (x,y,z) and return the results.
Definition: tensor_grid.h:1388

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