23 #ifndef O2SCL_TENSOR_GRID_H
24 #define O2SCL_TENSOR_GRID_H
36 #include <gsl/gsl_matrix.h>
37 #include <gsl/gsl_ieee_utils.h>
39 #include <o2scl/err_hnd.h>
40 #include <o2scl/interp.h>
41 #include <o2scl/tensor.h>
42 #include <o2scl/table3d.h>
52 template<
class vec_t,
class vec_
size_t>
55 template<
class vec_t,
class vec_
size_t>
60 #ifndef DOXYGEN_NO_O2NS
162 template<
class vec_t=std::vector<
double>,
163 class vec_
size_t=std::vector<
size_t> >
class tensor_grid :
164 public tensor<double,vec_t,vec_size_t> {
168 #ifndef DOXYGEN_INTERNAL
172 #ifdef O2SCL_NEVER_DEFINED
203 template<
class size_vec_t>
205 tensor<double,vec_t,vec_size_t>(rank,dim) {
209 for(
size_t i=0;i<this->rk;i++) {
211 O2SCL_ERR((((std::string)
"Requested zero size with non-zero ")+
212 "rank for index "+
szttos(i)+
" in tensor_grid::"+
213 "tensor_grid(size_t,size_vec_t)").c_str(),
223 tensor<double,vec_t,vec_size_t>() {
227 for(
size_t j=0;j<this->rk;j++) {
228 this->size.push_back(ugs[j].get_npoints());
229 tot*=ugs[j].get_npoints();
231 this->data.resize(tot);
249 if (this->rk>0 && grid_set) {
251 for(
size_t i=0;i<this->rk;i++) {
255 if (tot2!=grid.size()) {
256 O2SCL_ERR2(
"Value grid_set is true but grid vector size ",
257 "is wrong in tensor_grid::is_valid().",
262 if (!grid_set && grid.size()>0) {
263 O2SCL_ERR2(
"Value grid_set is false but grid vector size ",
264 "is not zero in tensor_grid::is_valid().",
305 template<
class vec2_t>
306 void set_val(
const vec2_t &grdp,
double val) {
309 vec_size_t index(this->rk);
310 for(
size_t i=0;i<this->rk;i++) {
311 index[i]=lookup_grid(i,grdp[i]);
316 for(
size_t i=1;i<this->rk;i++) {
334 template<
class vec2_t,
class vec3_t>
335 void set_val(
const vec2_t &grdp,
double val, vec3_t &closest) {
338 vec_size_t index(this->rk);
339 for(
size_t i=0;i<this->rk;i++) {
340 index[i]=lookup_grid_val(i,grdp[i],closest[i]);
345 for(
size_t i=1;i<this->rk;i++) {
360 template<
class vec2_t>
double get_val(
const vec2_t &gridp) {
363 vec_size_t index(this->rk);
364 for(
size_t i=0;i<this->rk;i++) {
365 index[i]=lookup_grid(i,gridp[i]);
370 for(
size_t i=1;i<this->rk;i++) {
376 return this->data[ix];
385 template<
class vec2_t,
class vec3_t>
386 double get_val(
const vec2_t &gridp, vec3_t &closest) {
389 vec_size_t index(this->rk);
390 for(
size_t i=0;i<this->rk;i++) {
391 index[i]=lookup_grid_val(i,gridp[i],closest[i]);
396 for(
size_t i=1;i<this->rk;i++) {
402 return this->data[ix];
423 template<
class size_vec2_t>
424 void resize(
size_t rank,
const size_vec2_t &dim) {
427 for(
size_t i=0;i<rank;i++) {
429 O2SCL_ERR((((std::string)
"Requested zero size with non-zero ")+
430 "rank for index "+
szttos(i)+
" in tensor_grid::"+
437 this->size.resize(this->rk);
444 this->data.resize(0);
448 for(
size_t i=0;i<this->rk;i++) {
449 this->size[i]=dim[i];
452 this->data.resize(tot);
485 template<
class vec2_t>
488 O2SCL_ERR2(
"Tried to set grid for empty tensor in ",
489 "tensor_grid::set_grid_packed().",
exc_einval);
492 for(
size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
494 for(
size_t i=0;i<ngrid;i++) {
503 template<
class vec_vec_t>
506 O2SCL_ERR2(
"Tried to set grid for empty tensor in ",
510 for(
size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
513 for(
size_t i=0;i<this->rk;i++) {
514 for(
size_t j=0;j<this->size[i];j++) {
515 grid[k]=grid_vecs[i][j];
527 for(
size_t i=0;i<this->rk;i++) ngrid+=this->size[i];
530 for(
size_t i=0;i<this->rk;i++) {
531 for(
size_t j=0;j<this->size[i];j++) {
542 template<
class vec2_t>
544 if (grid_set==
false) {
549 O2SCL_ERR2(
"Tried to set grid for empty tensor in ",
553 for(
size_t i=0;i<this->rk;i++) {
554 for(
size_t j=0;j<this->size[i];j++) {
566 template<
class vec2_t>
568 if (grid_set==
false) {
570 "tensor_grid::set_grid_i_func().",
exc_einval);
573 O2SCL_ERR2(
"Tried to set grid for empty tensor in ",
574 "tensor_grid::set_grid_i_func().",
exc_einval);
578 std::map<std::string,double> vars;
579 calc.
compile(func.c_str(),&vars);
582 for(
size_t i=0;i<this->rk;i++) {
583 for(
size_t j=0;j<this->size[i];j++) {
585 vars[
"i"]=((double)j);
586 grid[k]=calc.
eval(&vars);
601 O2SCL_ERR2(
"Tried to set grid for empty tensor in ",
605 for(
size_t i=0;i<this->rk;i++) ngrid+=ugs[i].get_npoints();
608 for(
size_t i=0;i<this->rk;i++) {
609 for(
size_t j=0;j<ugs[i].get_npoints();j++) {
623 template<
class rvec_t>
void copy_grid(
size_t i, rvec_t &v)
const {
624 v.resize(this->size[i]);
626 for(
size_t k=0;k<i;k++) istart+=this->size[k];
627 for(
size_t k=0;k<this->size[i];k++) {
636 O2SCL_ERR(
"Grid not set in tensor_grid::get_grid().",
641 " greater than or equal to rank, "+
szttos(this->rk)+
642 ", in tensor_grid::get_grid().").c_str(),
646 for(
size_t k=0;k<i;k++) istart+=this->size[k];
647 return grid[istart+j];
653 O2SCL_ERR(
"Grid not set in tensor_grid::get_grid().",
658 " greater than or equal to rank, "+
szttos(this->rk)+
659 ", in tensor_grid::get_grid().").c_str(),
663 for(
size_t k=0;k<i;k++) istart+=this->size[k];
677 " greater than or equal to rank, "+
szttos(this->rk)+
678 ", in tensor_grid::lookup_grid_val().").c_str(),
681 if (grid_set==
false) {
682 O2SCL_ERR(
"Grid not set in tensor_grid::lookup_grid_val().",
691 for(
size_t j=0;j<i;j++) istart+=this->size[j];
693 double min=fabs(grid[istart]-temp);
695 for(
size_t j=istart;j<istart+this->size[i];j++) {
696 if (fabs(grid[j]-temp)<min) {
698 min=fabs(grid[j]-temp);
708 return lookup_grid_val(i,val,val2);
721 template<
class vec2_t,
class size_vec2_t>
723 for(
size_t k=0;k<this->rk;k++) {
724 indices[k]=lookup_grid(k,vals[k]);
739 O2SCL_ERR(
"Grid not set in tensor_grid::lookup_grid_packed_val().",
745 ", in tensor_grid::lookup_grid_packed_val().").c_str(),
749 for(
size_t j=0;j<i;j++) istart+=this->size[j];
751 double min=fabs(grid[istart]-val);
753 for(
size_t j=istart;j<istart+this->size[i];j++) {
754 if (fabs(grid[j]-val)<min) {
756 min=fabs(grid[j]-val);
768 return lookup_grid_packed_val(i,val,val2);
777 template<
class size_vec2_t,
class vec2_t>
780 if (this->rk<1+ifix.size()) {
782 "tensor_grid::copy_slice_interp().",
785 if (ifix.size()!=vals.size()) {
786 O2SCL_ERR2(
"Mismatch between indices and values in ",
787 "tensor_grid::copy_slice_interp().",
792 size_t rank_new=this->rk-ifix.size();
795 std::vector<size_t> sz_new;
796 std::vector<std::vector<double> > grid_new;
797 for(
size_t i=0;i<this->rk;i++) {
799 for(
size_t j=0;j<ifix.size();j++) {
800 if (ifix[j]==i) found=
true;
803 sz_new.push_back(this->get_size(i));
804 std::vector<double> grid_temp;
805 for(
size_t j=0;j<this->get_size(i);j++) {
806 grid_temp.push_back(this->get_grid(i,j));
808 grid_new.push_back(grid_temp);
817 std::vector<size_t> ix_new(rank_new);
818 std::vector<double> point_old(this->rk);
821 for(
size_t i=0;i<tg_new.total_size();i++) {
824 tg_new.unpack_index(i,ix_new);
827 for(
size_t j=0;j<this->rk;j++) {
829 for(
size_t k=0;k<ifix.size();k++) {
830 if (ifix[k]==j) ix_found=k;
833 point_old[j]=this->get_grid(j,ix_new[j]);
835 point_old[j]=vals[ix_found];
853 void convert_table3d_sum
854 (
size_t ix_x,
size_t ix_y,
table3d &tab, std::string x_name=
"x",
855 std::string y_name=
"y", std::string slice_name=
"z")
const {
861 if (nx==0 && ny==0) {
863 if (x_name.length()==0) x_name=
"x";
864 if (y_name.length()==0) y_name=
"y";
868 std::vector<double> grid_x, grid_y;
869 copy_grid(ix_x,grid_x);
870 copy_grid(ix_y,grid_y);
871 tab.
set_xy(
"x",grid_x.size(),grid_x,
872 "y",grid_y.size(),grid_y);
878 if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
880 "tensor_grid::convert_table3d_sum().",
exc_einval);
885 std::vector<size_t> ix;
886 for(
size_t i=0;i<this->total_size();i++) {
887 this->unpack_index(i,ix);
888 tab.
set(ix[ix_x],ix[ix_y],slice_name,
889 tab.
get(ix[ix_x],ix[ix_y],slice_name)+
920 template<
class size_vec2_t>
922 table3d &tab, std::string slice_name=
"z")
const {
924 if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
925 O2SCL_ERR2(
"Either indices greater than rank or x and y ind",
926 "ices equal in tensor_grid::copy_table3d_align().",
935 if (nx!=this->size[ix_x] || ny!=this->size[ix_y]) {
937 "tensor_grid::copy_table3d_align().",
exc_einval);
945 for(
size_t i=0;i<nx;i++) {
946 for(
size_t j=0;j<ny;j++) {
949 double val=this->get(index);
950 tab.
set(i,j,slice_name,val);
960 template<
class size_vec2_t>
961 void copy_table3d_align_setxy
962 (
size_t ix_x,
size_t ix_y, size_vec2_t &index,
963 table3d &tab, std::string x_name=
"x", std::string y_name=
"y",
964 std::string slice_name=
"z")
const {
970 if (nx==0 && ny==0) {
972 if (x_name.length()==0) x_name=
"x";
973 if (y_name.length()==0) y_name=
"y";
977 std::vector<double> grid_x, grid_y;
978 copy_grid(ix_x,grid_x);
979 copy_grid(ix_y,grid_y);
980 tab.
set_xy(
"x",grid_x.size(),grid_x,
981 "y",grid_y.size(),grid_y);
986 copy_table3d_align(ix_x,ix_y,index,tab,slice_name);
1011 template<
class size_vec2_t>
1013 table3d &tab, std::string slice_name=
"z")
const {
1015 if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
1016 O2SCL_ERR2(
"Either indices greater than rank or x and y ",
1017 "indices equal in tensor_grid::copy_table3d_interp().",
1025 if (nx==0 && ny==0) {
1027 return copy_table3d_align(ix_x,ix_y,index,tab,slice_name);
1031 std::vector<double> vals(this->rk);
1032 for(
size_t i=0;i<this->rk;i++) {
1033 if (i!=ix_x && i!=ix_y) vals[i]=this->get_grid(i,index[i]);
1041 for(
size_t i=0;i<nx;i++) {
1042 for(
size_t j=0;j<ny;j++) {
1054 template<
class vec2_t>
1057 std::string slice_name=
"z",
1058 int verbose=0)
const {
1060 if (ix_x>=this->rk || ix_y>=this->rk || ix_x==ix_y) {
1061 O2SCL_ERR2(
"Either indices greater than rank or x and y ",
1062 "indices equal in tensor_grid::copy_table3d_interp().",
1065 if (values.size()!=this->rk) {
1066 O2SCL_ERR2(
"Values array not equal to rank ",
1067 "in tensor_grid::copy_table3d_interp_values().",
1085 for(
size_t i=0;i<nx;i++) {
1086 for(
size_t j=0;j<ny;j++) {
1091 std::cout <<
"At location values: ";
1092 for(
size_t k=0;k<values.size();k++) {
1093 std::cout << values[k] <<
" ";
1095 std::cout <<
"Interpolated to get: "
1096 << i <<
" " << j <<
" " << slice_name <<
" "
1112 template<
class vec2_t>
1113 void copy_table3d_interp_values_setxy
1114 (
size_t ix_x,
size_t ix_y, vec2_t &values,
table3d &tab,
1115 std::string x_name=
"x", std::string y_name=
"y",
1116 std::string slice_name=
"z")
const {
1122 if (nx==0 && ny==0) {
1124 if (x_name.length()==0) x_name=
"x";
1125 if (y_name.length()==0) y_name=
"y";
1129 std::vector<double> grid_x, grid_y;
1130 copy_grid(ix_x,grid_x);
1131 copy_grid(ix_y,grid_y);
1132 tab.
set_xy(
"x",grid_x.size(),grid_x,
1133 "y",grid_y.size(),grid_y);
1138 copy_table3d_interp_values(ix_x,ix_y,values,tab,slice_name);
1197 interp_t si(this->size[0],grid,this->data,itype);
1198 return si.eval(vals[0]);
1204 for(
size_t i=1;i<this->rk;i++) ss*=this->size[i];
1207 std::vector<vec_t> yvec(ss);
1208 std::vector<interp_t *> si(ss);
1209 for(
size_t i=0;i<ss;i++) {
1210 yvec[i].resize(this->size[0]);
1215 index_range_t size_new(this->size,
ub_range(1,this->rk));
1216 tdat.
resize(this->rk-1,size_new);
1219 data_range_t grid_new(grid,
ub_range(this->size[0],grid.size()));
1223 vec_size_t co(this->rk);
1224 for(
size_t i=0;i<this->rk;i++) co[i]=0;
1229 while(done==
false) {
1232 for(
size_t i=0;i<this->size[0];i++) {
1234 yvec[cnt][i]=this->get(co);
1237 si[cnt]=
new interp_t(this->size[0],grid,yvec[cnt],itype);
1239 index_range_t co2(co,
ub_range(1,this->rk));
1240 tdat.set(co2,si[cnt]->eval(vals[0]));
1246 for(
int j=((
int)this->rk)-1;j>0;j--) {
1247 if (co[j]>=this->size[j]) {
1254 if (cnt==ss) done=
true;
1262 for(
size_t i=0;i<ss;i++) {
1287 template<
class vec2_
size_t,
class vec3_
size_t,
class vec2_t>
1288 double interp_linear_partial
1289 (
const vec2_size_t &ix_to_interp,
1290 vec3_size_t &ix,
const vec2_t &val)
const {
1292 if (val.size()!=ix_to_interp.size()) {
1293 O2SCL_ERR2(
"Index and value list don't match in ",
1294 "tensor_grid::interp_linear_partial().",
1297 if (ix_to_interp.size()>this->get_rank() ||
1298 ix_to_interp.size()==0) {
1299 O2SCL_ERR2(
"Index list too large or too small in ",
1300 "tensor_grid::interp_linear_partial().",
1306 std::vector<size_t> loc(ix_to_interp.size());
1307 std::vector<double> gnew;
1308 for(
size_t i=0;i<ix_to_interp.size();i++) {
1309 size_t ixi=ix_to_interp[i];
1310 if (ixi>=this->get_rank()) {
1311 O2SCL_ERR2(
"Index to interpolate larger than tensor rank in ",
1312 "tensor_grid::interp_linear_partial().",
1315 std::vector<double> grid_one(this->size[ixi]);
1316 for(
size_t j=0;j<this->size[ixi];j++) {
1317 grid_one[j]=this->get_grid(ixi,j);
1320 loc[i]=sv.
find(val[i]);
1321 gnew.push_back(grid_one[loc[i]]);
1322 gnew.push_back(grid_one[loc[i]+1]);
1327 std::vector<size_t> snew(ix_to_interp.size());
1328 for(
size_t i=0;i<ix_to_interp.size();i++) {
1335 for(
size_t i=0;i<tnew.total_size();i++) {
1336 std::vector<size_t> index_new(ix_to_interp.size());
1337 tnew.unpack_index(i,index_new);
1338 for(
size_t j=0;j<ix_to_interp.size();j++) {
1339 ix[ix_to_interp[j]]=index_new[j]+loc[j];
1341 tnew.set(index_new,this->get(ix));
1370 std::vector<size_t> loc(this->rk);
1371 std::vector<double> gnew;
1372 for(
size_t i=0;i<this->rk;i++) {
1373 std::vector<double> grid_unpacked(this->size[i]);
1374 for(
size_t j=0;j<this->size[i];j++) {
1375 grid_unpacked[j]=grid[j+rgs];
1378 loc[i]=sv.
find(v[i]);
1379 gnew.push_back(grid_unpacked[loc[i]]);
1380 gnew.push_back(grid_unpacked[loc[i]+1]);
1386 std::vector<size_t> snew(this->rk);
1387 for(
size_t i=0;i<this->rk;i++) {
1394 for(
size_t i=0;i<tnew.total_size();i++) {
1395 std::vector<size_t> index_new(this->rk), index_old(this->rk);
1396 tnew.unpack_index(i,index_new);
1397 for(
size_t j=0;j<this->rk;j++) index_old[j]=index_new[j]+loc[j];
1398 tnew.set(index_new,this->get(index_old));
1416 template<
class vec2_t>
1420 return this->data[0]+(this->data[1]-this->data[0])/
1421 (grid[1]-grid[0])*(v[0]-grid[0]);
1424 size_t last=this->rk-1;
1425 double frac=(v[last]-get_grid(last,0))/
1426 (get_grid(last,1)-get_grid(last,0));
1434 for(
size_t i=0;i<tnew.total_size();i++) {
1435 std::vector<size_t> index(this->rk);
1436 tnew.unpack_index(i,index);
1437 index[this->rk-1]=0;
1438 double val_lo=this->get(index);
1439 index[this->rk-1]=1;
1440 double val_hi=this->get(index);
1441 tnew.set(index,val_lo+frac*(val_hi-val_lo));
1459 template<
class vec2_t,
class vec3_t>
1464 std::vector<size_t> loc(this->rk);
1465 std::vector<double> gnew;
1466 for(
size_t i=0;i<this->size[0];i++) {
1467 gnew.push_back(grid[i]);
1471 for(
size_t i=1;i<this->rk;i++) {
1472 std::vector<double> grid_unpacked(this->size[i]);
1473 for(
size_t j=0;j<this->size[i];j++) {
1474 grid_unpacked[j]=grid[j+rgs];
1477 loc[i]=sv.
find(v[i]);
1478 gnew.push_back(grid_unpacked[loc[i]]);
1479 gnew.push_back(grid_unpacked[loc[i]+1]);
1485 std::vector<size_t> snew(this->rk);
1486 snew[0]=this->size[0];
1487 for(
size_t i=1;i<this->rk;i++) {
1494 for(
size_t i=0;i<tnew.total_size();i++) {
1495 std::vector<size_t> index_new(this->rk), index_old(this->rk);
1496 tnew.unpack_index(i,index_new);
1497 for(
size_t j=0;j<this->rk;j++) {
1498 index_old[j]=index_new[j]+loc[j];
1500 tnew.set(index_new,this->get(index_old));
1520 template<
class vec2_t,
class vec3_t>
1524 size_t n=this->size[0];
1526 vec_size_t ix0(2), ix1(2);
1529 for(
size_t i=0;i<n;i++) {
1532 res[i]=this->get(ix0)+(this->get(ix1)-this->get(ix0))/
1533 (grid[n+1]-grid[n])*(v[1]-grid[n]);
1538 size_t last=this->rk-1;
1539 double frac=(v[last]-get_grid(last,0))/
1540 (get_grid(last,1)-get_grid(last,0));
1548 for(
size_t i=0;i<tnew.total_size();i++) {
1549 std::vector<size_t> index(this->rk);
1550 tnew.unpack_index(i,index);
1551 index[this->rk-1]=0;
1552 double val_lo=this->get(index);
1553 index[this->rk-1]=1;
1554 double val_hi=this->get(index);
1555 tnew.set(index,val_lo+frac*(val_hi-val_lo));
1574 template<
class vec2_t,
class vec3_t>
1577 size_t n=this->size[ifree];
1581 std::vector<size_t> map;
1582 map.push_back(ifree);
1583 for(
size_t i=0;i<this->rk;i++) {
1591 vec_size_t loc(this->rk);
1593 for(
size_t i=0;i<this->rk;i++) {
1594 vec_t grid_unpacked(this->size[i]);
1595 for(
size_t j=0;j<this->size[i];j++) {
1596 grid_unpacked[j]=grid[j+rgs];
1600 loc[i]=sv.
find(v[i]);
1606 std::vector<double> gnew, vnew;
1607 for(
size_t new_ix=0;new_ix<this->rk;new_ix++) {
1608 for(
size_t old_ix=0;old_ix<this->rk;old_ix++) {
1609 if (map[new_ix]==old_ix) {
1610 vnew.push_back(v[old_ix]);
1611 if (old_ix==ifree) {
1612 for(
size_t j=0;j<this->size[old_ix];j++) {
1613 gnew.push_back(this->get_grid(old_ix,j));
1616 gnew.push_back(this->get_grid(old_ix,loc[old_ix]));
1617 gnew.push_back(this->get_grid(old_ix,loc[old_ix]+1));
1627 std::vector<size_t> snew;
1629 for(
size_t i=0;i<this->rk;i++) {
1640 for(
size_t i=0;i<tnew.total_size();i++) {
1641 std::vector<size_t> index_new(this->rk), index_old(this->rk);
1642 tnew.unpack_index(i,index_new);
1643 for(
size_t j=0;j<this->rk;j++) {
1644 index_old[map[j]]=index_new[j]+loc[map[j]];
1646 tnew.set(index_new,this->get(index_old));
1656 template<
class vecf_t,
class vecf_
size_t>
friend void o2scl_hdf::hdf_output
1674 bool err_on_fail=
true)
const {
1677 size_t rank_old=this->rk;
1681 std::vector<size_t> size_new;
1685 std::vector<index_spec> spec_old(rank_old);
1686 std::vector<index_spec> spec_new;
1689 size_t n_sum_loop=1;
1694 std::vector<size_t> sum_sizes;
1697 std::vector<size_t> ix_to_interp;
1700 std::vector<double> new_grid;
1703 for(
size_t i=0;i<spec.size();i++) {
1706 if (spec[i].ix1>=rank_old) {
1708 O2SCL_ERR2(
"Index too large (index,reverse) in ",
1709 "tensor_grid::rearrange_and_copy().",
1713 std::cout <<
"Index too large (index,reverse) in "
1714 <<
"tensor_grid::rearrange_and_copy()."
1720 size_new.push_back(this->size[spec[i].ix1]);
1723 spec_old[spec[i].ix1]=
index_spec(spec[i].type,
1734 for(
size_t k=0;k<this->get_size(spec[i].ix1);k++) {
1735 new_grid.push_back(this->get_grid(spec[i].ix1,k));
1738 size_t nt=this->get_size(spec[i].ix1);
1739 for(
size_t k=0;k<nt;k++) {
1740 new_grid.push_back(this->get_grid(spec[i].ix1,nt-1-k));
1744 if (spec[i].ix1>=rank_old ||
1745 spec[i].ix2>=this->size[spec[i].ix1] ||
1746 spec[i].ix3>=this->size[spec[i].ix1]) {
1752 std::cout <<
"Index too large (range) in "
1753 <<
"tensor in tensor::rearrange_and_copy()."
1759 if (spec[i].ix3>spec[i].ix2) {
1760 size_new.push_back(spec[i].ix3-spec[i].ix2+1);
1762 size_new.push_back(spec[i].ix2-spec[i].ix3+1);
1766 spec_old[spec[i].ix1]=
1767 index_spec(spec[i].type,rank_new,spec[i].ix2,
1768 spec[i].ix3,spec[i].val1);
1771 spec[i].ix2,spec[i].ix3,spec[i].val1));
1774 if (spec[i].ix3>spec[i].ix2) {
1775 for(
size_t k=0;k<spec[i].ix3-spec[i].ix2+1;k++) {
1776 new_grid.push_back(this->get_grid(spec[i].ix1,k+spec[i].ix2));
1779 for(
size_t k=0;k<spec[i].ix2-spec[i].ix3+1;k++) {
1780 new_grid.push_back(this->get_grid(spec[i].ix1,k+spec[i].ix3));
1784 if (spec[i].ix1>=rank_old || spec[i].ix2>=rank_old) {
1787 "tensor_grid::rearrange_and_copy().",
1791 std::cout <<
"Index too large (trace) in "
1792 <<
"tensor_grid::rearrange_and_copy()."
1798 if (this->size[spec[i].ix1]<this->size[spec[i].ix2]) {
1799 n_sum_loop*=this->size[spec[i].ix1];
1800 sum_sizes.push_back(this->size[spec[i].ix1]);
1802 n_sum_loop*=this->size[spec[i].ix2];
1803 sum_sizes.push_back(this->size[spec[i].ix2]);
1807 spec_old[spec[i].ix1]=
index_spec(spec[i].type,spec[i].ix1,
1809 spec_old[spec[i].ix2]=
index_spec(spec[i].type,spec[i].ix2,
1812 if (spec[i].ix1>=rank_old) {
1815 "tensor_grid::rearrange_and_copy().",
1819 std::cout <<
"Index " << spec[i].ix1
1820 <<
" too large (sum) in "
1821 <<
"tensor_grid::rearrange_and_copy()."
1827 n_sum_loop*=this->size[spec[i].ix1];
1828 sum_sizes.push_back(this->size[spec[i].ix1]);
1829 spec_old[spec[i].ix1]=
index_spec(spec[i].type,
1834 if (spec[i].ix1>=rank_old ||
1835 spec[i].ix2>=this->size[spec[i].ix1]) {
1838 "tensor_grid::rearrange_and_copy().",
1842 std::cout <<
"Index " << spec[i].ix1 <<
" or "
1843 << spec[i].ix2 <<
" too large (fixed) in "
1844 <<
"tensor_grid::rearrange_and_copy()."
1852 spec_old[spec[i].ix1]=
index_spec(spec[i].type,rank_new,
1855 if (spec[i].ix1>=rank_old) {
1858 "tensor_grid::rearrange_and_copy().",
1862 std::cout <<
"Index " << spec[i].ix1
1863 <<
" too large (interp) in "
1864 <<
"tensor_grid::rearrange_and_copy()."
1872 spec_old[spec[i].ix1]=
index_spec(spec[i].type,
1877 ix_to_interp.push_back(spec[i].ix1);
1880 if (spec[i].ix1>=rank_old) {
1883 "tensor_grid::rearrange_and_copy().",
1887 std::cout <<
"Index " << spec[i].ix1
1888 <<
" too large (grid) in "
1889 <<
"tensor_grid::rearrange_and_copy()."
1898 if (spec[i].ix3==1) {
1903 rat=pow(spec[i].val2/spec[i].val1,1.0/((
double)spec[i].ix2));
1905 curr_grid_size=spec[i].ix2+1;
1908 curr_grid_size=((size_t)(log(spec[i].val2/spec[i].val1)/
1909 log(spec[i].val3)));
1911 size_new.push_back(curr_grid_size);
1912 for(
int k=0;k<curr_grid_size;k++) {
1914 new_grid.push_back(spec[i].val1);
1915 }
else if (curr_grid_size>1 && k==curr_grid_size-1) {
1916 new_grid.push_back(spec[i].val2);
1918 new_grid.push_back(spec[i].val1*pow(rat,((
double)k)));
1926 width=(spec[i].val2-spec[i].val1)/((
double)spec[i].ix2);
1928 curr_grid_size=spec[i].ix2+1;
1931 curr_grid_size=((size_t)(1+(spec[i].val2-
1932 spec[i].val1)/spec[i].val3));
1934 size_new.push_back(curr_grid_size);
1935 for(
int k=0;k<curr_grid_size;k++) {
1937 new_grid.push_back(spec[i].val1);
1938 }
else if (curr_grid_size>1 && k==curr_grid_size-1) {
1939 new_grid.push_back(spec[i].val2);
1941 new_grid.push_back(spec[i].val1+width*((
double)k));
1949 (spec[i].type,rank_new,spec[i].ix2,spec[i].ix3,
1950 spec[i].val1,spec[i].val2,spec[i].val3);
1952 (spec[i].type,spec[i].ix1,spec[i].ix2,
1953 spec[i].ix3,spec[i].val1,
1954 spec[i].val2,spec[i].val3));
1957 ix_to_interp.push_back(spec[i].ix1);
1963 O2SCL_ERR2(
"Index specification type not allowed in ",
1964 "tensor_grid::rearrange_and_copy()",
1968 std::cout <<
"Index specification type not allowed in "
1969 <<
"tensor_grid::rearrange_and_copy()." << std::endl;
1975 size_t n_sums=sum_sizes.size();
1981 "tensor_grid::rearrange_and_copy()",
1985 std::cout <<
"Zero new indices in "
1986 <<
"tensor_grid::rearrange_and_copy()." << std::endl;
1991 for(
size_t i=0;i<rank_old;i++) {
1994 O2SCL_ERR2(
"Not all indices accounted for in ",
1995 "tensor_grid::rearrange_and_copy()",
1999 std::cout <<
"Index " << i <<
" not accounted for in "
2000 <<
"tensor_grid::rearrange_and_copy()." << std::endl;
2009 std::cout <<
"Using a " << rank_old
2010 <<
" rank tensor to create a new "
2011 << rank_new <<
" rank tensor." << std::endl;
2014 for(
size_t i=0;i<rank_old;i++) {
2015 std::cout <<
"Old index " << i;
2017 std::cout <<
" is being remapped to new index "
2018 << spec_old[i].ix1 <<
"." << std::endl;
2020 std::cout <<
" is being remapped to new index "
2021 << spec_old[i].ix1 <<
" with a range from "
2022 << spec_old[i].ix2 <<
" to " << spec_old[i].ix3
2023 <<
"." << std::endl;
2025 std::cout <<
" is being reversed and remapped to new index "
2026 << spec_old[i].ix1 <<
"." << std::endl;
2028 std::cout <<
" is being traced with index "
2029 << spec_old[i].ix2 <<
"." << std::endl;
2031 std::cout <<
" is being summed." << std::endl;
2033 std::cout <<
" is being fixed to " << spec_old[i].ix2
2034 <<
"." << std::endl;
2036 std::cout <<
" is being interpolated from value "
2037 << spec_old[i].val1 <<
"." << std::endl;
2040 std::cout <<
" is being reinterpolated based on grid\n "
2041 << spec_old[i].val1 <<
" "
2042 << spec_old[i].val2 <<
" ";
2044 std::cout << spec_old[i].ix2;
2046 std::cout << spec_old[i].val3;
2048 if (spec_old[i].ix3==1) {
2049 std::cout <<
" (log)." << std::endl;
2051 std::cout <<
"." << std::endl;
2055 for(
size_t i=0;i<rank_new;i++) {
2056 std::cout <<
"New index " << i;
2058 std::cout <<
" was remapped from old index " << spec_new[i].ix1
2059 <<
"." << std::endl;
2061 std::cout <<
" was remapped from old index " << spec_new[i].ix1
2062 <<
" using range from " << spec_new[i].ix2 <<
" to "
2063 << spec_new[i].ix3 <<
"." << std::endl;
2065 std::cout <<
" was reversed and remapped from old index "
2066 << spec_new[i].ix1 <<
"." << std::endl;
2069 std::cout <<
" was obtained from grid\n "
2070 << spec_new[i].val1 <<
" "
2071 << spec_new[i].val2 <<
" ";
2073 std::cout << spec_new[i].ix2;
2075 std::cout << spec_new[i].val3;
2077 if (spec_new[i].ix3==1) {
2078 std::cout <<
" (log)." << std::endl;
2080 std::cout <<
"." << std::endl;
2090 std::cout <<
"New grid is: " << std::endl;
2091 for(
size_t k=0;k<rank_new;k++) {
2092 std::cout << k <<
" (" << t_new.get_size(k) <<
"): ";
2093 if (t_new.get_size(k)>3) {
2094 std::cout << t_new.
get_grid(k,0) <<
", ";
2095 std::cout << t_new.
get_grid(k,1) <<
" ... ";
2097 for(
int ell=0;ell<((int)t_new.get_size(k))-1;ell++) {
2098 std::cout << t_new.
get_grid(k,ell) <<
", ";
2101 std::cout << t_new.
get_grid(k,t_new.get_size(k)-1) << std::endl;
2103 std::cout <<
"Interpolations (" << n_interps <<
"): ";
2108 std::vector<size_t> ix_new(rank_new);
2109 std::vector<size_t> ix_old(rank_old);
2110 std::vector<size_t> sum_ix(n_sums);
2113 for(
size_t i=0;i<t_new.total_size();i++) {
2116 t_new.unpack_index(i,ix_new);
2119 std::vector<double> interp_vals;
2122 for(
size_t j=0;j<rank_old;j++) {
2124 ix_old[j]=ix_new[spec_old[j].ix1];
2126 if (spec_old[j].ix2<spec_old[j].ix3) {
2127 ix_old[j]=ix_new[spec_old[j].ix1]+spec_old[j].ix2;
2129 ix_old[j]=spec_old[j].ix2-ix_new[spec_old[j].ix1];
2132 ix_old[j]=this->get_size(j)-1-ix_new[spec_old[j].ix1];
2134 ix_old[j]=spec_old[j].ix2;
2136 interp_vals.push_back(spec_old[j].val1);
2138 if (spec_old[j].ix3==1) {
2139 double val=spec_old[j].val1*
2140 pow(spec_old[j].val3,ix_new[spec_old[j].ix1]);
2141 interp_vals.push_back(val);
2143 double val=spec_old[j].val1+
2144 ix_new[spec_old[j].ix1]*spec_old[j].val3;
2145 interp_vals.push_back(val);
2148 if (spec_old[j].ix3==1) {
2149 double val=spec_old[j].val1*
2150 pow(spec_old[j].val3,ix_new[spec_old[j].ix1]);
2151 interp_vals.push_back(val);
2153 double val=spec_old[j].val1+
2154 ix_new[spec_old[j].ix1]*spec_old[j].val3;
2155 interp_vals.push_back(val);
2163 for(
size_t j=0;j<n_sum_loop;j++) {
2167 size_t j2=j, sub_size;
2168 for(
size_t k=0;k<n_sums;k++) {
2173 for(
size_t kk=k+1;kk<n_sums;kk++) sub_size*=sum_sizes[kk];
2174 sum_ix[k]=j2/sub_size;
2175 j2-=sub_size*(j2/sub_size);
2180 std::cout <<
"n_sum_loop: " << n_sum_loop <<
" n_sums: "
2181 << n_sums <<
" sum_sizes: ";
2183 std::cout <<
"j: " << j <<
" sum_ix: ";
2189 for(
size_t k=0;k<rank_old;k++) {
2191 if (cnt>=sum_ix.size()) {
2192 std::cout <<
"X: " << cnt <<
" " << sum_ix.size()
2195 "tensor_grid::rearrange_and_copy()",
2198 ix_old[k]=sum_ix[cnt];
2201 spec_old[k].ix1<spec_old[k].ix2) {
2202 if (cnt>=sum_ix.size()) {
2203 std::cout <<
"X: " << cnt <<
" " << sum_ix.size()
2206 "tensor_grid::rearrange_and_copy()",
2209 ix_old[spec_old[k].ix1]=sum_ix[cnt];
2210 ix_old[spec_old[k].ix2]=sum_ix[cnt];
2216 std::cout <<
"Here old: ";
2218 std::cout <<
"Here new: ";
2222 val+=this->interp_linear_partial(ix_to_interp,ix_old,interp_vals);
2224 val+=this->get(ix_old);
2230 t_new.set(ix_new,val);
2242 template<
class vec_t=std::vector<
double>,
2256 this->data.resize(sz);
2257 this->grid_set=
false;
2260 #ifdef O2SCL_NEVER_DEFINED
2264 virtual ~tensor_grid1() {
2274 const double &
get(
size_t ix1)
const {
2280 void set(
size_t ix1,
double val) {
2290 <range_t,data_range_t,index_range_t>(&x);
2303 template<
class vec_t=std::vector<
double>,
2319 this->data.resize(tot);
2320 this->grid_set=
false;
2323 #ifdef O2SCL_NEVER_DEFINED
2327 virtual ~tensor_grid2() {
2331 double &
get(
size_t ix1,
size_t ix2) {
2332 size_t sz[2]={ix1,ix2};
2337 const double &
get(
size_t ix1,
size_t ix2)
const {
2338 size_t sz[2]={ix1,ix2};
2343 void set(
size_t ix1,
size_t ix2,
double val) {
2344 size_t sz[2]={ix1,ix2};
2353 double arr[2]={x,y};
2355 <range_t,data_range_t,index_range_t>(arr);
2360 double arr[2]={x,y};
2368 template<
class vec_t=std::vector<
double>,
2385 size_t tot=sz*sz2*sz3;
2386 this->data.resize(tot);
2387 this->grid_set=
false;
2390 #ifdef O2SCL_NEVER_DEFINED
2394 virtual ~tensor_grid3() {
2398 double &
get(
size_t ix1,
size_t ix2,
size_t ix3) {
2399 size_t sz[3]={ix1,ix2,ix3};
2404 const double &
get(
size_t ix1,
size_t ix2,
size_t ix3)
const {
2405 size_t sz[3]={ix1,ix2,ix3};
2410 void set(
size_t ix1,
size_t ix2,
size_t ix3,
double val) {
2411 size_t sz[3]={ix1,ix2, ix3};
2419 double interp(
double x,
double y,
double z) {
2420 double arr[3]={x,y,z};
2422 <range_t,data_range_t,index_range_t>(arr);
2427 double arr[3]={x,y,z};
2435 template<
class vec_t=std::vector<
double>,
2453 size_t tot=sz*sz2*sz3*sz4;
2454 this->data.resize(tot);
2455 this->grid_set=
false;
2458 #ifdef O2SCL_NEVER_DEFINED
2462 virtual ~tensor_grid4() {
2466 double &
get(
size_t ix1,
size_t ix2,
size_t ix3,
size_t ix4) {
2467 size_t sz[4]={ix1,ix2,ix3,ix4};
2472 const double &
get(
size_t ix1,
size_t ix2,
size_t ix3,
2474 size_t sz[4]={ix1,ix2,ix3,ix4};
2479 void set(
size_t ix1,
size_t ix2,
size_t ix3,
size_t ix4,
2481 size_t sz[4]={ix1,ix2,ix3,ix4};
2489 double interp(
double x,
double y,
double z,
double a) {
2490 double arr[4]={x,y,z,a};
2492 <range_t,data_range_t,index_range_t>(arr);
2497 double arr[4]={x,y,z,a};
2503 #ifndef DOXYGEN_NO_O2NS