1#ifndef __JTOOLS__JSPLINE__
2#define __JTOOLS__JSPLINE__
36 template<
class JOrdinate_t>
122 throw JNoValue(
"JSplineBounds: missing 1st derivative.");
136 throw JNoValue(
"JSplineBounds: missing 1st derivative.");
152 template<
class JOrdinate_t>
154 const JOrdinate_t fpAtXmax)
176 template<
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
178 public JCollection_t<JElement_t, JDistance_t>,
191 typedef typename collection_type::iterator
iterator;
204 const int numberOfElements = this->size();
208 if (numberOfElements > 2) {
217 const double dx = this->getDistance(i->getX(),
j->getX());
233 for (
iterator k = this->begin(), i = k++,
j = k++; k != this->end(); ++i, ++
j, ++k, ++index) {
243 const double sig = this->getDistance(x1, x2) / this->getDistance(x1, x3);
244 const double h = sig * buffer[index-1] + 2.0;
246 buffer[index] = (sig - 1.0) / h;
248 j->setU((y3 - y2) / this->getDistance(x2, x3) -
249 (y2 - y1) / this->getDistance(x1, x2));
251 j->setU((6.0 *
j->getU() / this->getDistance(x1, x3) - sig * i->getU()) / h);
260 index = numberOfElements - 2;
262 const double dx = this->getDistance(i->getX(),
j->getX());
267 i->setU((i->getU() - 0.5*
j->getU()) / (0.5*buffer[index] + 1.0));
278 index = numberOfElements - 2;
280 for ( ;
j != this->rend(); ++i, ++
j, --index) {
281 j->setU(
j->getU() + i->getU() * buffer[index]);
286 for (
iterator i = this->begin(); i != this->end(); ++i) {
314 template<
class JElement_t,
315 template<
class,
class>
class JCollection_t,
324 template<
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
330 public virtual JFunctional<typename JElement_t::abscissa_type,
331 typename JResultType<typename JElement_t::ordinate_type>::result_type>
372 if (this->size() > 1u) {
376 if ((p == this->begin() && this->getDistance(x, (p++)->getX()) > distance_type::precision) ||
377 (p == this->end() && this->getDistance((--p)->getX(), x) > distance_type::precision)) {
380 return this->getExceptionHandler().action();
384 std::ostringstream os;
386 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range "
387 <<
STREAM(
"?") << x <<
" <> "
388 <<
STREAM(
"?") << this->begin() ->getX() <<
' '
389 <<
STREAM(
"?") << this->rbegin()->getX();
397 const double dx = this->getDistance(p->getX(), q->getX());
398 const double a = this->getDistance(x, q->getX()) / dx;
399 const double b = 1.0 - a;
401 return (a * p->getY() + b * q->getY()
402 - a*b * ((a + 1.0)*p->getU() + (b + 1.0)*q->getU()) * dx*dx/6);
404 }
else if (this->size() == 1u && this->getDistance(x, this->begin()->getX()) <= distance_type::precision) {
406 return this->begin()->getY();
411 return this->getExceptionHandler().action();
415 std::ostringstream os;
417 os << __FILE__ <<
':' << __LINE__ <<
" not enough data " <<
STREAM(
"?") << x;
429 template<
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
435 public virtual JFunctional<typename JElement_t::abscissa_type,
436 JResultDerivative<typename JResultType<typename JElement_t::ordinate_type>::result_type> >
481 if (this->size() <= 1u) {
484 return this->getExceptionHandler().action();
488 std::ostringstream os;
490 os << __FILE__ <<
':' << __LINE__ <<
" not enough data " <<
STREAM(
"?") << x;
499 if ((p == this->begin() && this->getDistance(x, (p++)->getX()) > distance_type::precision) ||
500 (p == this->end() && this->getDistance((--p)->getX(), x) > distance_type::precision)) {
503 return this->getExceptionHandler().action();
507 std::ostringstream os;
509 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range "
510 <<
STREAM(
"?") << x <<
" <> "
511 <<
STREAM(
"?") << this->begin() ->getX() <<
' '
512 <<
STREAM(
"?") << this->rbegin()->getX();
520 const double dx = this->getDistance(p->getX(), q->getX());
521 const double a = this->getDistance(x, q->getX()) / dx;
522 const double b = 1.0 - a;
524 result.f = a * p->getY() + b * q->getY()
525 - a*b * ((a + 1.0)*p->getU() + (b + 1.0)*q->getU()) * dx*dx/6;
527 result.fp = (q->getY() - p->getY() + (p->getU()*(1.0 - 3*a*a) -
528 q->getU()*(1.0 - 3*b*b)) * dx*dx/6) / dx;
545 template<
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
551 public virtual JFunctional<typename JElement_t::abscissa_type,
552 JResultPDF<typename JResultType<typename JElement_t::ordinate_type>::result_type> >
592 if (this->size() >= 2u) {
594 collection_type::compile(bounds);
598 for (
iterator j = this->begin(), i =
j++;
j != this->end(); ++i, ++
j) {
600 const double dx = this->getDistance(i->getX(),
j->getX());
607 j->setIntegral(i->getIntegral() + v - w);
625 if (this->size() <= 1u) {
628 return this->getExceptionHandler().action();
632 std::ostringstream os;
634 os << __FILE__ <<
':' << __LINE__ <<
" not enough data " <<
STREAM(
"?") << x;
642 if (p == this->begin() && this->getDistance(x, (p++)->getX()) > distance_type::precision) {
646 result = this->getExceptionHandler().action();
651 result.V = this->rbegin()->getIntegral();
655 std::ostringstream os;
657 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range " <<
STREAM(
"?") << x <<
" < " <<
STREAM(
"?") << this->begin() ->getX();
664 }
else if (p == this->end() && this->getDistance((--p)->getX(), x) > distance_type::precision) {
668 result = this->getExceptionHandler().action();
672 result.v = this->rbegin()->getIntegral();
673 result.V = this->rbegin()->getIntegral();
677 std::ostringstream os;
679 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range " <<
STREAM(
"?") << x <<
" > " <<
STREAM(
"?") << this->rbegin()->getX();
689 const double dx = this->getDistance(p->getX(), q->getX());
690 const double a = this->getDistance(x, q->getX()) / dx;
691 const double b = 1.0 - a;
693 result.f = a * p->getY() + b * q->getY()
694 - a*b * ((a + 1.0)*p->getU() + (b + 1.0)*q->getU()) * dx*dx/6;
696 result.fp = (q->getY() - p->getY() + (p->getU()*(1.0 - 3*a*a) -
697 q->getU()*(1.0 - 3*b*b)) * dx*dx/6) / dx;
699 result.v = p->getIntegral()
700 + 0.5*dx * (p->getY() - 0.5*p->getU()*dx*dx/6)
701 - 0.5*dx * ((a*a*p->getY() - b*b*q->getY()) +
702 (p->getU() * a*a*(0.5*a*a - 1.0) -
703 q->getU() * b*b*(0.5*b*b - 1.0)) * dx*dx/6);
705 result.V = this->rbegin()->getIntegral();
727 template<
class JElement_t,
728 template<
class,
class>
class JCollection_t,
729 class JResult_t =
typename JElement_t::ordinate_type,
732 public JSplineFunction<JElement_t, JCollection_t, JResult_t, JDistance_t>,
733 public JFunction1D<typename JElement_t::abscissa_type, JResult_t>
769 template<
class JAbscissa_t,
class JOrdinate_t>
772 template<
template<
class,
class,
class>
class JMap_t>
779 template<
class JKey_t,
781 template<
class,
class,
class>
class JMap_t,
785 public JMap_t<JKey_t, JValue_t, JDistance_t>,
800 typedef typename collection_type::iterator
iterator;
852 for (
iterator i = this->begin(); i != this->end(); ++i) {
871 template<
class JElement_t,
872 template<
class,
class>
class JCollection_t,
875 inline typename JElement_t::ordinate_type
879 typedef typename JElement_t::ordinate_type ordinate_type;
884 if (input.getSize() > 1) {
886 output.
put(input.begin()->getX(), V);
888 for (const_iterator
j = input.begin(), i =
j++;
j != input.end(); ++i, ++
j) {
890 const double dx = input.getDistance(i->getX(),
j->getX());
891 const ordinate_type y = i->getY() +
j->getY();
892 const ordinate_type z = i->getU() +
j->getU();
893 const ordinate_type v = dx * 0.50 * y;
894 const ordinate_type w = dx * 0.25 * z*dx*dx/6;
898 output.
put(
j->getX(), V);
915 template<
class JElement_t,
916 template<
class,
class>
class JCollection_t,
918 inline typename JElement_t::ordinate_type
922 typedef typename JElement_t::ordinate_type ordinate_type;
926 if (input.getSize() > 1) {
928 for (const_iterator i = input.begin(); i != input.end(); ++i) {
929 output.
put(i->getX(), i->getIntegral());
932 return input.rbegin()->getIntegral();
This include file containes various data structures that can be used as specific return types for the...
Definition of zero value for any class.
Exception for division by zero.
Exception for a functional operation.
Exception for missing value.
Exception for accessing a value in a collection that is outside of its range.
static const JZero zero
Function object to assign zero value.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
JArgument< T >::argument_type argument_type
Auxiliary data structure for handling std::ostream.