23 #ifndef O2SCL_TENSOR_H 24 #define O2SCL_TENSOR_H 35 #include <boost/numeric/ublas/vector.hpp> 36 #include <boost/numeric/ublas/vector_proxy.hpp> 37 #include <boost/numeric/ublas/matrix.hpp> 38 #include <boost/numeric/ublas/matrix_proxy.hpp> 40 #include <gsl/gsl_matrix.h> 41 #include <gsl/gsl_ieee_utils.h> 43 #include <o2scl/err_hnd.h> 44 #include <o2scl/interp.h> 46 #ifndef DOXYGEN_NO_O2NS 119 template<
class vec_t=std::vector<
double>,
120 class vec_
size_t=std::vector<
size_t> >
class tensor {
124 #ifndef DOXYGEN_INTERNAL 153 template<
class size_vec_t>
154 tensor(
size_t rank,
const size_vec_t &dim) {
159 for(
size_t i=0;i<
rk;i++) {
162 O2SCL_ERR((((std::string)
"Requested zero size with non-zero ")+
163 "rank for index "+
szttos(i)+
164 " in tensor::tensor(size_t,size_vec_t)").c_str(),
170 for(
size_t i=0;i<
rk;i++) {
195 template<
class size_vec_t>
196 void set(
const size_vec_t &index,
double val) {
197 #if O2SCL_NO_RANGE_CHECK 202 if (index[0]>=size[0]) {
204 " greater than or equal to size[0]="+
205 szttos(size[0])+
" in tensor::set().").c_str(),
210 for(
size_t i=1;i<
rk;i++) {
211 #if O2SCL_NO_RANGE_CHECK 213 if (index[i]>=size[i]) {
215 szttos(index[i])+
" greater than or equal to size "+
216 szttos(size[i])+
" in tensor::set().").c_str(),
237 template<
class size_vec_t>
double &
get(
const size_vec_t &index) {
238 #if O2SCL_NO_RANGE_CHECK 243 if (index[0]>=size[0]) {
245 " greater than or equal to size[0]="+
246 szttos(size[0])+
" in tensor::get().").c_str(),
251 for(
size_t i=1;i<
rk;i++) {
252 #if O2SCL_NO_RANGE_CHECK 254 if (index[i]>=size[i]) {
256 szttos(index[i])+
" greater than or equal to size "+
257 szttos(size[i])+
" in tensor::get().").c_str(),
268 template<
class size_vec_t>
269 double const &
get(
const size_vec_t &index)
const {
270 #if O2SCL_NO_RANGE_CHECK 275 if (index[0]>=size[0]) {
277 " greater than or equal to size[0]="+
278 szttos(size[0])+
" in tensor::get() (const).").c_str(),
283 for(
size_t i=1;i<
rk;i++) {
284 #if O2SCL_NO_RANGE_CHECK 286 if (index[i]>=size[i]) {
288 szttos(index[i])+
" greater than or equal to size "+
289 szttos(size[i])+
" in tensor::get() (const).").c_str(),
300 typedef boost::numeric::ublas::vector_slice<
302 typedef boost::numeric::ublas::slice slice;
328 template<
class size_vec_t>
332 " greater than or equal to rank "+
szttos(rk)+
333 " in tensor::vector_slice()").c_str(),
339 for(
size_t i=1;i<
rk;i++) {
341 if (i!=ix) start+=index[i];
344 for(
size_t i=ix+1;i<
rk;i++) stride*=size[i];
345 return ubvector_slice(data,slice(start,stride,size[ix]));
349 #ifdef O2SCL_NEVER_DEFINED 357 template<
class size_vec_t> ubmatrix_array matrix_slice() {
373 template<
class size_vec_t>
374 void resize(
size_t rank,
const size_vec_t &dim) {
375 for(
size_t i=0;i<rank;i++) {
377 O2SCL_ERR((((std::string)
"Requested zero size with non-zero ")+
378 "rank for index "+
szttos(i)+
" in tensor::"+
388 for(
size_t i=0;i<
rk;i++) {
407 " greater than or equal to rank "+
szttos(rk)+
408 " in tensor::get_size()").c_str(),
430 for(
size_t i=0;i<
rk;i++) tot*=size[i];
438 template<
class size_vec_t>
444 if (index[0]>=size[0]) {
446 " greater than or equal to size[0]="+
447 szttos(size[0])+
" in tensor::pack_indices().").c_str(),
451 for(
size_t i=1;i<
rk;i++) {
452 if (index[i]>=size[i]) {
454 szttos(index[i])+
" greater than or equal to size "+
455 szttos(size[i])+
" in tensor::pack_indices().").c_str(),
465 template<
class size_vec_t>
469 " greater than or equal to total size"+
471 " in tensor::unpack_indices().").c_str(),
475 size_t ix2, sub_size;
476 for(
size_t i=0;i<
rk;i++) {
481 for(
size_t j=i+1;j<
rk;j++) sub_size*=size[j];
482 index[i]=ix/sub_size;
484 ix-=sub_size*(ix/sub_size);
502 size_t ix=o2scl::vector_min_index<vec_t,double>(
total_size(),
data);
510 void min(vec_size_t &index,
double &val) {
526 size_t ix=o2scl::vector_max_index<vec_t,double>(
total_size(),
data);
534 void max(vec_size_t &index,
double &val) {
550 size_t ix_min, ix_max;
551 o2scl::vector_minmax_index<vec_t,double>
561 void minmax(vec_size_t &index,
size_t &index_min,
double &
min,
562 size_t &index_max,
double &
max) {
563 size_t ix_min, ix_max;
576 template<
class vec_t=std::vector<
double>,
577 class vec_
size_t=std::vector<
size_t> >
class tensor1 :
578 public tensor<vec_t,vec_size_t> {
593 double &
get(
size_t ix) {
598 const double &
get(
size_t ix)
const {
603 void set(
size_t index,
double val)
611 template<
class size_vec_t>
612 void set(
const size_vec_t &index,
double val) {
631 template<
class vec_t=std::vector<
double>,
632 class vec_
size_t=std::vector<
size_t> >
class tensor2 :
633 public tensor<vec_t,vec_size_t> {
643 this->
size.resize(2);
647 this->
data.resize(tot);
651 double &
get(
size_t ix1,
size_t ix2) {
652 size_t sz[2]={ix1,ix2};
657 const double &
get(
size_t ix1,
size_t ix2)
const {
658 size_t sz[2]={ix1,ix2};
663 void set(
size_t ix1,
size_t ix2,
double val) {
664 size_t sz[2]={ix1,ix2};
674 template<
class size_vec_t>
675 void set(
const size_vec_t &index,
double val) {
682 {
return this->
data[ix*this->
size[1]+iy]; }
686 {
return this->
data[ix*this->
size[1]+iy]; }
691 template<
class vec_t=std::vector<
double>,
692 class vec_
size_t=std::vector<
size_t> >
class tensor3 :
693 public tensor<vec_t,vec_size_t> {
702 tensor<vec_t,vec_size_t>() {
704 this->
size.resize(3);
708 size_t tot=sz*sz2*sz3;
709 this->
data.resize(tot);
713 double &
get(
size_t ix1,
size_t ix2,
size_t ix3) {
714 size_t sz[3]={ix1,ix2,ix3};
719 const double &
get(
size_t ix1,
size_t ix2,
size_t ix3)
const {
720 size_t sz[3]={ix1,ix2,ix3};
725 void set(
size_t ix1,
size_t ix2,
size_t ix3,
double val) {
726 size_t sz[3]={ix1,ix2,ix3};
736 template<
class size_vec_t>
737 void set(
const size_vec_t &index,
double val) {
746 template<
class vec_t=std::vector<
double>,
747 class vec_
size_t=std::vector<
size_t> >
class tensor4 :
748 public tensor<vec_t,vec_size_t> {
756 tensor4(
size_t sz,
size_t sz2,
size_t sz3,
size_t sz4) :
757 tensor<vec_t,vec_size_t>() {
759 this->
size.resize(4);
764 size_t tot=sz*sz2*sz3*sz4;
765 this->
data.resize(tot);
769 double &
get(
size_t ix1,
size_t ix2,
size_t ix3,
size_t ix4) {
770 size_t sz[4]={ix1,ix2,ix3,ix4};
775 const double &
get(
size_t ix1,
size_t ix2,
size_t ix3,
777 size_t sz[4]={ix1,ix2,ix3,ix4};
782 void set(
size_t ix1,
size_t ix2,
size_t ix3,
size_t ix4,
784 size_t sz[4]={ix1,ix2,ix3,ix4};
794 template<
class size_vec_t>
795 void set(
const size_vec_t &index,
double val) {
801 #ifndef DOXYGEN_NO_O2NS tensor(size_t rank, const size_vec_t &dim)
Create a tensor of rank rank with sizes given in dim.
size_t pack_indices(const size_vec_t &index)
Pack the indices into a single vector index.
double max_value()
Compute the maximum value in the tensor.
const double & operator()(size_t ix, size_t iy) const
Get the element indexed by (ix1,ix2) (const version)
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
void max(vec_size_t &index, double &val)
Compute the index and value of the maximum value in the tensor and return the maximum.
void min(vec_size_t &index, double &val)
Compute the index of the minimum value in the tensor and return the minimum.
const double & operator[](size_t ix) const
Get an element using array-like indexing (const version)
void set(const size_vec_t &index, double val)
Set the element indexed by index to value val.
invalid argument supplied by user
double & operator()(size_t ix, size_t iy)
Get the element indexed by (ix1,ix2)
void unpack_indices(size_t ix, size_vec_t &index)
Unpack the single vector index into indices.
void max_index(vec_size_t &index)
Compute the index of the maximum value in the tensor.
void resize(size_t rank, const size_vec_t &dim)
Resize the tensor to rank rank with sizes given in dim.
void minmax(vec_size_t &index, size_t &index_min, double &min, size_t &index_max, double &max)
Compute the indices and values of the maximum and minimum in the tensor.
size_t get_size(size_t i) const
Returns the size of the ith index.
void min_index(vec_size_t &index)
Compute the index of the minimum value in the tensor.
size_t get_rank() const
Return the rank of the tensor.
tensor2()
Create an empty tensor.
vec_size_t size
A rank-sized vector of the sizes of each dimension.
const vec_t & get_data() const
Return the full data vector.
tensor3()
Create an empty tensor.
void minmax_value(double &min, double &max)
Compute the minimum and maximum values in the tensor.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
void set_all(double x)
Set all elements in a tensor to some fixed value.
tensor1()
Create an empty tensor.
void minmax_index(vec_size_t &index_min, vec_size_t &index_max)
Compute the indices of the minimum and maximum values in the tensor.
tensor1(size_t sz)
Create a rank 1 tensory of size sz.
tensor4()
Create an empty tensor.
ubvector_slice vector_slice(size_t ix, const size_vec_t &index)
Fix all but one index to create a vector.
size_t total_size() const
Returns the size of the tensor (the product of the sizes over every index)
tensor()
Create an empty tensor with zero rank.
const vec_size_t & get_size_arr() const
Return the full vector of sizes.
double & get(const size_vec_t &index)
Get the element indexed by index.
void clear()
Clear the tensor of all data and free allocated memory.
Invalid index for array or matrix.
tensor3(size_t sz, size_t sz2, size_t sz3)
Create a rank 3 tensor of size (sz,sz2,sz3)
double & operator[](size_t ix)
Get an element using array-like indexing.
std::string szttos(size_t x)
Convert a size_t to a string.
tensor2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
const double & operator()(size_t ix) const
Get an element using operator() (const version)
tensor4(size_t sz, size_t sz2, size_t sz3, size_t sz4)
Create a rank 4 tensor of size (sz,sz2,sz3,sz4)
double & operator()(size_t ix)
Get an element using operator()
Tensor class with arbitrary dimensions.
double min_value()
Compute the minimum value in the tensor.