Integrate with hints

include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp

@@ -76,7 +76,9 @@ class AdaptiveQuadratureBase
 	 * @return Approximation of the integral.
 	 */
 	template<class Function> LongScalar integrate(const Function& f, const Scalar& xmin, const Scalar& xmax);
-
+
+	template<class Function> LongScalar integrateWithHints(const Function& f, const std::span<const Scalar> mu, const Scalar& sigma);
+
 	/**
 	 * @brief Perform adaptive quadrature on (-inf, xmax].
 	 * 
@@ -168,6 +170,8 @@ class AdaptiveQuadratureBase
 
 	constexpr std::span<const Interval> getSubIntervals() const { return m_intervals; }
 private:
+	template<class Function> LongScalar adaptQuadrature(const Function& func);
+
 	std::vector<Interval>   m_intervals;
 	std::vector<LongScalar> m_subIntergrals;
 	std::vector<LongScalar> m_subIntergralsErr;

include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp

@@ -28,31 +28,64 @@ AdaptiveQuadratureBase<Derived>::AdaptiveQuadratureBase(const Size& maxIt, const
 }
 
 template<class Derived> template<class Function> 
-auto AdaptiveQuadratureBase<Derived>::integrate(const Function& f, const Scalar& xmin, const Scalar& xmax) -> LongScalar
-{	
-	using std::ceil;
+auto AdaptiveQuadratureBase<Derived>::adaptQuadrature(const Function& f) -> LongScalar
+{
 	using std::abs;
 	using std::isfinite;
 
 	using const_Iterator = typename std::vector<LongScalar>::const_iterator;
-	
+
 	m_hasConverged = false;
+
+	if (m_out) { fmt::print(m_out, "#Iteration integral estimated_error relative_tol absolute_tol\n"); }
+
+	for (m_it=0; m_it!=m_maxIt; ++m_it)
+	{
+		const LongScalar I   = getEstimatedIntegral();
+		const LongScalar err = getEstimatedError();
+
+		if (m_out) { fmt::print(m_out, "{} {:10.4e} {:10.4e} {:10.4e} {:10.4e}\n", m_it, Scalar(I), Scalar(err), Scalar(abs(I))*m_relativeTol, Scalar(m_absoluteTol)); }
+
+		if (not isfinite(I)) { return I; }
+		if (err < abs(I)*LongScalar(m_relativeTol) or err < LongScalar(m_absoluteTol)) { m_hasConverged = true; return I; }
+
+		// we find the interval over which the integral is the least accurate
+		const const_Iterator maxErrIt = std::ranges::max_element(m_subIntergralsErr);
+		const Size maxErrIdx = Size(std::ranges::distance(m_subIntergralsErr.begin(), maxErrIt));
+		// we split it in two
+		const auto& [a, b] = m_intervals[maxErrIdx];
+
+		Scalar midPoint = std::midpoint(a, b); // non-const because I want to move it when I do not need it.
+		// first interval
+		m_intervals[maxErrIdx] = Interval(a, midPoint);
+		std::tie(m_subIntergrals[maxErrIdx], m_subIntergralsErr[maxErrIdx]) = estimateIntegral(f, a, midPoint);
+		// second interval
+		const auto [res, estimatedErr] = estimateIntegral(f, midPoint, b);
+		m_intervals.emplace_back(std::move(midPoint), b);
+		m_subIntergrals.push_back(res);
+		m_subIntergralsErr.push_back(estimatedErr);
+	}
+	return getEstimatedIntegral();
+}
+
+template<class Derived> template<class Function> 
+auto AdaptiveQuadratureBase<Derived>::integrate(const Function& f, const Scalar& xmin, const Scalar& xmax) -> LongScalar
+{	
+	using std::ceil;
+
 	m_intervals.clear();
 	m_subIntergrals.clear();
 	m_subIntergralsErr.clear();
 
 	LongScalar res;
 	LongScalar estimatedErr;
 
-	if (m_out) { fmt::print(m_out, "#NumericalIntegrator addapting quadrature over [{}, {}]\n", xmin, xmax); }
-	if (m_out) { fmt::print(m_out, "#Iteration integral estimated_error relative_tol absolute_tol\n"); }
-
 	const Size N = Size(ceil(getMaxDeltaX(xmin, xmax)));
 
 	m_intervals.reserve(N);
 	m_subIntergrals.reserve(N);
 	m_subIntergralsErr.reserve(N);
-	
+
 	for (Size i=0; i!=N; ++i)
 	{
 		const Scalar x_i   = xmin + Scalar(i)*(xmax - xmin) / Scalar(N);
@@ -64,34 +97,89 @@ auto AdaptiveQuadratureBase<Derived>::integrate(const Function& f, const Scalar&
 		m_subIntergrals.push_back(res);
 		m_subIntergralsErr.push_back(estimatedErr);
 	}
+
+	if (m_out) { fmt::print(m_out, "#NumericalIntegrator adapting quadrature over [{}, {}]\n", xmin, xmax); }
+	return adaptQuadrature(f);
+}
+
+template<class Derived> template<class Function> 
+auto AdaptiveQuadratureBase<Derived>::integrateWithHints(const Function& f, const std::span<const Scalar> mu, const Scalar& sigma) -> LongScalar
+{
+	using std::abs;
+	using std::isfinite;
+	using std::swap;
+
+	using Iterator = typename std::vector<Interval>::iterator;
 
-	for (m_it=0; m_it!=m_maxIt; ++m_it)
+	constexpr Size scal = 3;
+
+	assert(sigma > Scalar{});
+
+	GaussLaguerreQuadrature<Scalar,LongScalar> gLaguerreQuad;
+
+	m_intervals.clear();
+	m_subIntergrals.clear();
+	m_subIntergralsErr.clear();
+
+	m_intervals.reserve(2*mu.size());
+
+	// First pass, we compute intervals near the peaks on which the function is not null.
+	for (const double& mu_i : mu)
 	{
-		const LongScalar I   = getEstimatedIntegral();
-		const LongScalar err = getEstimatedError();
+		Interval curr(mu_i - scal*sigma, mu_i + scal*sigma);
+
+		if (curr.first > curr.second) { swap(curr.first, curr.second); }
 
-		if (m_out) { fmt::print(m_out, "{} {:10.4e} {:10.4e} {:10.4e} {:10.4e}\n", m_it, Scalar(I), Scalar(err), Scalar(abs(I))*m_relativeTol, Scalar(m_absoluteTol)); }
+		LongScalar leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, curr.first);
+		while (isfinite(leftIntegral) and abs(leftIntegral) >= NumTraits<LongScalar>::epsilon)
+		{
+			curr.first -= scal*sigma;
+			leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, curr.first);
+			assert(curr.first <= curr.second);
+		}
 
-		if (not isfinite(I)) { return I; }
-		if (err < abs(I)*LongScalar(m_relativeTol) or err < LongScalar(m_absoluteTol)) { m_hasConverged = true; return I; }
+		if (not isfinite(leftIntegral)) { return NumTraits<LongScalar>::NaN; }
 
-		// we find the interval over which the integral is the least accurate
-		const const_Iterator maxErrIt = std::ranges::max_element(m_subIntergralsErr);
-		const Size maxErrIdx = Size(std::ranges::distance(m_subIntergralsErr.begin(), maxErrIt));
-		// we split it in two
-		const auto& [a, b] = m_intervals[maxErrIdx];
+		LongScalar rightIntegral = gLaguerreQuad.integrateRightInfinite(f, curr.second);
+		while (isfinite(rightIntegral) and abs(rightIntegral) >= NumTraits<LongScalar>::epsilon)
+		{
+			curr.second += scal*sigma;
+			rightIntegral = gLaguerreQuad.integrateRightInfinite(f, curr.second);
+			assert(curr.first <= curr.second);
+		}
+
+		if (not isfinite(rightIntegral)) { return NumTraits<LongScalar>::NaN; }
+
+		Iterator firstInterval = m_intervals.begin();
+		while (firstInterval != m_intervals.end() and firstInterval->second < curr.first) { ++firstInterval; }
+
+		Iterator boundInterval = firstInterval;
+		for ( ;boundInterval != m_intervals.end() and boundInterval->first <= curr.second; ++boundInterval)
+		{
+			curr.first  = std::min(curr.first,  boundInterval->first);
+			curr.second = std::max(curr.second, boundInterval->second);
+		}
+
+		const Iterator it = m_intervals.erase(firstInterval, boundInterval);
+		m_intervals.insert(it, curr);
+	}
+
+	LongScalar res;
+	LongScalar estimatedErr;	
+
+	m_subIntergrals.reserve(m_intervals.size());
+	m_subIntergralsErr.reserve(m_intervals.size());
+
+	for (const auto& [xmin, xmax] : m_intervals)
+	{
+		std::tie(res, estimatedErr) = estimateIntegral(f, xmin, xmax);
 
-		Scalar midPoint = std::midpoint(a, b); // non-const because I want to move it when I do not need it.
-		// first interval
-		m_intervals[maxErrIdx] = Interval(a, midPoint);
-		std::tie(m_subIntergrals[maxErrIdx], m_subIntergralsErr[maxErrIdx]) = estimateIntegral(f, a, midPoint);
-		// second interval
-		std::tie(res, estimatedErr) = estimateIntegral(f, midPoint, b);
-		m_intervals.emplace_back(std::move(midPoint), b);
 		m_subIntergrals.push_back(res);
 		m_subIntergralsErr.push_back(estimatedErr);
 	}
-	return getEstimatedIntegral();
+
+	if (m_out) { fmt::print(m_out, "#NumericalIntegrator adapting quadrature over {}\n", m_intervals); }
+	return adaptQuadrature(f); 
 }
 
 template<class Derived> template<class Function>
@@ -100,7 +188,6 @@ auto AdaptiveQuadratureBase<Derived>::integrateLeftInfinite(const Function& f, c
 	using std::isfinite;
 	using std::abs;
 
-
 	GaussLaguerreQuadrature<Scalar,LongScalar> gLaguerreQuad;
 
 	Scalar xmin = -1;