1 #ifndef __JTOOLS__JPOLINT__
2 #define __JTOOLS__JPOLINT__
37 template<
unsigned int N,
39 template<
class,
class>
class JCollection_t,
52 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
58 public JCollection_t<JElement_t, JDistance_t>,
59 public JFunction<typename JElement_t::abscissa_type,
60 typename JResultType<typename JElement_t::ordinate_type>::result_type>
73 typedef typename collection_type::iterator
iterator;
99 if (this->size() <= 1
u) {
100 return this->getExceptionHandler().action(
JFunctionalException(
"JPolintFunction<>::evaluate() not enough data."));
107 if ((p == this->begin() && this->
getDistance(x, (p++)->getX()) > distance_type::precision) ||
108 (p == this->end() && this->
getDistance((--p)->getX(), x) > distance_type::precision)) {
109 return this->getExceptionHandler().action(
JValueOutOfRange(
"JPolintFunction::evaluate() x out of range."));
115 const int n = std::min((
int) (N + 1), (
int) this->size());
117 for (
int i =
n/2; i != 0 && p != this->end(); --i, ++p) {}
118 for (
int i =
n ; i != 0 && p != this->begin(); --i, --p) {}
123 for (
int i = 0; i !=
n; ++p, ++i) {
129 if (fabs(
u[i]) < fabs(
u[
j])) {
139 for (
int m = 1; m !=
n; ++m) {
141 for (
int i = 0; i !=
n-m; ++i) {
143 const double ho =
u[ i ];
144 const double hp =
u[i+m];
145 const double dx = ho - hp;
172 mutable double u[N+1];
182 template<
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
188 public JCollection_t<JElement_t, JDistance_t>,
189 public JFunction<typename JElement_t::abscissa_type,
190 typename JResultType<typename JElement_t::ordinate_type>::result_type>
203 typedef typename collection_type::iterator
iterator;
230 return this->getExceptionHandler().action(
JFunctionalException(
"JPolintFunction<>::evaluate() not enough data."));
237 if ((p == this->begin() && this->
getDistance(x, (p++)->getX()) > distance_type::precision) ||
238 (p == this->end() && this->
getDistance((--p)->getX(), x) > distance_type::precision)) {
240 return this->getExceptionHandler().action(
JValueOutOfRange(
"JPolintFunction::evaluate() x out of range."));
267 template<
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
273 public JCollection_t<JElement_t, JDistance_t>,
274 public JFunction<typename JElement_t::abscissa_type,
275 typename JResultType<typename JElement_t::ordinate_type>::result_type>
288 typedef typename collection_type::iterator
iterator;
314 if (this->size() <= 1
u) {
315 return this->getExceptionHandler().action(
JFunctionalException(
"JPolintFunction<>::evaluate() not enough data."));
322 if ((p == this->begin() && this->
getDistance(x, (p++)->getX()) > distance_type::precision) ||
323 (p == this->end() && this->
getDistance((--p)->getX(), x) > distance_type::precision)) {
325 return this->getExceptionHandler().action(
JValueOutOfRange(
"JPolintFunction::evaluate() x out of range."));
333 const double dx = this->
getDistance(p->getX(), q->getX());
334 const double a = this->
getDistance(x, q->getX()) / dx;
335 const double b = 1.0 - a;
370 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
376 public JCollection_t<JElement_t, JDistance_t>,
377 public JFunction<typename JElement_t::abscissa_type,
378 JResultPDF<typename JResultType<typename JElement_t::ordinate_type>::result_type> >
391 typedef typename collection_type::iterator
iterator;
419 if (this->size() <= 1
u) {
420 return this->getExceptionHandler().action(
JFunctionalException(
"JPolintFunction<>::evaluate() not enough data."));
427 if ((p == this->begin() && this->
getDistance(x, (p++)->getX()) > distance_type::precision) ||
428 (p == this->end() && this->
getDistance((--p)->getX(), x) > distance_type::precision)) {
430 return this->getExceptionHandler().action(
JValueOutOfRange(
"JPolintFunction::evaluate() x out of range."));
436 const int n = std::min((
int) (N + 1), (
int) this->size());
438 for (
int i =
n/2; i != 0 && p != this->end(); --i, ++p) {}
439 for (
int i =
n ; i != 0 && p != this->begin(); --i, --p) {}
444 for (
int i = 0; i !=
n; ++p, ++i) {
450 w[i][2] =
v[i][2] = p->getIntegral();
452 if (fabs(
u[i]) < fabs(
u[
j])) {
461 result.V = this->rbegin()->getIntegral();
465 for (
int m = 1; m !=
n; ++m) {
467 for (
int i = 0; i !=
n-m; ++i) {
469 const double ho =
u[ i ];
470 const double hp =
u[i+m];
471 const double dx = ho - hp;
473 for (
int k = 0; k != 3; ++k) {
474 r[k] = (
v[i+1][k] -
w[i][k]) / dx;
479 v[i][1] = ho *
r[1] -
r[0];
480 w[i][1] = hp *
r[1] -
r[0];
485 if (2*(
j+1) <
n - m) {
517 this->begin()->setIntegral(V);
519 for (
iterator j = this->begin(), i =
j++;
j != this->end(); ++i, ++
j) {
526 const abscissa_type x = 0.5 * (xmax + xmin + m->getX() * (xmax - xmin));
538 mutable double u[N+1];
550 template<
unsigned int N,
552 template<
class,
class>
class JCollection_t,
563 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t,
unsigned int M>
569 public JCollection_t<JElement_t, JDistance_t>,
575 typedef JCollection_t<JElement_t, JDistance_t> collection_type;
577 typedef typename collection_type::abscissa_type abscissa_type;
578 typedef typename collection_type::ordinate_type ordinate_type;
579 typedef typename collection_type::value_type value_type;
580 typedef typename collection_type::distance_type distance_type;
582 typedef typename collection_type::const_iterator const_iterator;
583 typedef typename collection_type::const_reverse_iterator const_reverse_iterator;
584 typedef typename collection_type::iterator iterator;
585 typedef typename collection_type::reverse_iterator reverse_iterator;
587 typedef typename JResultType<ordinate_type>::result_type data_type;
588 typedef JFunction<abscissa_type, JResultPolynome<M, data_type> > function_type;
590 typedef typename function_type::argument_type argument_type;
591 typedef typename function_type::result_type result_type;
593 using JFunctional<>::compile;
609 result_type evaluate(const argument_type* pX) const
611 if (this->size() <= N) {
612 THROW(JFunctionalException, "JPolintFunction<>::evaluate() not enough data.");
615 const argument_type x = *pX;
617 const_iterator p = this->lower_bound(x);
619 if ((p == this->begin() && this->
getDistance(x, (p++)->getX()) > distance_type::precision) ||
620 (p == this->end() && this->
getDistance((--p)->getX(), x) > distance_type::precision)) {
628 const int n = std::min((
int) (N + 1), (
int) this->size());
630 for (
int i =
n/2; i != 0 && p != this->end(); --i, ++p) {}
631 for (
int i =
n ; i != 0 && p != this->begin(); --i, --p) {}
636 for (
int i = 0; i !=
n; ++p, ++i) {
642 for (
unsigned int k = 1; k != M+1; ++k) {
646 if (fabs(
u[i]) < fabs(
u[
j])) {
652 for (
unsigned int k = 0; k != M+1; ++k) {
658 for (
int m = 1; m !=
n; ++m) {
660 for (
int i = 0; i !=
n-m; ++i) {
662 const double ho =
u[ i ];
663 const double hp =
u[i+m];
664 const double dx = ho - hp;
666 for (
int k = 0; k != M+1; ++k) {
667 r[k] = (
v[i+1][k] -
w[i][k]) / dx;
673 for (
int k = 1; k != M+1; ++k) {
674 v[i][k] = ho *
r[k] - k *
r[k-1];
675 w[i][k] = hp *
r[k] - k *
r[k-1];
679 if (2*(
j+1) <
n - m) {
681 for (
int k = 0; k != M+1; ++k) {
687 for (
int k = 0; k != M+1; ++k) {
706 mutable double u[N+1];
707 mutable data_type
v[N+1][M+1];
708 mutable data_type
w[N+1][M+1];
709 mutable data_type
r[M+1];
711 mutable result_type
result;
718 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t,
unsigned int M>
727 JResultPolynome<M, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
729 public JFunction<typename JElement_t::abscissa_type,
730 JResultPolynome<N, typename JResultType<typename JElement_t::ordinate_type>::result_type> >
747 typedef typename collection_type::iterator
iterator;
774 return collection_type::evaluate(
pX);
777 return this->getExceptionHandler().action(error);
786 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
795 JResultPolynome<1, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
797 public JFunction<typename JElement_t::abscissa_type,
798 JResultDerivative<typename JResultType<typename JElement_t::ordinate_type>::result_type> >
815 typedef typename collection_type::iterator
iterator;
844 return collection_type::evaluate(
pX);
847 return this->getExceptionHandler().action(error);
856 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
865 JResultPolynome<2, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
867 public JFunction<typename JElement_t::abscissa_type,
868 JResultHesse<typename JResultType<typename JElement_t::ordinate_type>::result_type> >
885 typedef typename collection_type::iterator
iterator;
914 return collection_type::evaluate(
pX);
917 return this->getExceptionHandler().action(error);
928 template<
unsigned int N,
930 template<
class,
class>
class JCollection_t,
931 class JResult_t =
typename JElement_t::ordinate_type,
934 public JPolintFunction<N, JElement_t, JCollection_t, JResult_t, JDistance_t>,
935 public JFunction1D<typename JElement_t::abscissa_type, JResult_t>
969 template<
unsigned int N,
972 template<
class,
class,
class>
class JMap_t,
977 JElement2D<JKey_t, JValue_t>,
978 JMapCollection<JMap_t>::template collection_type,
998 typedef typename JPolintFunction_t::iterator
iterator;
1023 template<
unsigned int N,
1025 template<
class,
class>
class JCollection_t,
1028 inline typename JElement_t::ordinate_type
1047 template<
unsigned int N,
1049 template<
class,
class>
class JCollection_t,
1052 inline typename JElement_t::ordinate_type
1057 typedef typename JElement_t::abscissa_type abscissa_type;
1058 typedef typename JElement_t::ordinate_type ordinate_type;
1063 if (input.getSize() > 1) {
1065 output.
put(input.begin()->getX(), V);
1069 for (const_iterator
j = input.begin(), i =
j++;
j != input.end(); ++i, ++
j) {
1071 const abscissa_type xmin = i->getX();
1072 const abscissa_type xmax =
j->getX();
1076 const abscissa_type x = 0.5 * (xmax + xmin + m->getX() * (xmax - xmin));
1077 const ordinate_type
v = 0.5 * (xmax - xmin) * m->getY() *
get_value(input(x));
1082 output.
put(xmax, V);
1100 template<
class JElement_t,
1101 template<
class,
class>
class JCollection_t,
1104 inline typename JElement_t::ordinate_type
1109 typedef typename JElement_t::ordinate_type ordinate_type;
1114 if (input.getSize() > 1) {
1116 output.
put(input.begin()->getX(), V);
1118 for (const_iterator
j = input.begin(), i =
j++;
j != input.end(); ++i, ++
j) {
1120 V += input.getDistance(i->getX(),
j->getX()) *
j->getY();
1122 output.
put(
j->getX(), V);
1140 template<
class JElement_t,
1141 template<
class,
class>
class JCollection_t,
1144 inline typename JElement_t::ordinate_type
1149 typedef typename JElement_t::ordinate_type ordinate_type;
1154 if (input.getSize() > 1) {
1156 output.
put(input.begin()->getX(), V);
1158 for (const_iterator
j = input.begin(), i =
j++;
j != input.end(); ++i, ++
j) {
1160 V += 0.5 * input.getDistance(i->getX(),
j->getX()) * (i->getY() +
j->getY());
1162 output.
put(
j->getX(), V);