1 #ifndef __JTOOLS__JPOLINT__
2 #define __JTOOLS__JPOLINT__
26 namespace JPP {
using namespace JTOOLS; }
39 template<
unsigned int N,
41 template<
class,
class>
class JCollection_t,
54 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
60 public JCollection_t<JElement_t, JDistance_t>,
61 public virtual JFunction<typename JElement_t::abscissa_type,
62 typename JResultType<typename JElement_t::ordinate_type>::result_type>
75 typedef typename collection_type::iterator
iterator;
103 if (this->size() > 1
u) {
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)) {
111 return this->getExceptionHandler().action();
115 std::ostringstream os;
117 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range "
118 <<
STREAM(
"?") << x <<
" <> "
119 <<
STREAM(
"?") << this->begin() ->getX() <<
' '
120 <<
STREAM(
"?") << this->rbegin()->getX();
129 const int n = std::min((
int) (
N + 1), (
int) this->size());
131 for (
int i = n/2;
i != 0 && p != this->end(); --
i, ++p) {}
132 for (
int i = n ;
i != 0 && p != this->begin(); --
i, --p) {}
137 for (
int i = 0;
i !=
n; ++p, ++
i) {
143 if (fabs(
u[
i]) < fabs(
u[j])) {
153 for (
int m = 1; m !=
n; ++m) {
155 for (
int i = 0;
i != n-m; ++
i) {
157 const double ho =
u[
i ];
158 const double hp =
u[
i+m];
159 const double dx = ho - hp;
177 }
else if (this->size() == 1
u && this->
getDistance(x, this->begin()->getX()) <= distance_type::precision) {
184 return this->getExceptionHandler().action();
188 std::ostringstream os;
190 os << __FILE__ <<
':' << __LINE__ <<
" not enough data " <<
STREAM(
"?") <<
x;
206 mutable double u[
N+1];
216 template<
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
222 public JCollection_t<JElement_t, JDistance_t>,
223 public virtual JFunction<typename JElement_t::abscissa_type,
224 typename JResultType<typename JElement_t::ordinate_type>::result_type>
237 typedef typename collection_type::iterator
iterator;
265 if (this->size() > 1
u) {
269 if ((p == this->begin() && this->
getDistance(x, (p++)->getX()) > distance_type::precision) ||
270 (p == this->end() && this->
getDistance((--p)->getX(), x) > distance_type::precision)) {
273 return this->getExceptionHandler().action();
277 std::ostringstream os;
279 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range "
280 <<
STREAM(
"?") << x <<
" <> "
281 <<
STREAM(
"?") << this->begin() ->getX() <<
' '
282 <<
STREAM(
"?") << this->rbegin()->getX();
298 }
else if (this->size() == 1
u && this->
getDistance(x, this->begin()->getX()) <= distance_type::precision) {
305 return this->getExceptionHandler().action();
309 std::ostringstream os;
311 os << __FILE__ <<
':' << __LINE__ <<
" not enough data " <<
STREAM(
"?") <<
x;
331 template<
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
337 public JCollection_t<JElement_t, JDistance_t>,
338 public virtual JFunction<typename JElement_t::abscissa_type,
339 typename JResultType<typename JElement_t::ordinate_type>::result_type>
352 typedef typename collection_type::iterator
iterator;
380 if (this->size() > 1
u) {
384 if ((p == this->begin() && this->
getDistance(x, (p++)->getX()) > distance_type::precision) ||
385 (p == this->end() && this->
getDistance((--p)->getX(), x) > distance_type::precision)) {
388 return this->getExceptionHandler().action();
392 std::ostringstream os;
394 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range "
395 <<
STREAM(
"?") << x <<
" <> "
396 <<
STREAM(
"?") << this->begin() ->getX() <<
' '
397 <<
STREAM(
"?") << this->rbegin()->getX();
408 const double dx = this->
getDistance(p->getX(), q->getX());
410 const double b = 1.0 -
a;
422 }
else if (this->size() == 1
u && this->
getDistance(x, this->begin()->getX()) <= distance_type::precision) {
429 return this->getExceptionHandler().action();
433 std::ostringstream os;
435 os << __FILE__ <<
':' << __LINE__ <<
" not enough data " <<
STREAM(
"?") <<
x;
466 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t>
472 public JCollection_t<JElement_t, JDistance_t>,
473 public virtual JFunction<typename JElement_t::abscissa_type,
474 JResultPDF<typename JResultType<typename JElement_t::ordinate_type>::result_type> >
487 typedef typename collection_type::iterator
iterator;
517 if (this->size() <= 1
u) {
520 return this->getExceptionHandler().action();
524 std::ostringstream os;
526 os << __FILE__ <<
':' << __LINE__ <<
" not enough data " <<
STREAM(
"?") <<
x;
534 if (p == this->begin() && this->
getDistance(x, (p++)->getX()) > distance_type::precision) {
538 result = this->getExceptionHandler().action();
543 result.V = this->rbegin()->getIntegral();
547 std::ostringstream os;
549 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range " <<
STREAM(
"?") << x <<
" < " <<
STREAM(
"?") << this->begin() ->getX();
556 }
else if (p == this->end() && this->
getDistance((--p)->getX(), x) > distance_type::precision) {
560 result = this->getExceptionHandler().action();
564 result.v = this->rbegin()->getIntegral();
565 result.V = this->rbegin()->getIntegral();
569 std::ostringstream os;
571 os << __FILE__ <<
':' << __LINE__ <<
" abscissa out of range " <<
STREAM(
"?") << x <<
" > " <<
STREAM(
"?") << this->rbegin()->getX();
582 const int n = std::min((
int) (
N + 1), (
int) this->size());
584 for (
int i = n/2;
i != 0 && p != this->end(); --
i, ++p) {}
585 for (
int i = n ;
i != 0 && p != this->begin(); --
i, --p) {}
590 for (
int i = 0;
i !=
n; ++p, ++
i) {
596 w[
i][2] =
v[
i][2] = p->getIntegral();
598 if (fabs(
u[
i]) < fabs(
u[j])) {
607 result.V = this->rbegin()->getIntegral();
611 for (
int m = 1; m !=
n; ++m) {
613 for (
int i = 0;
i != n-m; ++
i) {
615 const double ho =
u[
i ];
616 const double hp =
u[
i+m];
617 const double dx = ho - hp;
619 for (
int k = 0;
k != 3; ++
k) {
625 v[
i][1] = ho * r[1] - r[0];
626 w[
i][1] = hp * r[1] - r[0];
631 if (2*(j+1) < n - m) {
663 this->begin()->setIntegral(V);
665 for (
iterator j = this->begin(),
i =
j++;
j != this->end(); ++
i, ++
j) {
684 mutable double u[
N+1];
696 template<
unsigned int N,
698 template<
class,
class>
class JCollection_t,
709 template<
unsigned int N,
class JElement_t,
template<
class,
class>
class JCollection_t,
class JDistance_t,
unsigned int M>
715 public JCollection_t<JElement_t, JDistance_t>,
721 typedef JCollection_t<JElement_t, JDistance_t> collection_type;
723 typedef typename collection_type::abscissa_type abscissa_type;
724 typedef typename collection_type::ordinate_type ordinate_type;
725 typedef typename collection_type::value_type value_type;
726 typedef typename collection_type::distance_type distance_type;
728 typedef typename collection_type::const_iterator const_iterator;
729 typedef typename collection_type::const_reverse_iterator const_reverse_iterator;
730 typedef typename collection_type::iterator iterator;
731 typedef typename collection_type::reverse_iterator reverse_iterator;
733 typedef typename JResultType<ordinate_type>::result_type data_type;
734 typedef JFunction<abscissa_type, JResultPolynome<M, data_type> > function_type;
736 typedef typename function_type::argument_type argument_type;
737 typedef typename function_type::result_type result_type;
739 using JFunctional<>::compile;
755 result_type evaluate(const argument_type* pX) const
757 const argument_type x = *pX;
759 if (this->size() <= N) {
761 std::ostringstream os;
763 os << __FILE__ << ':' << __LINE__ << " not enough data " << STREAM("?") << x;
765 throw JFunctionalException(os.str());
768 const_iterator p = this->lower_bound(x);
770 if ((p == this->begin() && this->getDistance(x, (p++)->getX()) > distance_type::precision) ||
771 (p == this->end() && this->getDistance((--p)->getX(), x) > distance_type::precision)) {
773 std::ostringstream os;
775 os << __FILE__ << ':' << __LINE__ << " abscissa out of range "
776 << STREAM("?") << x << " <> "
777 << STREAM("?") << this->begin() ->getX() << ' '
778 << STREAM("?") << this->rbegin()->getX();
780 throw JValueOutOfRange(os.str());
786 const int n = std::min((int) (N + 1), (int) this->size());
788 for (int i = n/2; i != 0 && p != this->end(); --i, ++p) {}
789 for (int i = n ; i != 0 && p != this->begin(); --i, --p) {}
794 for (int i = 0; i != n; ++p, ++i) {
796 u[i] = this->getDistance(x, p->getX());
798 w[i][0] = v[i][0] = JFunctional<argument_type, data_type>::getValue(p->getY(), pX);
800 for (unsigned int k = 1; k != M+1; ++k) {
801 w[i][k] = v[i][k] = JMATH::zero;
804 if (fabs(u[i]) < fabs(u[j])) {
810 for (unsigned int k = 0; k != M+1; ++k) {
811 result.y[k] = v[j][k];
816 for (int m = 1; m != n; ++m) {
818 for (int i = 0; i != n-m; ++i) {
820 const double ho = u[ i ];
821 const double hp = u[i+m];
822 const double dx = ho - hp;
824 for (int k = 0; k != M+1; ++k) {
825 r[k] = (v[i+1][k] - w[i][k]) / dx;
831 for (int k = 1; k != M+1; ++k) {
832 v[i][k] = ho * r[k] - k * r[k-1];
833 w[i][k] = hp * r[k] - k * r[k-1];
837 if (2*(j+1) < n - m) {
839 for (int k = 0; k != M+1; ++k) {
840 result.y[k] += v[j+1][k];
845 for (int k = 0; k != M+1; ++k) {
846 result.y[k] += w[j][k];
860 virtual void do_compile() override
865 mutable double u[N+1];
866 mutable data_type v[N+1][M+1];
867 mutable data_type w[N+1][M+1];
868 mutable data_type r[M+1];
870 mutable result_type result;
877 template<unsigned int N, class JElement_t, template<class, class> class JCollection_t, class JDistance_t, unsigned int M>
878 class JPolintFunction<N,
881 JResultPolynome<M, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
883 public JPolintCollection<N,
886 JResultPolynome<M, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
888 public virtual JFunction<typename JElement_t::abscissa_type,
889 JResultPolynome<N, typename JResultType<typename JElement_t::ordinate_type>::result_type> >
893 typedef JPolintCollection<N,
896 JResultPolynome<M, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
897 JDistance_t> collection_type;
899 typedef typename collection_type::abscissa_type abscissa_type;
900 typedef typename collection_type::ordinate_type ordinate_type;
901 typedef typename collection_type::value_type value_type;
902 typedef typename collection_type::distance_type distance_type;
904 typedef typename collection_type::const_iterator const_iterator;
905 typedef typename collection_type::const_reverse_iterator const_reverse_iterator;
906 typedef typename collection_type::iterator iterator;
907 typedef typename collection_type::reverse_iterator reverse_iterator;
909 typedef typename JResultType<ordinate_type>::result_type data_type;
910 typedef JFunction<abscissa_type, JResultPolynome<M, data_type> > function_type;
912 typedef typename function_type::argument_type argument_type;
913 typedef typename function_type::result_type result_type;
914 typedef typename function_type::JExceptionHandler exceptionhandler_type;
930 virtual result_type evaluate(const argument_type* pX) const override
933 return collection_type::evaluate(pX);
935 catch(const JException& error) {
936 return this->getExceptionHandler().action(error);
945 template<unsigned int N, class JElement_t, template<class, class> class JCollection_t, class JDistance_t>
946 class JPolintFunction<N,
949 JResultDerivative<typename JResultType<typename JElement_t::ordinate_type>::result_type>,
951 public JPolintCollection<N,
954 JResultPolynome<1, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
956 public virtual JFunction<typename JElement_t::abscissa_type,
957 JResultDerivative<typename JResultType<typename JElement_t::ordinate_type>::result_type> >
961 typedef JPolintCollection<N,
964 JResultPolynome<1, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
965 JDistance_t> collection_type;
967 typedef typename collection_type::abscissa_type abscissa_type;
968 typedef typename collection_type::ordinate_type ordinate_type;
969 typedef typename collection_type::value_type value_type;
970 typedef typename collection_type::distance_type distance_type;
972 typedef typename collection_type::const_iterator const_iterator;
973 typedef typename collection_type::const_reverse_iterator const_reverse_iterator;
974 typedef typename collection_type::iterator iterator;
975 typedef typename collection_type::reverse_iterator reverse_iterator;
977 typedef typename JResultType<ordinate_type>::result_type data_type;
978 typedef JFunction<abscissa_type, JResultDerivative<data_type> > function_type;
980 typedef typename function_type::argument_type argument_type;
981 typedef typename function_type::result_type result_type;
982 typedef typename function_type::JExceptionHandler exceptionhandler_type;
984 using JFunctional<>::compile;
1000 virtual result_type evaluate(const argument_type* pX) const override
1003 return collection_type::evaluate(pX);
1005 catch(const JException& error) {
1006 return this->getExceptionHandler().action(error);
1015 template<unsigned int N, class JElement_t, template<class, class> class JCollection_t, class JDistance_t>
1016 class JPolintFunction<N,
1019 JResultHesse<typename JResultType<typename JElement_t::ordinate_type>::result_type>,
1021 public JPolintCollection<N,
1024 JResultPolynome<2, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
1026 public virtual JFunction<typename JElement_t::abscissa_type,
1027 JResultHesse<typename JResultType<typename JElement_t::ordinate_type>::result_type> >
1031 typedef JPolintCollection<N,
1034 JResultPolynome<2, typename JResultType<typename JElement_t::ordinate_type>::result_type>,
1035 JDistance_t> collection_type;
1037 typedef typename collection_type::abscissa_type abscissa_type;
1038 typedef typename collection_type::ordinate_type ordinate_type;
1039 typedef typename collection_type::value_type value_type;
1040 typedef typename collection_type::distance_type distance_type;
1042 typedef typename collection_type::const_iterator const_iterator;
1043 typedef typename collection_type::const_reverse_iterator const_reverse_iterator;
1044 typedef typename collection_type::iterator iterator;
1045 typedef typename collection_type::reverse_iterator reverse_iterator;
1047 typedef typename JResultType<ordinate_type>::result_type data_type;
1048 typedef JFunction<abscissa_type, JResultHesse<data_type> > function_type;
1050 typedef typename function_type::argument_type argument_type;
1051 typedef typename function_type::result_type result_type;
1052 typedef typename function_type::JExceptionHandler exceptionhandler_type;
1054 using JFunctional<>::compile;
1070 virtual result_type evaluate(const argument_type* pX) const override
1073 return collection_type::evaluate(pX);
1075 catch(const JException& error) {
1076 return this->getExceptionHandler().action(error);
1087 template<unsigned int N,
1089 template<class, class> class JCollection_t,
1090 class JResult_t = typename JElement_t::ordinate_type,
1091 class JDistance_t = JDistance<typename JElement_t::abscissa_type> >
1092 class JPolintFunction1D :
1093 public JPolintFunction<N, JElement_t, JCollection_t, JResult_t, JDistance_t>,
1094 public virtual JFunction1D<typename JElement_t::abscissa_type, JResult_t>
1098 typedef JCollection_t<JElement_t, JDistance_t> collection_type;
1100 typedef typename collection_type::abscissa_type abscissa_type;
1101 typedef typename collection_type::ordinate_type ordinate_type;
1102 typedef typename collection_type::value_type value_type;
1103 typedef typename collection_type::distance_type distance_type;
1105 typedef typename collection_type::const_iterator const_iterator;
1106 typedef typename collection_type::const_reverse_iterator const_reverse_iterator;
1107 typedef typename collection_type::iterator iterator;
1108 typedef typename collection_type::reverse_iterator reverse_iterator;
1110 typedef JFunction1D<abscissa_type, JResult_t> function_type;
1112 typedef typename function_type::argument_type argument_type;
1113 typedef typename function_type::result_type result_type;
1114 typedef typename function_type::JExceptionHandler exceptionhandler_type;
1130 template<class JAbscissa_t, class JOrdinate_t>
1133 template<template<class, class, class> class JMap_t>
1134 struct JMapCollection;
1140 template<unsigned int N,
1143 template<class, class, class> class JMap_t,
1145 class JDistance_t = JDistance<JKey_t> >
1147 public JPolintFunction<N,
1148 JElement2D<JKey_t, JValue_t>,
1149 JMapCollection<JMap_t>::template collection_type,
1155 typedef JElement2D<JKey_t, JValue_t> element_type;
1156 typedef JPolintFunction<N,
1158 JMapCollection<JMap_t>::template collection_type,
1160 JDistance_t> JPolintFunction_t;
1162 typedef typename JPolintFunction_t::abscissa_type abscissa_type;
1163 typedef typename JPolintFunction_t::ordinate_type ordinate_type;
1164 typedef typename JPolintFunction_t::value_type value_type;
1165 typedef typename JPolintFunction_t::distance_type distance_type;
1167 typedef typename JPolintFunction_t::const_iterator const_iterator;
1168 typedef typename JPolintFunction_t::const_reverse_iterator const_reverse_iterator;
1169 typedef typename JPolintFunction_t::iterator iterator;
1170 typedef typename JPolintFunction_t::reverse_iterator reverse_iterator;
1172 typedef typename JPolintFunction_t::argument_type argument_type;
1173 typedef typename JPolintFunction_t::result_type result_type;
1174 typedef typename JPolintFunction_t::JExceptionHandler exceptionhandler_type;
1194 template<unsigned int N,
1196 template<class, class> class JCollection_t,
1199 inline typename JElement_t::ordinate_type
1200 integrate(const JPolintFunction1D<N, JElement_t, JCollection_t, JResult_t, JDistance_t>& input,
1201 typename JMappable<JElement_t>::map_type& output)
1203 return integrate(input, output, JLANG::JBool<N == 0 || N == 1>());
1218 template<unsigned int N,
1220 template<class, class> class JCollection_t,
1223 inline typename JElement_t::ordinate_type
1224 integrate(const JPolintFunction1D<N, JElement_t, JCollection_t, JResult_t, JDistance_t>& input,
1225 typename JMappable<JElement_t>::map_type& output,
1226 const JLANG::JBool<false>& option)
1228 typedef typename JElement_t::abscissa_type abscissa_type;
1229 typedef typename JElement_t::ordinate_type ordinate_type;
1230 typedef typename JPolintFunction1D<N, JElement_t, JCollection_t, JResult_t, JDistance_t>::const_iterator const_iterator;
1232 ordinate_type V(JMATH::zero);
1234 if (input.getSize() > 1) {
1236 output.put(input.begin()->getX(), V);
1238 const JGaussLegendre engine(N);
1240 for (const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
1242 const abscissa_type xmin = i->getX();
1243 const abscissa_type xmax = j->getX();
1245 for (JGaussLegendre::const_iterator m = engine.begin(); m != engine.end(); ++m) {
1247 const abscissa_type x = 0.5 * (xmax + xmin + m->getX() * (xmax - xmin));
1248 const ordinate_type v = 0.5 * (xmax - xmin) * m->getY() * get_value(input(x));
1253 output.put(xmax, V);
1271 template<class JElement_t,
1272 template<class, class> class JCollection_t,
1275 inline typename JElement_t::ordinate_type
1276 integrate(const JPolintFunction1D<0, JElement_t, JCollection_t, JResult_t, JDistance_t>& input,
1277 typename JMappable<JElement_t>::map_type& output,
1278 const JLANG::JBool<true>& option)
1280 typedef typename JElement_t::ordinate_type ordinate_type;
1281 typedef typename JPolintFunction1D<0, JElement_t, JCollection_t, JResult_t, JDistance_t>::const_iterator const_iterator;
1283 ordinate_type V(JMATH::zero);
1285 if (input.getSize() > 1) {
1287 output.put(input.begin()->getX(), V);
1289 for (const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
1291 V += input.getDistance(i->getX(), j->getX()) * j->getY();
1293 output.put(j->getX(), V);
1311 template<class JElement_t,
1312 template<class, class> class JCollection_t,
1315 inline typename JElement_t::ordinate_type
1316 integrate(const JPolintFunction1D<1, JElement_t, JCollection_t, JResult_t, JDistance_t>& input,
1317 typename JMappable<JElement_t>::map_type& output,
1318 const JLANG::JBool<true>& option)
1320 typedef typename JElement_t::ordinate_type ordinate_type;
1321 typedef typename JPolintFunction1D<1, JElement_t, JCollection_t, JResult_t, JDistance_t>::const_iterator const_iterator;
1323 ordinate_type V(JMATH::zero);
1325 if (input.getSize() > 1) {
1327 output.put(input.begin()->getX(), V);
1329 for (const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
1331 V += 0.5 * input.getDistance(i->getX(), j->getX()) * (i->getY() + j->getY());
1333 output.put(j->getX(),
V);
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Exception for a functional operation.
The elements in a collection are sorted according to their abscissa values and a given distance opera...
This include file containes various data structures that can be used as specific return types for the...
static const JZero zero
Function object to assign zero value.
Generation of compiler error.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
size_t getSize(T(&array)[N])
Get size of c-array.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Auxiliary data structure for handling std::ostream.
Exception for numerical precision error.
Auxiliary classes for numerical integration.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Exception for accessing a value in a collection that is outside of its range.