Skip to content

Commit

Permalink
Remove boost accumulator dependence. Revert back to simple summation.
Browse files Browse the repository at this point in the history
  • Loading branch information
yjian012 committed Jan 8, 2024
1 parent 377f578 commit dacef8b
Showing 1 changed file with 24 additions and 29 deletions.
53 changes: 24 additions & 29 deletions ql/math/integrals/discreteintegrals.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,10 +18,6 @@
*/

#include <ql/math/integrals/discreteintegrals.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/sum.hpp>

using namespace boost::accumulators;

namespace QuantLib {

Expand All @@ -31,13 +27,13 @@ namespace QuantLib {
const Size n = f.size();
QL_REQUIRE(n == x.size(), "inconsistent size");

accumulator_set<Real, features<tag::sum> > acc;
Real sum = 0.0;

for (Size i=0; i < n-1; ++i) {
acc((x[i+1]-x[i])*(f[i]+f[i+1]));
sum += (x[i+1]-x[i])*(f[i]+f[i+1]);
}

return 0.5*sum(acc);
return 0.5*sum;
}

Real DiscreteSimpsonIntegral::operator()(
Expand All @@ -46,7 +42,7 @@ namespace QuantLib {
const Size n = f.size();
QL_REQUIRE(n == x.size(), "inconsistent size");

accumulator_set<Real, features<tag::sum> > acc;
Real sum = 0.0;

for (Size j=0; j < n-2; j+=2) {
const Real dxj = x[j+1]-x[j];
Expand All @@ -58,61 +54,60 @@ namespace QuantLib {
const Real beta = dd*dd;
const Real gamma = dxj*(2*dxjp1-dxj);

acc(k*(alpha*f[j]+beta*f[j+1]+gamma*f[j+2]));
sum += k*(alpha*f[j]+beta*f[j+1]+gamma*f[j+2]);
}
if ((n & 1) == 0U) {
acc(0.5*(x[n-1]-x[n-2])*(f[n-1]+f[n-2]));
sum += 0.5*(x[n-1]-x[n-2])*(f[n-1]+f[n-2]);
}

return sum(acc);
return sum;
}

Real DiscreteTrapezoidIntegrator::integrate(
const ext::function<Real (Real)>& f, Real a, Real b) const {
const Size n=maxEvaluations()-1;
const Real d=(b-a)/n;

accumulator_set<Real, features<tag::sum> > acc;
Real sum = f(a)/2;

acc(f(a)/2);
for(Size i=0;i<n-1;++i) {
a+=d;
acc(f(a));
a += d;
sum += f(a);
}
acc(f(b)/2);
sum += f(b)/2;

increaseNumberOfEvaluations(maxEvaluations());

return d*sum(acc);
return d * sum;
}

Real DiscreteSimpsonIntegrator::integrate(
const ext::function<Real (Real)>& f, Real a, Real b) const {
const Size n=maxEvaluations()-1;
const Real d=(b-a)/n, d2=d*2;

Real sum=0.0, x=a+d;
Real sum = 0.0, x = a + d;
for(Size i=1;i<n;i+=2) {//to time 4
sum+=f(x);
x+=d2;
sum += f(x);
x += d2;
}
sum*=2;
sum *= 2;

x=a+d2;
x = a+d2;
for(Size i=2;i<n-1;i+=2) {//to time 2
sum+=f(x);
x+=d2;
sum += f(x);
x += d2;
}
sum*=2;
sum *= 2;

sum+=f(a);
sum += f(a);
if(n&1)
sum+=1.5*f(b)+2.5*f(b-d);
sum += 1.5*f(b)+2.5*f(b-d);
else
sum+=f(b);
sum += f(b);

increaseNumberOfEvaluations(maxEvaluations());

return d/3*sum;
return d/3 * sum;
}
}

0 comments on commit dacef8b

Please sign in to comment.