46 #ifndef O2SCL_INTERP_H
47 #define O2SCL_INTERP_H
56 #include <boost/numeric/ublas/vector.hpp>
57 #include <boost/numeric/ublas/vector_proxy.hpp>
59 #include <o2scl/search_vec.h>
60 #include <o2scl/tridiag.h>
61 #include <o2scl/vector.h>
63 #ifndef DOXYGEN_NO_O2NS
109 #ifdef O2SCL_NEVER_DEFINED
113 #ifndef DOXYGEN_INTERNAL
136 double integ_eval(
double ai,
double bi,
double ci,
double di,
double xi,
137 double a,
double b)
const {
142 double bterm=0.5*bi*r12;
143 double cterm=ci*(r1*r1+r2*r2+r1*r2)/3.0;
144 double dterm=0.25*di*r12*(r1*r1+r2*r2);
146 return (b-a)*(ai+bterm+cterm+dterm);
157 virtual ~interp_base() {
168 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y)=0;
171 virtual double eval(
double x0)
const=0;
179 virtual double deriv(
double x0)
const=0;
184 virtual double deriv2(
double x0)
const=0;
187 virtual double integ(
double a,
double b)
const=0;
190 virtual const char *type()
const=0;
192 #ifndef DOXYGEN_INTERNAL
213 #ifdef O2SCL_NEVER_DEFINED
223 virtual ~interp_linear() {}
226 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
227 if (size<this->min_size) {
229 " than "+
szttos(this->min_size)+
" in interp_linear::"+
232 this->svx.set_vec(size,x);
240 virtual double eval(
double x0)
const {
243 size_t index=this->svx.find_const(x0,cache);
245 double x_lo=(*this->px)[index];
246 double x_hi=(*this->px)[index+1];
247 double y_lo=(*this->py)[index];
248 double y_hi=(*this->py)[index+1];
251 return y_lo+(x0-x_lo)/dx*(y_hi-y_lo);
255 virtual double deriv(
double x0)
const {
258 size_t index=this->svx.find_const(x0,cache);
260 double x_lo=(*this->px)[index];
261 double x_hi=(*this->px)[index+1];
262 double y_lo=(*this->py)[index];
263 double y_hi=(*this->py)[index+1];
278 virtual double integ(
double a,
double b)
const {
281 size_t i, index_a, index_b;
284 if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
285 ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
292 index_a=this->svx.find_const(a,cache);
293 index_b=this->svx.find_const(b,cache);
296 for(i=index_a; i<=index_b; i++) {
298 double x_lo=(*this->px)[i];
299 double x_hi=(*this->px)[i+1];
300 double y_lo=(*this->py)[i];
301 double y_hi=(*this->py)[i+1];
306 if (i == index_a || i == index_b) {
307 double x1=(i == index_a) ? a : x_lo;
308 double x2=(i == index_b) ? b : x_hi;
309 double D=(y_hi-y_lo)/dx;
310 result += (
x2-
x1)*(y_lo+0.5*D*((
x2-x_lo)+(
x1-x_lo)));
312 result += 0.5*dx*(y_lo+y_hi);
316 std::string str=((std::string)
"Interval of length zero ")+
319 " in interp_linear::integ().";
324 if (flip) result=-result;
329 virtual const char *
type()
const {
return "interp_linear"; }
331 #ifndef DOXYGEN_INTERNAL
350 #ifdef O2SCL_NEVER_DEFINED
360 virtual ~interp_nearest_neigh() {}
363 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
364 if (size<this->min_size) {
366 " than "+
szttos(this->min_size)+
" in interp_nearest_neigh::"+
369 this->svx.set_vec(size,x);
377 virtual double eval(
double x0)
const {
379 size_t index=this->svx.ordered_lookup(x0);
380 return (*this->py)[index];
384 virtual double deriv(
double x0)
const {
396 virtual double integ(
double a,
double b)
const {
401 virtual const char *
type()
const {
return "interp_nearest_neigh"; }
403 #ifndef DOXYGEN_INTERNAL
428 #ifdef O2SCL_NEVER_DEFINED
435 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
436 typedef boost::numeric::ublas::vector_range<ubvector>
ubvector_range;
437 typedef boost::numeric::ublas::slice slice;
438 typedef boost::numeric::ublas::range range;
440 #ifndef DOXYGEN_INTERNAL
457 size_t index,
double &b,
double &c2,
double &d)
const {
459 double c_i=c_array[index];
460 double c_ip1=c_array[index+1];
461 b=(dy/dx)-dx*(c_ip1+2.0*c_i)/3.0;
463 d=(c_ip1-c_i)/(3.0*dx);
484 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
486 if (size<this->min_size) {
488 " than "+
szttos(this->min_size)+
" in interp_cspline::"+
492 if (size!=this->sz) {
496 offdiag.resize(size);
504 this->svx.set_vec(size,xa);
509 size_t num_points=size;
510 size_t max_index=num_points-1;
511 size_t sys_size=max_index-1;
516 for (i=0; i < sys_size; i++) {
517 double h_i=xa[i+1]-xa[i];
518 double h_ip1=xa[i+2]-xa[i+1];
519 double ydiff_i=ya[i+1]-ya[i];
520 double ydiff_ip1=ya[i+2]-ya[i+1];
521 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
522 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
524 diag[i]=2.0*(h_ip1+h_i);
525 g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
535 ubvector_range cp1(c,range(1,c.size()));
538 (diag,offdiag,g,cp1,sys_size,p4m);
544 virtual double eval(
double x0)
const {
547 size_t index=this->svx.find_const(x0,cache);
549 double x_lo=(*this->px)[index];
550 double x_hi=(*this->px)[index+1];
553 double y_lo=(*this->py)[index];
554 double y_hi=(*this->py)[index+1];
557 double b_i, c_i, d_i;
559 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
561 return y_lo+delx*(b_i+delx*(c_i+delx*d_i));
565 virtual double deriv(
double x0)
const {
568 size_t index=this->svx.find_const(x0,cache);
570 double x_lo=(*this->px)[index];
571 double x_hi=(*this->px)[index+1];
574 double y_lo=(*this->py)[index];
575 double y_hi=(*this->py)[index+1];
578 double b_i, c_i, d_i;
580 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
582 return b_i+delx*(2.0*c_i+3.0*d_i*delx);
591 size_t index=this->svx.find_const(x0,cache);
593 double x_lo=(*this->px)[index];
594 double x_hi=(*this->px)[index+1];
597 double y_lo=(*this->py)[index];
598 double y_hi=(*this->py)[index+1];
601 double b_i, c_i, d_i;
603 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
605 return 2.0*c_i+6.0*d_i*delx;
609 virtual double integ(
double a,
double b)
const {
611 size_t i, index_a, index_b;
614 if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
615 ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
623 index_a=this->svx.find_const(a,cache);
624 index_b=this->svx.find_const(b,cache);
628 for(i=index_a; i<=index_b; i++) {
630 double x_lo=(*this->px)[i];
631 double x_hi=(*this->px)[i+1];
632 double y_lo=(*this->py)[i];
633 double y_hi=(*this->py)[i+1];
638 double b_i, c_i, d_i;
639 coeff_calc(c,dy,dx,i,b_i,c_i,d_i);
640 if (i == index_a || i == index_b) {
641 double x1=(i == index_a) ? a : x_lo;
642 double x2=(i == index_b) ? b : x_hi;
643 result += this->integ_eval(y_lo,b_i,c_i,d_i,x_lo,
x1,
x2);
645 result += dx*(y_lo+dx*(0.5*b_i+
646 dx*(c_i/3.0+0.25*d_i*dx)));
649 std::string str=((std::string)
"Interval of length zero ")+
652 " in interp_cspline::integ().";
658 if (flip) result*=-1.0;
664 virtual const char *
type()
const {
return "interp_cspline"; }
666 #ifndef DOXYGEN_INTERNAL
686 #ifdef O2SCL_NEVER_DEFINED
698 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
699 typedef boost::numeric::ublas::vector_range<ubvector>
ubvector_range;
700 typedef boost::numeric::ublas::slice slice;
701 typedef boost::numeric::ublas::range range;
711 virtual const char *
type()
const {
return "interp_cspline_peri"; }
715 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
717 if (size<this->min_size) {
719 " than "+
szttos(this->min_size)+
" in interp_cspline"+
723 if (size!=this->sz) {
724 this->c.resize(size);
725 this->g.resize(size);
726 this->diag.resize(size);
727 this->offdiag.resize(size);
735 this->svx.set_vec(size,xa);
740 size_t num_points=size;
742 size_t max_index=num_points-1;
744 size_t sys_size=max_index;
750 double h0=xa[1]-xa[0];
751 double h1=xa[2]-xa[1];
752 double h2=xa[3]-xa[2];
753 double A=2.0*(h0+h1);
758 gx[0]=3.0*((ya[2]-ya[1])/h1-(ya[1]-ya[0])/h0);
759 gx[1]=3.0*((ya[1]-ya[2])/h2-(ya[2]-ya[1])/h1);
761 det=3.0*(h0+h1)*(h0+h1);
762 this->c[1]=( A*gx[0]-B*gx[1])/det;
763 this->c[2]=(-B*gx[0]+A*gx[1])/det;
764 this->c[0]=this->c[2];
770 for (i=0; i < sys_size-1; i++) {
771 double h_i=xa[i+1]-xa[i];
772 double h_ip1=xa[i+2]-xa[i+1];
773 double ydiff_i=ya[i+1]-ya[i];
774 double ydiff_ip1=ya[i+2]-ya[i+1];
775 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
776 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
777 this->offdiag[i]=h_ip1;
778 this->diag[i]=2.0*(h_ip1+h_i);
779 this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
784 double h_i=xa[i+1]-xa[i];
785 double h_ip1=xa[1]-xa[0];
786 double ydiff_i=ya[i+1]-ya[i];
787 double ydiff_ip1=ya[1]-ya[0];
788 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
789 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
790 this->offdiag[i]=h_ip1;
791 this->diag[i]=2.0*(h_ip1+h_i);
792 this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
795 ubvector_range cp1(this->c,range(1,this->c.size()));
798 (this->diag,this->offdiag,this->g,cp1,sys_size,p5m);
799 this->c[0]=this->c[max_index];
806 #ifndef DOXYGEN_INTERNAL
835 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
836 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
837 typedef boost::numeric::ublas::slice slice;
838 typedef boost::numeric::ublas::range range;
840 #ifndef DOXYGEN_INTERNAL
856 for(
size_t i=0;i<this->
sz-1;i++) {
858 double NE=fabs(umx[3+i]-umx[2+i])+fabs(umx[1+i]-umx[i]);
865 double h_i=(*this->
px)[i+1]-(*this->
px)[i];
866 double NE_next=fabs(umx[4+i]-umx[3+i])+
867 fabs(umx[2+i]-umx[1+i]);
868 double alpha_i=fabs(umx[1+i]-umx[i])/NE;
871 if (NE_next == 0.0) {
874 alpha_ip1=fabs(umx[2+i]-umx[1+i])/NE_next;
875 tL_ip1=(1.0-alpha_ip1)*umx[2+i]+alpha_ip1*umx[3+i];
877 b[i]=(1.0-alpha_i)*umx[1+i]+alpha_i*umx[2+i];
878 c[i]=(3.0*umx[2+i]-2.0*b[i]-tL_ip1)/h_i;
879 d[i]=(b[i]+tL_ip1-2.0*umx[2+i])/(h_i*h_i);
900 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
904 " than "+
szttos(this->min_size)+
" in interp_akima::"+
908 if (size!=this->
sz) {
923 ubvector_range m(um,range(2,um.size()));
925 for (i=0;i<=size-2;i++) {
926 m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
929 um[0]=3.0*m[0]-2.0*m[1];
931 m[this->
sz-1]=2.0*m[size-2]-m[size-3];
932 m[size]=3.0*m[size-2]-2.0*m[size-3];
940 virtual double eval(
double x0)
const {
943 size_t index=this->
svx.find_const(x0,cache);
945 double x_lo=(*this->
px)[index];
951 return (*this->
py)[index]+delx*(bb+delx*(cc+dd*delx));
955 virtual double deriv(
double x0)
const {
958 size_t index=this->
svx.find_const(x0,cache);
960 double x_lo=(*this->
px)[index];
966 return bb+delx*(2.0*cc+3.0*dd*delx);
975 size_t index=this->
svx.find_const(x0,cache);
977 double x_lo=(*this->
px)[index];
982 return 2.0*cc+6.0*dd*delx;
986 virtual double integ(
double aa,
double bb)
const {
988 size_t i, index_a, index_b;
991 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && aa>bb) ||
992 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && aa<bb)) {
1000 index_a=this->
svx.find_const(aa,cache);
1001 index_b=this->
svx.find_const(bb,cache);
1005 for(i=index_a; i<=index_b; i++) {
1007 double x_lo=(*this->
px)[i];
1008 double x_hi=(*this->
px)[i+1];
1009 double y_lo=(*this->
py)[i];
1010 double dx=x_hi-x_lo;
1014 if (i==index_a || i==index_b) {
1015 double x1=(i==index_a) ? aa : x_lo;
1016 double x2=(i==index_b) ? bb : x_hi;
1019 result+=dx*(y_lo+dx*(0.5*b[i]+dx*(c[i]/3.0+0.25*d[i]*dx)));
1023 double y_hi=(*this->
py)[i+1];
1024 std::string str=((std::string)
"Interval of length zero ")+
1027 " in interp_akima::integ().";
1032 if (flip) result*=-1.0;
1038 virtual const char *
type()
const {
return "interp_akima"; }
1040 #ifndef DOXYGEN_INTERNAL
1062 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1063 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
1064 typedef boost::numeric::ublas::slice slice;
1065 typedef boost::numeric::ublas::range range;
1076 virtual const char *
type()
const {
return "interp_akima_peri"; }
1079 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
1083 " than "+
szttos(this->min_size)+
" in interp_akima"+
1087 if (size!=this->
sz) {
1088 this->b.resize(size);
1089 this->c.resize(size);
1090 this->d.resize(size);
1091 this->um.resize(size+4);
1102 ubvector_range m(this->um,range(2,this->um.size()));
1105 for (
size_t i=0;i<=size-2;i++) {
1106 m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
1109 this->um[0]=m[this->
sz-3];
1110 this->um[1]=m[this->
sz-2];
1119 #ifndef DOXYGEN_INTERNAL
1139 #ifdef O2SCL_NEVER_DEFINED
1146 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1147 typedef boost::numeric::ublas::vector_range<ubvector>
ubvector_range;
1148 typedef boost::numeric::ublas::slice slice;
1149 typedef boost::numeric::ublas::range range;
1151 #ifndef DOXYGEN_INTERNAL
1167 if ((x < 0 && y > 0) || (x > 0 && y < 0)) {
1187 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
1189 if (size<this->min_size) {
1191 " than "+
szttos(this->min_size)+
" in interp_steffen::"+
1195 if (size!=this->sz) {
1200 y_prime.resize(size);
1207 this->svx.set_vec(size,xa);
1214 double h0=(xa[1]-xa[0]);
1215 double s0=(ya[1]-ya[0]) / h0;
1221 for (
size_t i=1; i < (size-1); i++) {
1226 double hi=(xa[i+1]-xa[i]);
1227 double him1=(xa[i]-xa[i-1]);
1230 double si=(ya[i+1]-ya[i]) / hi;
1231 double sim1=(ya[i]-ya[i-1]) / him1;
1234 pi=(sim1*hi + si*him1) / (him1 + hi);
1237 y_prime[i]=(copysign(1.0,sim1)+copysign(1.0,si))*
1238 std::min(fabs(sim1),std::min(fabs(si),0.5*fabs(
pi)));
1249 y_prime[size-1]=(ya[size-1]-ya[size-2])/
1250 (xa[size-1]-xa[size-2]);
1253 for (
size_t i=0; i < (size-1); i++) {
1254 double hi=(xa[i+1]-xa[i]);
1255 double si=(ya[i+1]-ya[i]) / hi;
1258 a[i]=(y_prime[i] + y_prime[i+1]-2*si) / hi / hi;
1259 b[i]=(3*si-2*y_prime[i]-y_prime[i+1]) / hi;
1268 virtual double eval(
double x0)
const {
1271 size_t index=this->svx.find_const(x0,cache);
1272 double x_lo=(*this->px)[index];
1273 double delx=x0-x_lo;
1276 double y = d[index]+delx*(c[index]+delx*(b[index]+delx*a[index]));
1285 size_t index=this->svx.find_const(x0,cache);
1286 double x_lo=(*this->px)[index];
1287 double delx=x0-x_lo;
1289 return c[index]+delx*(2.0*b[index]+delx*3.0*a[index]);
1298 size_t index=this->svx.find_const(x0,cache);
1299 double x_lo=(*this->px)[index];
1300 double delx=x0-x_lo;
1302 return 2.0*b[index]+delx*6.0*a[index];
1306 virtual double integ(
double al,
double bl)
const {
1308 size_t i, index_a, index_b;
1311 if (((*this->px)[0]<(*this->px)[this->sz-1] && al>bl) ||
1312 ((*this->px)[0]>(*this->px)[this->sz-1] && al<bl)) {
1320 index_a=this->svx.find_const(al,cache);
1321 index_b=this->svx.find_const(bl,cache);
1325 for(i=index_a; i<=index_b; i++) {
1327 double x_lo=(*this->px)[i];
1328 double x_hi=(*this->px)[i+1];
1329 double y_lo=(*this->py)[i];
1330 double y_hi=(*this->py)[i+1];
1331 double dx=x_hi-x_lo;
1335 double x1=(i == index_a) ? al-x_lo : 0.0;
1336 double x2=(i == index_b) ? bl-x_lo : x_hi-x_lo;
1342 std::string str=((std::string)
"Interval of length zero ")+
1345 " in interp_steffen::integ().";
1351 if (flip) result*=-1.0;
1357 virtual const char *
type()
const {
return "interp_steffen"; }
1359 #ifndef DOXYGEN_INTERNAL
1398 #ifdef O2SCL_NEVER_DEFINED
1427 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
1430 if (size<this->min_size) {
1432 " than "+
szttos(this->min_size)+
" in interp_monotonic::"+
1437 this->svx.set_vec(size,x);
1440 if (this->sz!=size) {
1442 Delta.resize(size-1);
1443 alpha.resize(size-1);
1444 beta.resize(size-1);
1455 for(
size_t i=0;i<size-1;i++) {
1456 Delta[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
1458 m[i]=(Delta[i]+Delta[i-1])/2.0;
1462 m[size-1]=Delta[size-2];
1465 for(
size_t i=0;i<size-1;i++) {
1473 for(
size_t i=0;i<size-1;i++) {
1474 alpha[i]=m[i]/Delta[i];
1475 beta[i]=m[i+1]/Delta[i];
1479 for(
size_t i=0;i<size-1;i++) {
1480 double norm2=alpha[i]*alpha[i]+beta[i]*beta[i];
1482 double tau=3.0/sqrt(norm2);
1483 m[i]=tau*alpha[i]*Delta[i];
1484 m[i+1]=tau*beta[i]*Delta[i];
1492 virtual double eval(
double x0)
const {
1495 size_t index=this->svx.find_const(x0,cache);
1497 double x_lo=(*this->px)[index];
1498 double x_hi=(*this->px)[index+1];
1499 double y_lo=(*this->py)[index];
1500 double y_hi=(*this->py)[index+1];
1502 double t=(x0-x_lo)/h;
1503 double t2=t*t, t3=t2*t;
1505 double h00=2.0*t3-3.0*t2+1.0;
1506 double h10=t3-2.0*t2+t;
1507 double h01=-2.0*t3+3.0*t2;
1509 double interp=y_lo*h00+h*m[index]*h10+y_hi*h01+h*m[index+1]*h11;
1518 size_t index=this->svx.find_const(x0,cache);
1520 double x_lo=(*this->px)[index];
1521 double x_hi=(*this->px)[index+1];
1522 double y_lo=(*this->py)[index];
1523 double y_hi=(*this->py)[index+1];
1525 double t=(x0-x_lo)/h;
1528 double dh00=6.0*t2-6.0*t;
1529 double dh10=3.0*t2-4.0*t+1.0;
1530 double dh01=-6.0*t2+6.0*t;
1531 double dh11=3.0*t2-2.0*t;
1532 double deriv=(y_lo*dh00+h*m[index]*dh10+y_hi*dh01+
1533 h*m[index+1]*dh11)/h;
1544 size_t index=this->svx.find_const(x0,cache);
1546 double x_lo=(*this->px)[index];
1547 double x_hi=(*this->px)[index+1];
1548 double y_lo=(*this->py)[index];
1549 double y_hi=(*this->py)[index+1];
1551 double t=(x0-x_lo)/h;
1553 double ddh00=12.0*t-6.0;
1554 double ddh10=6.0*t-4.0;
1555 double ddh01=-12.0*t+6.0;
1556 double ddh11=6.0*t-2.0;
1557 double deriv2=(y_lo*ddh00+h*m[index]*ddh10+y_hi*ddh01+
1558 h*m[index+1]*ddh11)/h/h;
1564 virtual double integ(
double a,
double b)
const {
1566 size_t i, index_a, index_b;
1569 if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
1570 ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
1578 index_a=this->svx.find_const(a,cache);
1579 index_b=this->svx.find_const(b,cache);
1583 for(i=index_a; i<=index_b; i++) {
1585 double x_hi=(*this->px)[i+1];
1586 double x_lo=(*this->px)[i];
1587 double y_lo=(*this->py)[i];
1588 double y_hi=(*this->py)[i+1];
1600 double t=(x_hi-x_lo)/h;
1601 double t2=t*t, t3=t2*t, t4=t3*t;
1603 double ih00=t4/2.0-t3+t;
1604 double ih10=t4/4.0-2.0*t3/3.0+t2/2.0;
1605 double ih01=-t4/2.0+t3;
1606 double ih11=t4/4.0-t3/3.0;
1607 double intres=h*(y_lo*ih00+h*m[i]*ih10+y_hi*ih01+
1612 std::string str=((std::string)
"Interval of length zero ")+
1615 " in interp_monotonic::integ().";
1621 if (flip) result*=-1.0;
1627 virtual const char *
type()
const {
return "interp_monotonic"; }
1629 #ifndef DOXYGEN_INTERNAL
1653 template<
class vec_t=boost::numeric::ublas::vector<
double>,
1656 #ifndef DOXYGEN_INTERNAL
1686 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1699 virtual double eval(
const double x0,
size_t n,
const vec_t &x,
1702 return itp->eval(x0);
1710 virtual double deriv(
const double x0,
size_t n,
const vec_t &x,
1713 return itp->deriv(x0);
1720 virtual double deriv2(
const double x0,
size_t n,
const vec_t &x,
1723 return itp->deriv2(x0);
1730 virtual double integ(
const double x1,
const double x2,
size_t n,
1731 const vec_t &x,
const vec2_t &y) {
1756 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1763 #ifndef DOXYGEN_INTERNAL
1794 template<
class vec_t=boost::numeric::ublas::vector<
double>,
1798 #ifndef DOXYGEN_INTERNAL
1823 const vec2_t &y,
size_t interp_type=
itp_cspline) {
1826 O2SCL_ERR((((std::string)
"Vector endpoints equal (value=")+
1848 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1850 "interp_vec::interp_vec().").c_str(),
exc_einval);
1866 void set(
size_t n,
const vec_t &x,
const vec2_t &y) {
1873 void set(
size_t n,
const vec_t &x,
1874 const vec2_t &y,
size_t interp_type) {
1877 O2SCL_ERR((((std::string)
"Vector endpoints equal (value=")+
1900 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1920 virtual double eval(
const double x0)
const {
1922 O2SCL_ERR(
"No vector set in interp_vec::eval().",
1925 return itp->eval(x0);
1931 O2SCL_ERR(
"No vector set in interp_vec::operator().",
1934 return itp->eval(x0);
1938 virtual double deriv(
const double x0)
const {
1940 O2SCL_ERR(
"No vector set in interp_vec::deriv().",
1943 return itp->deriv(x0);
1949 virtual double deriv2(
const double x0)
const {
1951 O2SCL_ERR(
"No vector set in interp_vec::deriv2().",
1954 return itp->deriv2(x0);
1958 virtual double integ(
const double x1,
const double x2)
const {
1960 O2SCL_ERR(
"No vector set in interp_vec::integ().",
1968 return "interp_vec";
1971 #ifndef DOXYGEN_INTERNAL
1988 public interp<double[n]> {
1994 :
interp<double[n]>(interp_type) {}
2014 size_t interp_type) :
2042 (
double level,
size_t n, vec_t &x, vec2_t &y) {
2045 O2SCL_ERR2(
"Need at least two data points in ",
2052 if (y[0]==level) count++;
2055 for(
size_t i=0;i<n-1;i++) {
2057 if ((y[i]<level && y[i+1]>=level) ||
2058 (y[i]>level && y[i+1]<=level)) {
2072 template<
class ovec_t,
class vec2_t>
2076 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2078 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(((
double)i));
2085 template<
class ovec_t,
class vec2_t>
2089 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2091 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv2(((
double)i));
2098 template<
class vec_t,
class vec2_t,
class vec3_t>
2102 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(vx[i]);
2109 template<
class vec_t,
class vec2_t,
class vec3_t>
2113 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(vx[i]);
2119 template<
class ovec_t>
2122 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2124 return oi.
integ(0.0,((
double)(n-1)));
2140 template<
class vec_t,
class vec2_t>
2148 double total=si.
integ(x[0],x[n-1],n,x,y);
2156 template<
class vec_t,
class vec2_t,
class vec3_t>
2165 for(
size_t i=1;i<n;i++) {
2166 iy[i]=si.
integ(x[0],x[i],n,x,y);
2175 template<
class ovec_t>
2177 ovec_t &v,
size_t interp_type) {
2179 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2187 template<
class vec_t,
class vec2_t>
2189 const vec2_t &y,
double x2,
2196 double total=si.
integ(x[0],
x2,n,x,y);
2225 (
double level,
size_t n, vec_t &x, vec2_t &y, std::vector<double> &locs) {
2228 O2SCL_ERR2(
"Need at least two data points in ",
2237 locs.push_back(x[0]);
2241 for(
size_t i=0;i<n-1;i++) {
2243 if ((y[i]<level && y[i+1]>level) ||
2244 (y[i]>level && y[i+1]<level)) {
2247 double x0=x[i]+(x[i+1]-x[i])*(level-y[i])/(y[i+1]-y[i]);
2249 }
else if (y[i+1]==level) {
2250 locs.push_back(x[i+1]);
2295 (
double sum,
size_t n, vec_t &x, vec2_t &y,
double &lev,
2296 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2299 O2SCL_ERR2(
"Need at least two data points in ",
2306 std::vector<double>
x2, y2;
2308 if (boundaries==1) {
2310 std::cout <<
"Fix left boundary to zero." << std::endl;
2314 x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2316 for(
size_t i=0;i<n;i++) {
2321 }
else if (boundaries==2) {
2323 std::cout <<
"Fix right boundary to zero." << std::endl;
2327 for(
size_t i=0;i<n;i++) {
2331 x2[n]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2334 }
else if (boundaries==3) {
2336 std::cout <<
"Fix both boundaries to zero." << std::endl;
2340 x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2342 for(
size_t i=0;i<n;i++) {
2346 x2[n+1]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2351 std::cout <<
"No boundary extrapolation." << std::endl;
2355 for(
size_t i=0;i<n;i++) {
2365 vector_sort<ubvector,double>(ysort.size(),ysort);
2371 std::vector<double> ylist;
2372 for(
size_t i=0;i<ysort.size()-1;i++) {
2373 if (ysort[i]!=ysort[i+1]) {
2374 ylist.push_back((ysort[i+1]+3.0*ysort[i])/4.0);
2375 ylist.push_back((ysort[i+1]*3.0+ysort[i])/4.0);
2385 std::vector<double> xi, yi;
2389 for(
size_t k=0;k<ylist.size();k++) {
2390 double lev_tmp=ylist[k];
2391 std::vector<double> locs;
2393 if ((locs.size()%2)!=0) {
2396 std::cout << k <<
" " << lev_tmp <<
" " << 0.0 <<
" "
2397 << locs.size() <<
" (fail)" << std::endl;
2400 double sum_temp=0.0;
2401 for(
size_t i=0;i<locs.size()/2;i++) {
2402 double x0=locs[2*i];
2403 double x1=locs[2*i+1];
2406 xi.push_back(sum_temp);
2407 yi.push_back(lev_tmp);
2409 std::cout << k <<
" " << lev_tmp <<
" " << sum_temp <<
" "
2410 << locs.size() <<
" ";
2411 for(
size_t i=0;i<locs.size();i++) {
2412 std::cout << locs[i] <<
" ";
2414 std::cout << std::endl;
2426 lev=itp2.
eval(sum,xi.size(),xi,yi);
2434 (
size_t n, vec_t &x, vec2_t &y,
double intl, std::vector<double> &locs,
2435 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2440 std::cout <<
"Total integral: " << total << std::endl;
2444 std::cout <<
"Target integral: " << intl << std::endl;
2449 boundaries,verbose,err_on_fail);
2452 O2SCL_ERR2(
"Failed to find a level which enclosed the ",
2453 "specified integral in vector_region_int().",
2459 std::cout <<
"Level from vector_invert: " << lev << std::endl;
2465 std::cout <<
"Locations from vector_find_level: ";
2466 for(
size_t i=0;i<locs.size();i++) {
2467 std::cout << locs[i];
2468 if (i!=locs.size()-1) {
2472 std::cout << std::endl;
2480 (
size_t n, vec_t &x, vec2_t &y,
double frac, std::vector<double> &locs,
2481 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2484 verbose,err_on_fail);
2494 (
size_t n, vec_t &x, vec2_t &y,
double frac,
double &low,
double &high,
2495 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2497 std::vector<double> locs;
2499 verbose,err_on_fail);
2500 if (locs.size()==0 || ret!=0) {
2502 O2SCL_ERR2(
"Zero level crossings or vector_region_fracint() ",
2521 (
size_t n, vec_t &x, vec2_t &y,
double frac,
double &low,
double &high,
2522 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2524 std::vector<double> locs;
2526 verbose,err_on_fail);
2527 if (locs.size()==0 || ret!=0) {
2529 O2SCL_ERR2(
"Zero level crossings or vector_region_int() ",
2547 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t>
2549 vec3_t &new_x, vec4_t &new_y,
size_t n_pts,
2550 size_t interp_type) {
2552 if (x.size()!=y.size()) {
2553 O2SCL_ERR2(
"The x and y vectors must have the same size ",
2557 O2SCL_ERR2(
"Number of points must be at least 2 ",
2569 new_y.resize(n_pts);
2571 for(
size_t i=0;i<n_pts;i++) {
2572 new_y[i]=itp.
eval(new_x[i]);
2582 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t>
2584 vec3_t &new_x, vec4_t &new_y,
size_t n_pts,
2585 size_t interp_type1,
size_t interp_type2,
2586 double acc=1.0e-4) {
2588 if (x.size()!=y.size()) {
2589 O2SCL_ERR2(
"The x and y vectors must have the same size ",
2593 O2SCL_ERR2(
"Number of points must be at least 2 ",
2605 std::vector<double> new_y2(n_pts);
2606 new_y.resize(n_pts);
2609 for(
size_t i=0;i<n_pts;i++) {
2610 new_y[i]=itp1.
eval(new_x[i]);
2611 new_y2[i]=itp2.
eval(new_x[i]);
2614 double max_y, min_y;
2616 for(
size_t i=0;i<n_pts;i++) {
2617 if (fabs(new_y2[i]-new_y[i])/(max_y-min_y)>acc) {
2629 template<
class vec_t,
class vec2_t>
2632 static const size_t n2=11;
2635 std::vector<double> xn, yn;
2639 double min_y, max_y;
2641 for(
size_t i=0;i<n2;i++) {
2642 yn[i]=(yn[i]-min_y)/(max_y-min_y);
2648 for(
size_t i=0;i<n2;i++) {
2649 double ty=((double)i)/10.0;
2650 chi2+=pow(yn[i]-ty,2.0);
2665 template<
class vec_t,
class vec2_t>
2667 if (x.size()!=y.size()) {
2668 O2SCL_ERR2(
"The x and y vectors must have the same size ",
2675 O2SCL_ERR2(
"Vector size must be at least 2 ",
2684 bool x_positive=
true;
2685 bool y_positive=
true;
2686 for(
size_t i=0;i<n;i++) {
2687 if (x[i]<=0.0) x_positive=
false;
2688 if (y[i]<=0.0) y_positive=
false;
2691 if (x_positive==
false && y_positive==
false)
return;
2696 std::vector<double> lx(n), ly(n);
2698 for(
size_t i=0;i<n;i++) {
2703 for(
size_t i=0;i<n;i++) {
2714 if (chi2_xy<chi2_x && chi2_xy<chi2_y && chi2_xy<chi2) {
2717 }
else if (chi2_x<chi2_xy && chi2_x<chi2_y && chi2_x<chi2) {
2719 }
else if (chi2_y<chi2_xy && chi2_y<chi2_x && chi2_y<chi2) {
2724 if (chi2_x<chi2) log_x=
true;
2728 if (chi2_y<chi2) log_y=
true;
2739 template<
class vec_t,
class vec2_t,
class data_t>
2741 size_t factor,
size_t interp_type=
itp_linear) {
2744 data.resize((n-1)*factor+1);
2745 for (
size_t j=0;j<n-1;j++) {
2746 for(
size_t k=0;k<factor;k++) {
2747 data[j*factor+k]=iv.
eval(index[j]+
2748 ((data_t)k)/((data_t)factor)*
2749 (index[j+1]-index[j]));
2752 data[data.size()-1]=copy[copy.size()-1];
2765 template<
class vec_t>
2767 std::vector<double> x(y.size());
2768 for(
size_t i=0;i<y.size();i++) {
2777 #ifndef DOXYGEN_NO_O2NS