tensor_grid.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, 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, and the tensor interpolation function to copy to
579  the slice named \c slice_name in the table3d object \c tab .
580 
581  If the table3d object does not currently have a grid set, then
582  the grid is automatically set to be the same as that stored in
583  the tensor_grid object associated with ranks \c ix_x and \c
584  iy_y. If the \ref o2scl::table3d object does have a grid set,
585  then the values returned by \ref o2scl::table3d::get_nx() and
586  \ref o2scl::table3d::get_ny() must be equal to the size of the
587  tensor in indices \c ix_x and ix_y, respectively.
588 
589  This currently requires a copy of all the tensor data
590  into the table3d object.
591  */
592  template<class size_vec2_t>
593  void copy_slice_align(size_t ix_x, size_t ix_y, size_vec2_t &index,
594  table3d &tab, std::string slice_name) {
595 
596  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
597  O2SCL_ERR2("Either indices greater than rank or x and y ind",
598  "ices equal in tensor_grid::copy_slice_align().",
599  exc_efailed);
600  }
601 
602  // Get current table3d grid
603  size_t nx, ny;
604  tab.get_size(nx,ny);
605 
606  if (nx==0 && ny==0) {
607 
608  // If there's no grid, just create it
609  std::vector<double> gx, gy;
610  for(size_t i=0;i<this->size[ix_x];i++) {
611  gx.push_back(this->get_grid(ix_x,i));
612  }
613  for(size_t i=0;i<this->size[ix_y];i++) {
614  gy.push_back(this->get_grid(ix_y,i));
615  }
616  nx=gx.size();
617  ny=gy.size();
618  tab.set_xy("x",nx,gx,"y",ny,gy);
619  }
620 
621  // Check that the grids are commensurate
622  if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
623  O2SCL_ERR2("Grids not commensurate in ",
624  "tensor_grid::copy_slice_align().",exc_einval);
625  }
626 
627  // Create slice if not already present
628  size_t is;
629  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
630 
631  // Copy over data
632  for(size_t i=0;i<nx;i++) {
633  for(size_t j=0;j<ny;j++) {
634  index[ix_x]=i;
635  index[ix_y]=j;
636  double val=this->get(index);
637  tab.set(i,j,slice_name,val);
638  }
639  }
640 
641  return;
642  }
643 
644  /** \brief Copy to a slice in a table3d object using interpolation
645 
646  This function uses the grid associated with indices \c ix_x
647  and \c ix_y, and the tensor interpolation function to copy the
648  tensor information to the slice named \c slice_name in the
649  table3d object \c tab .
650 
651  \note This function uses the \ref tensor_grid::interp_linear()
652  for the interpolation.
653  */
654  template<class size_vec2_t>
655  void copy_slice_interp(size_t ix_x, size_t ix_y, size_vec2_t &index,
656  table3d &tab, std::string slice_name) {
657 
658  if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
659  O2SCL_ERR2("Either indices greater than rank or x and y ",
660  "indices equal in tensor_grid::copy_slice_interp().",
661  exc_efailed);
662  }
663 
664  // Get current table3d grid
665  size_t nx, ny;
666  tab.get_size(nx,ny);
667 
668  if (nx==0 && ny==0) {
669  // If there's no grid, then just use the aligned version
670  return copy_slice_align(ix_x,ix_y,index,tab,slice_name);
671  }
672 
673  // Create vector of values to interpolate with
674  std::vector<double> vals(this->rk);
675  for(size_t i=0;i<this->rk;i++) {
676  if (i!=ix_x && i!=ix_y) vals[i]=this->get_grid(i,index[i]);
677  }
678 
679  // Create slice if not already present
680  size_t is;
681  if (!tab.is_slice(slice_name,is)) tab.new_slice(slice_name);
682 
683  // Loop through the table grid to perform the interpolation
684  for(size_t i=0;i<nx;i++) {
685  for(size_t j=0;j<ny;j++) {
686  vals[ix_x]=tab.get_grid_x(i);
687  vals[ix_y]=tab.get_grid_y(j);
688  tab.set(i,j,slice_name,this->interp_linear(vals));
689  }
690  }
691 
692  return;
693  }
694  //@}
695 
696  /// \name Clear method
697  //@{
698  /// Clear the tensor of all data and free allocated memory
699  void clear() {
700  grid.resize(0);
701  grid_set=false;
703  return;
704  }
705  //@}
706 
707  /// \name Interpolation
708  //@{
709  /// Set interpolation type for \ref interpolate()
710  void set_interp_type(size_t interp_type) {
711  itype=interp_type;
712  return;
713  }
714 
715  /** \brief Interpolate values \c vals into the tensor,
716  returning the result
717 
718  This is a quick and dirty implementation of n-dimensional
719  interpolation by recursive application of the 1-dimensional
720  routine from \ref interp_vec, using the base interpolation
721  object specified in the template parameter \c base_interp_t.
722  This will be very slow for sufficiently large data sets.
723 
724  \note This function requires a range objects to obtain ranges
725  of vector objects. In ublas, this is done with
726  <tt>ublas::vector_range</tt> objects, so this function will
727  certainly work for \ref tensor_grid objects built on ublas
728  vector types. There is no corresponding <tt>std::range</tt>,
729  but you may be able to use either <tt>ublas::vector_range</tt>
730  or <tt>Boost.Range</tt> in order to call this function with
731  \ref tensor_grid objects built on <tt>std::vector</tt>.
732  However, this is not fully tested at the moment.
733 
734  \future It should be straightforward to improve the scaling of
735  this algorithm significantly by creating a "window" of local
736  points around the point of interest. This could be done easily
737  by constructing an initial subtensor. However, this should
738  probably be superceded by a more generic alternative which
739  avoids explicit use of the 1-d interpolation types.
740  */
741  template<class range_t=ub_range,
742  class data_range_t=ubvector_range,
743  class index_range_t=ubvector_size_t_range>
744  double interpolate(double *vals) {
745 
746  typedef interp_vec<vec_t> interp_t;
747 
748  if (this->rk==1) {
749 
750  interp_t si(this->size[0],grid,this->data,itype);
751  return si.eval(vals[0]);
752 
753  } else {
754 
755  // Get total number of interpolations at this level
756  size_t ss=1;
757  for(size_t i=1;i<this->rk;i++) ss*=this->size[i];
758 
759  // Create space for y vectors and interpolators
760  std::vector<vec_t> yvec(ss);
761  std::vector<interp_t *> si(ss);
762  for(size_t i=0;i<ss;i++) {
763  yvec[i].resize(this->size[0]);
764  }
765 
766  // Create space for interpolation results
767  tensor_grid tdat;
768  index_range_t size_new(this->size,ub_range(1,this->rk));
769  tdat.resize(this->rk-1,size_new);
770 
771  // Set grid for temporary tensor
772  data_range_t grid_new(grid,ub_range(this->size[0],grid.size()));
773  tdat.set_grid_packed(grid_new);
774 
775  // Create starting coordinate and counter
776  vec_size_t co(this->rk);
777  for(size_t i=0;i<this->rk;i++) co[i]=0;
778  size_t cnt=0;
779 
780  // Loop over every interpolation
781  bool done=false;
782  while(done==false) {
783 
784  // Fill yvector with the appropriate data
785  for(size_t i=0;i<this->size[0];i++) {
786  co[0]=i;
787  yvec[cnt][i]=this->get(co);
788  }
789 
790  si[cnt]=new interp_t(this->size[0],grid,yvec[cnt],itype);
791 
792  index_range_t co2(co,ub_range(1,this->rk));
793  tdat.set(co2,si[cnt]->eval(vals[0]));
794 
795  // Go to next interpolation
796  cnt++;
797  co[this->rk-1]++;
798  // carry if necessary
799  for(int j=((int)this->rk)-1;j>0;j--) {
800  if (co[j]>=this->size[j]) {
801  co[j]=0;
802  co[j-1]++;
803  }
804  }
805 
806  // Test if done
807  if (cnt==ss) done=true;
808 
809  // End of while loop
810  }
811 
812  // Now call the next level of interpolation
813  double res=tdat.interpolate(vals+1);
814 
815  for(size_t i=0;i<ss;i++) {
816  delete si[i];
817  }
818 
819  return res;
820  }
821  }
822 
823  /** \brief Perform a linear interpolation of \c v into the
824  function implied by the tensor and grid
825 
826  This performs multi-dimensional linear interpolation (or
827  extrapolation) It works by first using \ref o2scl::search_vec
828  to find the interval containing (or closest to) the specified
829  point in each direction and constructing the corresponding
830  hypercube of size \f$ 2^{\mathrm{rank}} \f$ containing \c v.
831  It then calls \ref interp_linear_power_two() to perform the
832  interpolation in that hypercube.
833 
834  \future This starts with a small copy, which can be eliminated
835  by creating a new version of interp_linear_power_two
836  which accepts an offset vector parameter so that the
837  first interpolation is offset. Remaining interpolations
838  don't need to be offset because the tensor has to be
839  created from the previous interpolation round.
840  */
841  template<class vec2_size_t> double interp_linear(vec2_size_t &v) {
842 
843  // Find the the corner of the hypercube containing v
844  size_t rgs=0;
845  std::vector<size_t> loc(this->rk);
846  std::vector<double> gnew;
847  for(size_t i=0;i<this->rk;i++) {
848  std::vector<double> grid_unpacked(this->size[i]);
849  for(size_t j=0;j<this->size[i];j++) {
850  grid_unpacked[j]=grid[j+rgs];
851  }
852  search_vec<std::vector<double> > sv(this->size[i],grid_unpacked);
853  loc[i]=sv.find(v[i]);
854  gnew.push_back(grid_unpacked[loc[i]]);
855  gnew.push_back(grid_unpacked[loc[i]+1]);
856  rgs+=this->size[i];
857  }
858 
859  // Now construct a 2^{rk}-sized tensor containing only that
860  // hypercube
861  std::vector<size_t> snew(this->rk);
862  for(size_t i=0;i<this->rk;i++) {
863  snew[i]=2;
864  }
865  tensor_grid tnew(this->rk,snew);
866  tnew.set_grid_packed(gnew);
867 
868  // Copy over the relevant data
869  for(size_t i=0;i<tnew.total_size();i++) {
870  std::vector<size_t> index_new(this->rk), index_old(this->rk);
871  tnew.unpack_indices(i,index_new);
872  for(size_t j=0;j<this->rk;j++) index_old[j]=index_new[j]+loc[j];
873  tnew.set(index_new,this->get(index_old));
874  }
875 
876  // Now use interp_power_two()
877  return tnew.interp_linear_power_two(v);
878  }
879 
880  /** \brief Perform linear interpolation assuming that all
881  indices can take only two values
882 
883  This function works by recursively slicing the hypercube of
884  size \f$ 2^{\mathrm{rank}} \f$ into a hypercube of size \f$
885  2^{\mathrm{rank-1}} \f$ performing linear interpolation for
886  each pair of points.
887 
888  \note This is principally a function for internal use
889  by \ref interp_linear().
890  */
891  template<class vec2_size_t>
892  double interp_linear_power_two(vec2_size_t &v) {
893 
894  if (this->rk==1) {
895  return this->data[0]+(this->data[1]-this->data[0])/
896  (grid[1]-grid[0])*(v[0]-grid[0]);
897  }
898 
899  size_t last=this->rk-1;
900  double frac=(v[last]-get_grid(last,0))/
901  (get_grid(last,1)-get_grid(last,0));
902 
903  // Create new size vector and grid
904  tensor_grid tnew(this->rk-1,this->size);
905  tnew.set_grid_packed(grid);
906 
907  // Create data in new tensor, removing the last index through
908  // linear interpolation
909  for(size_t i=0;i<tnew.total_size();i++) {
910  std::vector<size_t> index(this->rk);
911  tnew.unpack_indices(i,index);
912  index[this->rk-1]=0;
913  double val_lo=this->get(index);
914  index[this->rk-1]=1;
915  double val_hi=this->get(index);
916  tnew.set(index,val_lo+frac*(val_hi-val_lo));
917  }
918 
919  // Recursively interpolate the smaller tensor
920  return tnew.interp_linear_power_two(v);
921  }
922 
923  /** \brief Perform a linear interpolation of <tt>v[1]</tt>
924  to <tt>v[n-1]</tt> resulting in a vector
925 
926  \note The type <tt>vec2_t</tt> for the vector <tt>res</tt>
927  must have a <tt>resize()</tt> method.
928 
929  This performs multi-dimensional linear interpolation (or
930  extrapolation) in the last <tt>n-1</tt> indices of the
931  rank-<tt>n</tt> tensor leaving the first index free and places
932  the results in the vector \c res.
933  */
934  template<class vec2_size_t, class vec2_t>
935  void interp_linear_vec0(vec2_size_t &v, vec2_t &res) {
936 
937  // Find the the corner of the hypercube containing v
938  size_t rgs=0;
939  std::vector<size_t> loc(this->rk);
940  std::vector<double> gnew;
941  for(size_t i=0;i<this->size[0];i++) {
942  gnew.push_back(grid[i]);
943  }
944  rgs=this->size[0];
945  loc[0]=0;
946  for(size_t i=1;i<this->rk;i++) {
947  std::vector<double> grid_unpacked(this->size[i]);
948  for(size_t j=0;j<this->size[i];j++) {
949  grid_unpacked[j]=grid[j+rgs];
950  }
951  search_vec<std::vector<double> > sv(this->size[i],grid_unpacked);
952  loc[i]=sv.find(v[i]);
953  gnew.push_back(grid_unpacked[loc[i]]);
954  gnew.push_back(grid_unpacked[loc[i]+1]);
955  rgs+=this->size[i];
956  }
957 
958  // Now construct a n*2^{rk-1}-sized tensor containing only that
959  // hypercube
960  std::vector<size_t> snew(this->rk);
961  snew[0]=this->size[0];
962  for(size_t i=1;i<this->rk;i++) {
963  snew[i]=2;
964  }
965  tensor_grid tnew(this->rk,snew);
966  tnew.set_grid_packed(gnew);
967 
968  // Copy over the relevant data
969  for(size_t i=0;i<tnew.total_size();i++) {
970  std::vector<size_t> index_new(this->rk), index_old(this->rk);
971  tnew.unpack_indices(i,index_new);
972  for(size_t j=0;j<this->rk;j++) {
973  index_old[j]=index_new[j]+loc[j];
974  }
975  tnew.set(index_new,this->get(index_old));
976  }
977 
978  // Now use interp_power_two_vec0()
979  tnew.interp_linear_power_two_vec0(v,res);
980 
981  return;
982  }
983 
984  /** \brief Perform linear interpolation assuming that the last
985  <tt>n-1</tt> indices can take only two values
986 
987  This function performs linear interpolation assuming that the
988  last <tt>n-1</tt> indices can take only two values and placing
989  the result into <tt>res</tt>.
990 
991  \note The type <tt>vec2_t</tt> for the vector <tt>res</tt>
992  must have a <tt>resize()</tt> method. This is principally a
993  function for internal use by \ref interp_linear_vec0().
994  */
995  template<class vec2_size_t, class vec2_t>
996  void interp_linear_power_two_vec0(vec2_size_t &v, vec2_t &res) {
997 
998  if (this->rk==2) {
999  size_t n=this->size[0];
1000  res.resize(n);
1001  vec_size_t ix0(2), ix1(2);
1002  ix0[1]=0;
1003  ix1[1]=1;
1004  for(size_t i=0;i<n;i++) {
1005  ix0[0]=i;
1006  ix1[0]=i;
1007  res[i]=this->get(ix0)+(this->get(ix1)-this->get(ix0))/
1008  (grid[n+1]-grid[n])*(v[1]-grid[n]);
1009  }
1010  return;
1011  }
1012 
1013  size_t last=this->rk-1;
1014  double frac=(v[last]-get_grid(last,0))/
1015  (get_grid(last,1)-get_grid(last,0));
1016 
1017  // Create new size vector and grid
1018  tensor_grid tnew(this->rk-1,this->size);
1019  tnew.set_grid_packed(grid);
1020 
1021  // Create data in new tensor, removing the last index through
1022  // linear interpolation
1023  for(size_t i=0;i<tnew.total_size();i++) {
1024  std::vector<size_t> index(this->rk);
1025  tnew.unpack_indices(i,index);
1026  index[this->rk-1]=0;
1027  double val_lo=this->get(index);
1028  index[this->rk-1]=1;
1029  double val_hi=this->get(index);
1030  tnew.set(index,val_lo+frac*(val_hi-val_lo));
1031  }
1032 
1033  // Recursively interpolate the smaller tensor
1034  tnew.interp_linear_power_two_vec0(v,res);
1035 
1036  return;
1037  }
1038 
1039  /** \brief Perform a linear interpolation of <tt>v</tt>
1040  into the tensor leaving one index free resulting in a vector
1041 
1042  This performs multi-dimensional linear interpolation (or
1043  extrapolation) in the last <tt>n-1</tt> indices of the
1044  rank-<tt>n</tt> tensor leaving the first index free and places
1045  the results in the vector \c res.
1046 
1047  \future This function could be more efficient.
1048  */
1049  template<class vec2_size_t, class vec2_t>
1050  void interp_linear_vec(vec2_size_t &v, size_t ifree, vec2_t &res) {
1051 
1052  size_t n=this->size[ifree];
1053 
1054  // This function uses interp_linear_power_two_vec0(), so it
1055  // works by remapping the indices. This defines the remapping.
1056  std::vector<size_t> map;
1057  map.push_back(ifree);
1058  for(size_t i=0;i<this->rk;i++) {
1059  if (i!=ifree) {
1060  map.push_back(i);
1061  }
1062  }
1063 
1064  // Find the the corner of the hypercube containing v
1065  size_t rgs=0;
1066  vec_size_t loc(this->rk);
1067  loc[ifree]=0;
1068  for(size_t i=0;i<this->rk;i++) {
1069  vec_t grid_unpacked(this->size[i]);
1070  for(size_t j=0;j<this->size[i];j++) {
1071  grid_unpacked[j]=grid[j+rgs];
1072  }
1073  search_vec<vec_t> sv(this->size[i],grid_unpacked);
1074  if (i!=ifree) {
1075  loc[i]=sv.find(v[i]);
1076  }
1077  rgs+=this->size[i];
1078  }
1079 
1080  // Compute the remapped grid and interpolating vector
1081  std::vector<double> gnew, vnew;
1082  for(size_t new_ix=0;new_ix<this->rk;new_ix++) {
1083  for(size_t old_ix=0;old_ix<this->rk;old_ix++) {
1084  if (map[new_ix]==old_ix) {
1085  vnew.push_back(v[old_ix]);
1086  if (old_ix==ifree) {
1087  for(size_t j=0;j<this->size[old_ix];j++) {
1088  gnew.push_back(this->get_grid(old_ix,j));
1089  }
1090  } else {
1091  gnew.push_back(this->get_grid(old_ix,loc[old_ix]));
1092  gnew.push_back(this->get_grid(old_ix,loc[old_ix]+1));
1093  }
1094  }
1095  }
1096  }
1097 
1098  // Now construct a n*2^{rk-1}-sized tensor containing only the
1099  // hypercube needed to do the interpolation
1100 
1101  // Specify the size of each rank
1102  std::vector<size_t> snew;
1103  snew.push_back(n);
1104  for(size_t i=0;i<this->rk;i++) {
1105  if (i!=ifree) {
1106  snew.push_back(2);
1107  }
1108  }
1109 
1110  // Create the tensor and set the grid
1111  tensor_grid tnew(this->rk,snew);
1112  tnew.set_grid_packed(gnew);
1113 
1114  // Copy over the relevant data
1115  for(size_t i=0;i<tnew.total_size();i++) {
1116  std::vector<size_t> index_new(this->rk), index_old(this->rk);
1117  tnew.unpack_indices(i,index_new);
1118  for(size_t j=0;j<this->rk;j++) {
1119  index_old[map[j]]=index_new[j]+loc[map[j]];
1120  }
1121  tnew.set(index_new,this->get(index_old));
1122  }
1123 
1124  // Now use interp_power_two_vec()
1125  tnew.interp_linear_power_two_vec0(vnew,res);
1126 
1127  return;
1128  }
1129  //@}
1130 
1131  template<class vecf_t, class vecf_size_t> friend void o2scl_hdf::hdf_output
1133  std::string name);
1134 
1135  template<class vecf_t, class vecf_size_t> friend void o2scl_hdf::hdf_input
1137  std::string name);
1138 
1139  };
1140 
1141  /** \brief Rank 1 tensor with a grid
1142 
1143  \future Make rank-specific get_val and set_val functions?
1144  */
1145  template<class vec_t=std::vector<double>,
1146  class vec_size_t=std::vector<size_t> > class tensor_grid1 :
1147  public tensor_grid<vec_t,vec_size_t> {
1148 
1149  public:
1150 
1151  /// Create an empty tensor
1152  tensor_grid1() : tensor_grid<vec_t,vec_size_t>() {}
1153 
1154  /// Create a rank 2 tensor of size \c (sz,sz2,sz3)
1155  tensor_grid1(size_t sz) : tensor_grid<vec_t,vec_size_t>() {
1156  this->rk=1;
1157  this->size.resize(1);
1158  this->size[0]=sz;
1159  this->data.resize(sz);
1160  this->grid_set=false;
1161  }
1162 
1163 #ifdef O2SCL_NEVER_DEFINED
1164  }{
1165 #endif
1166 
1167  virtual ~tensor_grid1() {
1168  }
1169 
1170  /// Get the element indexed by \c (ix1)
1171  double &get(size_t ix1) {
1172  size_t sz[1]={ix1};
1174  }
1175 
1176  /// Get the element indexed by \c (ix1)
1177  const double &get(size_t ix1) const {
1178  size_t sz[1]={ix1};
1180  }
1181 
1182  /// Set the element indexed by \c (ix1) to value \c val
1183  void set(size_t ix1, double val) {
1184  size_t sz[1]={ix1};
1186  }
1187 
1188  /// Interpolate \c x and return the results
1189  template<class range_t=ub_range, class data_range_t=ubvector_range,
1190  class index_range_t=ubvector_size_t_range>
1191  double interp(double x) {
1192  return tensor_grid<vec_t,vec_size_t>::template interpolate
1193  <range_t,data_range_t,index_range_t>(&x);
1194  }
1195 
1196  /// Interpolate \c x and return the results
1197  double interp_linear(double x) {
1198  double arr[1]={x};
1200  }
1201 
1202  };
1203 
1204  /** \brief Rank 2 tensor with a grid
1205  */
1206  template<class vec_t=std::vector<double>,
1207  class vec_size_t=std::vector<size_t> > class tensor_grid2 :
1208  public tensor_grid<vec_t,vec_size_t> {
1209 
1210  public:
1211 
1212  /// Create an empty tensor
1213  tensor_grid2() : tensor_grid<vec_t,vec_size_t>() {}
1214 
1215  /// Create a rank 2 tensor of size \c (sz,sz2)
1216  tensor_grid2(size_t sz, size_t sz2) : tensor_grid<vec_t,vec_size_t>() {
1217  this->rk=2;
1218  this->size.resize(2);
1219  this->size[0]=sz;
1220  this->size[1]=sz2;
1221  size_t tot=sz*sz2;
1222  this->data.resize(tot);
1223  this->grid_set=false;
1224  }
1225 
1226 #ifdef O2SCL_NEVER_DEFINED
1227  }{
1228 #endif
1229 
1230  virtual ~tensor_grid2() {
1231  }
1232 
1233  /// Get the element indexed by \c (ix1,ix2)
1234  double &get(size_t ix1, size_t ix2) {
1235  size_t sz[2]={ix1,ix2};
1237  }
1238 
1239  /// Get the element indexed by \c (ix1,ix2)
1240  const double &get(size_t ix1, size_t ix2) const {
1241  size_t sz[2]={ix1,ix2};
1243  }
1244 
1245  /// Set the element indexed by \c (ix1,ix2) to value \c val
1246  void set(size_t ix1, size_t ix2, double val) {
1247  size_t sz[2]={ix1,ix2};
1249  return;
1250  }
1251 
1252  /// Interpolate \c (x,y) and return the results
1253  template<class range_t=ub_range, class data_range_t=ubvector_range,
1254  class index_range_t=ubvector_size_t_range>
1255  double interp(double x, double y) {
1256  double arr[2]={x,y};
1257  return tensor_grid<vec_t,vec_size_t>::template interpolate
1258  <range_t,data_range_t,index_range_t>(arr);
1259  }
1260 
1261  /// Interpolate \c (x,y) and return the results
1262  double interp_linear(double x, double y) {
1263  double arr[2]={x,y};
1265  }
1266 
1267  };
1268 
1269  /** \brief Rank 3 tensor with a grid
1270  */
1271  template<class vec_t=std::vector<double>,
1272  class vec_size_t=std::vector<size_t> > class tensor_grid3 :
1273  public tensor_grid<vec_t,vec_size_t> {
1274 
1275  public:
1276 
1277  /// Create an empty tensor
1278  tensor_grid3() : tensor_grid<vec_t,vec_size_t>() {}
1279 
1280  /// Create a rank 3 tensor of size \c (sz,sz2,sz3)
1281  tensor_grid3(size_t sz, size_t sz2, size_t sz3) :
1282  tensor_grid<vec_t,vec_size_t>() {
1283  this->rk=3;
1284  this->size.resize(3);
1285  this->size[0]=sz;
1286  this->size[1]=sz2;
1287  this->size[2]=sz3;
1288  size_t tot=sz*sz2*sz3;
1289  this->data.resize(tot);
1290  this->grid_set=false;
1291  }
1292 
1293 #ifdef O2SCL_NEVER_DEFINED
1294  }{
1295 #endif
1296 
1297  virtual ~tensor_grid3() {
1298  }
1299 
1300  /// Get the element indexed by \c (ix1,ix2,ix3)
1301  double &get(size_t ix1, size_t ix2, size_t ix3) {
1302  size_t sz[3]={ix1,ix2,ix3};
1304  }
1305 
1306  /// Get the element indexed by \c (ix1,ix2,ix3)
1307  const double &get(size_t ix1, size_t ix2, size_t ix3) const {
1308  size_t sz[3]={ix1,ix2,ix3};
1310  }
1311 
1312  /// Set the element indexed by \c (ix1,ix2,ix3) to value \c val
1313  void set(size_t ix1, size_t ix2, size_t ix3, double val) {
1314  size_t sz[3]={ix1,ix2, ix3};
1316  return;
1317  }
1318 
1319  /// Interpolate \c (x,y,z) and return the results
1320  template<class range_t=ub_range, class data_range_t=ubvector_range,
1321  class index_range_t=ubvector_size_t_range>
1322  double interp(double x, double y, double z) {
1323  double arr[3]={x,y,z};
1324  return tensor_grid<vec_t,vec_size_t>::template interpolate
1325  <range_t,data_range_t,index_range_t>(arr);
1326  }
1327 
1328  /// Interpolate \c (x,y,z) and return the results
1329  double interp_linear(double x, double y, double z) {
1330  double arr[3]={x,y,z};
1332  }
1333 
1334  };
1335 
1336  /** \brief Rank 4 tensor with a grid
1337  */
1338  template<class vec_t=std::vector<double>,
1339  class vec_size_t=std::vector<size_t> > class tensor_grid4 :
1340  public tensor_grid<vec_t,vec_size_t> {
1341 
1342  public:
1343 
1344  /// Create an empty tensor
1345  tensor_grid4() : tensor_grid<vec_t,vec_size_t>() {}
1346 
1347  /// Create a rank 4 tensor of size \c (sz,sz2,sz3,sz4)
1348  tensor_grid4(size_t sz, size_t sz2, size_t sz3,
1349  size_t sz4) : tensor_grid<vec_t,vec_size_t>() {
1350  this->rk=4;
1351  this->size.resize(4);
1352  this->size[0]=sz;
1353  this->size[1]=sz2;
1354  this->size[2]=sz3;
1355  this->size[3]=sz4;
1356  size_t tot=sz*sz2*sz3*sz4;
1357  this->data.resize(tot);
1358  this->grid_set=false;
1359  }
1360 
1361 #ifdef O2SCL_NEVER_DEFINED
1362  }{
1363 #endif
1364 
1365  virtual ~tensor_grid4() {
1366  }
1367 
1368  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1369  double &get(size_t ix1, size_t ix2, size_t ix3, size_t ix4) {
1370  size_t sz[4]={ix1,ix2,ix3,ix4};
1372  }
1373 
1374  /// Get the element indexed by \c (ix1,ix2,ix3,ix4)
1375  const double &get(size_t ix1, size_t ix2, size_t ix3,
1376  size_t ix4) const {
1377  size_t sz[4]={ix1,ix2,ix3,ix4};
1379  }
1380 
1381  /// Set the element indexed by \c (ix1,ix2,ix3,ix4) to value \c val
1382  void set(size_t ix1, size_t ix2, size_t ix3, size_t ix4,
1383  double val) {
1384  size_t sz[4]={ix1,ix2,ix3,ix4};
1386  return;
1387  }
1388 
1389  /// Interpolate \c (x,y,z,a) and return the results
1390  template<class range_t=ub_range, class data_range_t=ubvector_range,
1391  class index_range_t=ubvector_size_t_range>
1392  double interp(double x, double y, double z, double a) {
1393  double arr[4]={x,y,z,a};
1394  return tensor_grid<vec_t,vec_size_t>::template interpolate
1395  <range_t,data_range_t,index_range_t>(arr);
1396  }
1397 
1398  /// Interpolate \c (x,y,z,a) and return the results
1399  double interp_linear(double x, double y, double z, double a) {
1400  double arr[4]={x,y,z,a};
1402  }
1403 
1404  };
1405 
1406 #ifndef DOXYGEN_NO_O2NS
1407 }
1408 #endif
1409 
1410 #endif
1411 
1412 
1413 
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:996
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:1146
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:744
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.
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:1207
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:841
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:710
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:1272
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:1281
double interp(double x)
Interpolate x and return the results.
Definition: tensor_grid.h:1191
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:1322
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:1392
size_t find(const double x0) const
Search an increasing or decreasing vector for the interval containing x0
Definition: search_vec.h:137
double interp_linear(double x, double y)
Interpolate (x,y) and return the results.
Definition: tensor_grid.h:1262
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:2384
tensor_grid2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
Definition: tensor_grid.h:1216
double interp(double x, double y)
Interpolate (x,y) and return the results.
Definition: tensor_grid.h:1255
#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:593
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:1345
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:892
Rank 4 tensor with a grid.
Definition: tensor_grid.h:1339
#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:1050
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:935
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:1399
void clear()
Clear the tensor of all data and free allocated memory.
Definition: tensor_grid.h:699
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.
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:79
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:96
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:1213
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:655
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:1348
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:126
tensor_grid3()
Create an empty tensor.
Definition: tensor_grid.h:1278
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:1152
double interp_linear(double x)
Interpolate x and return the results.
Definition: tensor_grid.h:1197
Linear.
Definition: interp.h:70
tensor_grid1(size_t sz)
Create a rank 2 tensor of size (sz,sz2,sz3)
Definition: tensor_grid.h:1155
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:1329

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