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> 45 #include <o2scl/table3d.h> 46 #include <o2scl/misc.h> 48 #ifndef DOXYGEN_NO_O2NS 111 template<
class data_t=
double,
class vec_t=std::vector<data_t>,
112 class vec_
size_t=std::vector<
size_t> >
class tensor {
116 #ifndef DOXYGEN_INTERNAL 145 template<
class size_vec_t>
146 tensor(
size_t rank,
const size_vec_t &dim) {
151 for(
size_t i=0;i<
rk;i++) {
154 O2SCL_ERR((((std::string)
"Requested zero size with non-zero ")+
155 "rank for index "+
szttos(i)+
156 " in tensor::tensor(size_t,size_vec_t)").c_str(),
162 for(
size_t i=0;i<
rk;i++) {
179 if (data.size()!=0) {
180 O2SCL_ERR2(
"Rank is zero but the data vector has non-zero size ",
185 if (rk!=size.size()) {
186 O2SCL_ERR2(
"Rank does not match size vector size ",
192 for(
size_t i=0;i<
rk;i++) tot*=size[i];
194 O2SCL_ERR2(
"One entry in the size vector is zero ",
197 if (tot!=data.size()) {
198 O2SCL_ERR2(
"Product of size vector entries does not match data ",
245 template<
class size_vec_t>
246 void set(
const size_vec_t &index, data_t val) {
247 #if O2SCL_NO_RANGE_CHECK 252 if (index[0]>=size[0]) {
254 " greater than or equal to size[0]="+
255 szttos(size[0])+
" in tensor::set().").c_str(),
260 for(
size_t i=1;i<
rk;i++) {
261 #if O2SCL_NO_RANGE_CHECK 263 if (index[i]>=size[i]) {
265 szttos(index[i])+
" greater than or equal to size "+
266 szttos(size[i])+
" in tensor::set().").c_str(),
286 if (data.size()!=dat.size()) {
287 O2SCL_ERR2(
"Size of new vector does not equal tensor size in ",
298 template<
class size_vec_t> data_t &
get(
const size_vec_t &index) {
299 #if O2SCL_NO_RANGE_CHECK 304 if (index[0]>=size[0]) {
306 " greater than or equal to size[0]="+
307 szttos(size[0])+
" in tensor::get().").c_str(),
312 for(
size_t i=1;i<
rk;i++) {
313 #if O2SCL_NO_RANGE_CHECK 315 if (index[i]>=size[i]) {
317 szttos(index[i])+
" greater than or equal to size "+
318 szttos(size[i])+
" in tensor::get().").c_str(),
329 template<
class size_vec_t>
330 data_t
const &
get(
const size_vec_t &index)
const {
331 #if O2SCL_NO_RANGE_CHECK 336 if (index[0]>=size[0]) {
338 " greater than or equal to size[0]="+
339 szttos(size[0])+
" in tensor::get() (const).").c_str(),
344 for(
size_t i=1;i<
rk;i++) {
345 #if O2SCL_NO_RANGE_CHECK 347 if (index[i]>=size[i]) {
349 szttos(index[i])+
" greater than or equal to size "+
350 szttos(size[i])+
" in tensor::get() (const).").c_str(),
361 typedef boost::numeric::ublas::vector_slice<
363 typedef boost::numeric::ublas::slice slice;
389 template<
class size_vec_t>
393 " greater than or equal to rank "+
szttos(rk)+
394 " in tensor::vector_slice()").c_str(),
400 for(
size_t i=1;i<
rk;i++) {
402 if (i!=ix) start+=index[i];
405 for(
size_t i=ix+1;i<
rk;i++) stride*=size[i];
406 return ubvector_slice(data,slice(start,stride,size[ix]));
410 #ifdef O2SCL_NEVER_DEFINED 418 template<
class size_vec_t> ubmatrix_array matrix_slice() {
434 template<
class size_vec_t>
435 void resize(
size_t rank,
const size_vec_t &dim) {
436 for(
size_t i=0;i<rank;i++) {
438 O2SCL_ERR((((std::string)
"Requested zero size with non-zero ")+
439 "rank for index "+
szttos(i)+
" in tensor::"+
449 for(
size_t i=0;i<
rk;i++) {
468 " greater than or equal to rank "+
szttos(rk)+
469 " in tensor::get_size()").c_str(),
491 for(
size_t i=0;i<
rk;i++) tot*=size[i];
499 template<
class size_vec_t>
505 if (index[0]>=size[0]) {
507 " greater than or equal to size[0]="+
508 szttos(size[0])+
" in tensor::pack_indices().").c_str(),
512 for(
size_t i=1;i<
rk;i++) {
513 if (index[i]>=size[i]) {
515 szttos(index[i])+
" greater than or equal to size "+
516 szttos(size[i])+
" in tensor::pack_indices().").c_str(),
526 template<
class size_vec_t>
530 " greater than or equal to total size"+
532 " in tensor::unpack_index().").c_str(),
536 size_t ix2, sub_size;
537 for(
size_t i=0;i<
rk;i++) {
542 for(
size_t j=i+1;j<
rk;j++) sub_size*=size[j];
543 index[i]=ix/sub_size;
545 ix-=sub_size*(ix/sub_size);
563 size_t ix=o2scl::vector_min_index<vec_t,data_t>(
total_size(),
data);
571 void min(vec_size_t &index, data_t &val) {
587 size_t ix=o2scl::vector_max_index<vec_t,data_t>(
total_size(),
data);
595 void max(vec_size_t &index, data_t &val) {
611 size_t ix_min, ix_max;
612 o2scl::vector_minmax_index<vec_t,data_t>
622 void minmax(vec_size_t &index,
size_t &index_min, data_t &
min,
623 size_t &index_max, data_t &
max) {
624 size_t ix_min, ix_max;
635 if (rk==0)
return 0.0;
637 for(
size_t i=0;i<data.size();i++) {
650 (
size_t ix_x,
size_t ix_y,
table3d &tab, std::string x_name=
"x",
651 std::string y_name=
"y", std::string slice_name=
"z") {
657 if (nx==0 && ny==0) {
659 if (x_name.length()==0) x_name=
"x";
660 if (y_name.length()==0) y_name=
"y";
664 std::vector<double> grid_x(size[ix_x]), grid_y(size[ix_y]);
665 for(
size_t i=0;i<size[ix_x];i++) {
666 grid_x[i]=((double)i);
668 for(
size_t i=0;i<size[ix_y];i++) {
669 grid_y[i]=((double)i);
671 tab.
set_xy(
"x",grid_x.size(),grid_x,
672 "y",grid_y.size(),grid_y);
678 if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
680 "tensor_grid::convert_table3d_sum().",
exc_einval);
685 std::vector<size_t> ix;
688 tab.
set(ix[ix_x],ix[ix_y],slice_name,
689 tab.
get(ix[ix_x],ix[ix_y],slice_name)+
699 #ifdef O2SCL_NEVER_DEFINED 707 template<
class size_vec2_t,
class vec2_t>
708 tensor<> sum_slice(size_vec2_t &isum) {
710 if (isum.size()==0) {
712 "tensor::copy_slice_interp().",
715 if (this->rk<1+isum.size()) {
717 "tensor::copy_slice_interp().",
722 size_t rank_new=this->rk-isum.size();
725 std::vector<size_t> mapping;
726 std::vector<size_t> sz_new;
728 for(
size_t i=0;i<this->
rk;i++) {
730 for(
size_t j=0;j<isum.size();j++) {
731 if (isum[j]==i) found=
true;
734 sz_new.push_back(this->
get_size(i));
735 mapping.push_back(i);
748 std::vector<size_t> ix_old(this->rk), ix_new(rank_new);
750 for(
size_t j=0;j<rank_new;j++) {
751 ix_new=ix_old[mapping[j]];
756 tg_new.
set(ix_new,tg_new.
get(ix_new)+this->data[k]);
763 #ifdef O2SCL_NEVER_DEFINED 774 static const size_t index=1;
775 static const size_t fixed=2;
776 static const size_t sum=3;
777 static const size_t contract=4;
778 static const size_t reverse=5;
779 static const size_t range=6;
782 static const size_t interp=7;
783 static const size_t grid=8;
785 index_spec(
size_t typ,
size_t i1,
size_t i2,
double v) {
794 index_spec &ix_index(
size_t ix) {
795 return index_spec(index_spec::index,ix,0,0.0);
798 index_spec &ix_fixed(
size_t ix,
size_t ix2) {
799 return index_spec(index_spec::fixed,ix,ix2,0.0);
802 index_spec &ix_sum(
size_t ix) {
803 return index_spec(index_spec::sum,ix,0,0.0);
806 index_spec &ix_contract(
size_t ix,
size_t ix2) {
807 return index_spec(index_spec::contract,ix,ix2,0.0);
810 index_spec &ix_reverse(
size_t ix) {
811 return index_spec(index_spec::reverse,ix,0,0.0);
814 index_spec &ix_interp(
size_t ix,
double v) {
815 return index_spec(index_spec::interp,ix,0,v);
818 index_spec &ix_range(
size_t ix,
double v) {
819 return index_spec(index_spec::range,ix,0,v);
822 index_spec &ix_grid(
size_t ix,
double v) {
823 return index_spec(index_spec::grid,ix,0,v);
834 tensor<> copy(std::vector<index_spec> &spec) {
836 size_t rank_old=this->
rk;
839 std::vector<size_t> size_new;
841 size_t n_inner_loop=1;
843 for(
size_t i=0;i<spec.size();i++) {
844 if (spec[i].type==index_spec::index ||
845 spec[i].type==index_spec::reverse) {
849 size_new.push_back(this->size[spec[i].ix1]);
850 }
else if (spec_type==index_spec::contract) {
855 if (size[spec[i].ix1]<size[spec[i].ix2]) {
856 n_inner_loop*=size[spec[i].ix1];
858 n_inner_loop*=size[spec[i].ix2];
860 }
else if (spec[i].type==index_spec::sum) {
863 n_inner_loop*=size[spec[i].ix1];
864 }
else if (spec[i].type==index_spec::fixed) {
871 O2SCL_ERR(
"Zero new indices in tensor::copy().",
874 if (p.
valid()==
false) {
875 O2SCL_ERR(
"Not all indices accounted for in tensor::copy().",
883 std::vector<size_t> ix_new(rank_new);
884 std::vector<size_t> ix_old(rank_old);
894 for(
size_t j=0;j<n_inner_loop;j++) {
898 tg_new.set(ix_new,val);
910 template<
class data_t=
double,
class vec_t=std::vector<data_t>,
911 class vec_
size_t=std::vector<
size_t> >
class tensor1 :
912 public tensor<data_t,vec_t,vec_size_t> {
934 "tensor1::is_valid().",
945 data_t &
get(
size_t ix) {
950 const data_t &
get(
size_t ix)
const {
955 void set(
size_t index, data_t val)
963 template<
class size_vec_t>
964 void set(
const size_vec_t &index, data_t val) {
975 const data_t &
operator[](
size_t ix)
const {
return this->data[ix]; }
981 const data_t &
operator()(
size_t ix)
const {
return this->data[ix]; }
988 template<
class data_t=
double,
class vec_t=std::vector<data_t>,
989 class vec_
size_t=std::vector<
size_t> >
class tensor2 :
990 public tensor<data_t,vec_t,vec_size_t> {
1004 this->data.resize(tot);
1013 if (this->rk!=0 && this->rk!=2) {
1015 "tensor2::is_valid().",
1026 data_t &
get(
size_t ix1,
size_t ix2) {
1027 size_t sz[2]={ix1,ix2};
1032 const data_t &
get(
size_t ix1,
size_t ix2)
const {
1033 size_t sz[2]={ix1,ix2};
1038 void set(
size_t ix1,
size_t ix2, data_t val) {
1039 size_t sz[2]={ix1,ix2};
1049 template<
class size_vec_t>
1050 void set(
const size_vec_t &index, data_t val) {
1057 {
return this->data[ix*this->size[1]+iy]; }
1061 {
return this->data[ix*this->size[1]+iy]; }
1067 template<
class data_t=
double,
class vec_t=std::vector<data_t>,
1068 class vec_
size_t=std::vector<
size_t> >
class tensor3 :
1069 public tensor<data_t,vec_t,vec_size_t> {
1078 tensor<data_t,vec_t,vec_size_t>() {
1084 size_t tot=sz*sz2*sz3;
1085 this->data.resize(tot);
1094 if (this->rk!=0 && this->rk!=3) {
1096 "tensor3::is_valid().",
1107 data_t &
get(
size_t ix1,
size_t ix2,
size_t ix3) {
1108 size_t sz[3]={ix1,ix2,ix3};
1113 const data_t &
get(
size_t ix1,
size_t ix2,
size_t ix3)
const {
1114 size_t sz[3]={ix1,ix2,ix3};
1119 void set(
size_t ix1,
size_t ix2,
size_t ix3, data_t val) {
1120 size_t sz[3]={ix1,ix2,ix3};
1130 template<
class size_vec_t>
1131 void set(
const size_vec_t &index, data_t val) {
1141 template<
class data_t=
double,
class vec_t=std::vector<data_t>,
1142 class vec_
size_t=std::vector<
size_t> >
class tensor4 :
1143 public tensor<data_t,vec_t,vec_size_t> {
1151 tensor4(
size_t sz,
size_t sz2,
size_t sz3,
size_t sz4) :
1152 tensor<data_t,vec_t,vec_size_t>() {
1159 size_t tot=sz*sz2*sz3*sz4;
1160 this->data.resize(tot);
1169 if (this->rk!=0 && this->rk!=4) {
1171 "tensor4::is_valid().",
1182 data_t &
get(
size_t ix1,
size_t ix2,
size_t ix3,
size_t ix4) {
1183 size_t sz[4]={ix1,ix2,ix3,ix4};
1188 const data_t &
get(
size_t ix1,
size_t ix2,
size_t ix3,
1190 size_t sz[4]={ix1,ix2,ix3,ix4};
1195 void set(
size_t ix1,
size_t ix2,
size_t ix3,
size_t ix4,
1197 size_t sz[4]={ix1,ix2,ix3,ix4};
1207 template<
class size_vec_t>
1208 void set(
const size_vec_t &index, data_t val) {
1217 template<
class tensor_t>
1218 void tensor_out(std::ostream &os, tensor_t &t,
bool pretty=
true) {
1222 size_t rk=t.get_rank();
1223 os <<
"rank: " << rk <<
" sizes: ";
1224 for(
size_t i=0;i<
rk;i++) {
1225 os << t.get_size(i) <<
" ";
1228 auto &data=t.get_data();
1229 std::vector<size_t> ix(rk);
1230 std::vector<std::string> sv, sv_out;
1231 for(
size_t i=0;i<t.total_size();i++) {
1232 t.unpack_index(i,ix);
1233 std::string tmp=
"(";
1234 for(
size_t j=0;j<
rk;j++) {
1245 for(
size_t k=0;k<sv_out.size();k++) {
1246 os << sv_out[k] << std::endl;
1251 size_t rk=t.get_rank();
1253 for(
size_t i=0;i<
rk;i++) {
1254 os << t.get_size(i) <<
" ";
1257 auto &data=t.get_data();
1258 for(
size_t i=0;i<t.total_size();i++) {
1259 os << data[i] <<
" ";
1260 if (i%10==9) os << std::endl;
1269 #ifndef DOXYGEN_NO_O2NS void resize(size_t rank, const size_vec_t &dim)
Resize the tensor to rank rank with sizes given in dim.
void unpack_index(size_t ix, size_vec_t &index)
Unpack the single vector index into indices.
data_t min_value()
Compute the minimum value in the 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.
tensor3(size_t sz, size_t sz2, size_t sz3)
Create a rank 3 tensor of size (sz,sz2,sz3)
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
tensor1(size_t sz)
Create a rank 1 tensory of size sz.
void minmax(vec_size_t &index, size_t &index_min, data_t &min, size_t &index_max, data_t &max)
Compute the indices and values of the maximum and minimum in the tensor.
tensor4(size_t sz, size_t sz2, size_t sz3, size_t sz4)
Create a rank 4 tensor of size (sz,sz2,sz3,sz4)
tensor1()
Create an empty tensor.
sanity check failed - shouldn't happen
double & get(size_t ix, size_t iy, std::string name)
Get element in slice name at location ix,iy
size_t get_rank() const
Return the rank of the tensor.
invalid argument supplied by user
const data_t & operator[](size_t ix) const
Get an element using array-like indexing (const version)
tensor2()
Create an empty tensor.
A class for representing permutations.
void clear()
Clear the tensor of all data and free allocated memory.
ubvector_slice vector_slice(size_t ix, const size_vec_t &index)
Fix all but one index to create a vector.
void minmax_value(data_t &min, data_t &max)
Compute the minimum and maximum values in the tensor.
data_t max_value()
Compute the maximum value in the tensor.
void set_slice_all(std::string name, double val)
Set all of the values in slice name to val.
tensor2(size_t sz, size_t sz2)
Create a rank 2 tensor of size (sz,sz2)
size_t total_size() const
Returns the size of the tensor (the product of the sizes over every index)
void convert_table3d_sum(size_t ix_x, size_t ix_y, table3d &tab, std::string x_name="x", std::string y_name="y", std::string slice_name="z")
Convert to a o2scl::table3d object by summing over all but two indices.
void is_valid() const
Check that the o2scl::tensor4 object is valid.
void max(vec_size_t &index, data_t &val)
Compute the index and value of the maximum value in the tensor and return the maximum.
void min_index(vec_size_t &index)
Compute the index of the minimum value in the tensor.
void get_size(size_t &nx, size_t &ny) const
Get the size of the slices.
void is_valid() const
Check that the o2scl::tensor3 object is valid.
data_t & get(const size_vec_t &index)
Get the element indexed by index.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
tensor4()
Create an empty tensor.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
void max_index(vec_size_t &index)
Compute the index of the maximum value in the tensor.
data_t & operator()(size_t ix, size_t iy)
Get the element indexed by (ix1,ix2)
void is_valid() const
Check that the o2scl::tensor object is valid.
const vec_size_t & get_size_arr() const
Return the full vector of sizes.
const data_t & operator()(size_t ix) const
Get an element using operator() (const version)
void screenify(size_t nin, const string_arr_t &in_cols, std::vector< std::string > &out_cols, size_t max_size=80)
Reformat the columns for output of width size.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
data_t & operator[](size_t ix)
Get an element using array-like indexing.
const vec_t & get_data() const
Return the full data vector.
size_t pack_indices(const size_vec_t &index)
Pack the indices into a single vector index.
void is_valid() const
Check that the o2scl::tensor2 object is valid.
double total_sum() const
Return the sum over every element in the tensor.
void min(vec_size_t &index, data_t &val)
Compute the index of the minimum value in the tensor and return the minimum.
tensor3()
Create an empty tensor.
data_t & operator()(size_t ix)
Get an element using operator()
tensor(size_t rank, const size_vec_t &dim)
Create a tensor of rank rank with sizes given in dim.
const data_t & operator()(size_t ix, size_t iy) const
Get the element indexed by (ix1,ix2) (const version)
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.
A data structure containing one or more slices of two-dimensional data points defined on a grid...
void set(const size_vec_t &index, data_t val)
Set the element indexed by index to value val.
void tensor_out(std::ostream &os, tensor_t &t, bool pretty=true)
Output a tensor to a stream.
size_t get_size(size_t i) const
Returns the size of the ith index.
void swap_data(vec_t &dat)
Swap the data vector.
Invalid index for array or matrix.
tensor()
Create an empty tensor with zero rank.
vec_size_t size
A rank-sized vector of the sizes of each dimension.
bool valid() const
Check to see that a permutation is valid.
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.
std::string szttos(size_t x)
Convert a size_t to a string.
The default vector type from uBlas.
Interpolation class for general vectors.
Tensor class with arbitrary dimensions.
void is_valid() const
Check that the o2scl::tensor1 object is valid.
void set_all(data_t x)
Set all elements in a tensor to some fixed value.