I've wrapped the 'NumericalIntegration' C++ library in Haskell. Here is the latest version of the package (the version on Hackage is older).
Here is the main part of the C++ code :
class Integrand {
private:
std::function<double(double)> f;
public:
Integrand(std::function<double(double)>& f_) : f(f_) {}
double operator()(const double x) const { return f(x); }
};
double integration(double f(double),
double lower,
double upper,
double relError,
int subdiv,
double* errorEstimate,
int* errorCode) {
// Define the integrand.
std::function<double(double)> f_ = [&](double x) { return f(x); };
Integrand integrand(f_);
// Define the integrator.
Eigen::Integrator<double> integrator(subdiv);
// Define a quadrature rule.
Eigen::Integrator<double>::QuadratureRule rule =
Eigen::Integrator<double>::GaussKronrod201;
// Define the desired absolute error.
double absError = 0.0;
// Integrate.
double result = integrator.quadratureAdaptive(integrand, lower, upper,
absError, relError, rule);
*errorEstimate = integrator.estimatedError();
*errorCode = integrator.errorCode();
return result;
}
And here is the main part of the Haskell code:
foreign import ccall safe "wrapper" funPtr
:: (Double -> Double) -> IO (FunPtr (Double -> Double))
foreign import ccall safe "integration" c_integration
:: FunPtr (Double -> Double) -> Double -> Double -> Double -> Int
-> Ptr Double -> Ptr Int -> IO Double
-- | Numerical integration.
integration :: (Double -> Double) -- ^ integrand
-> Double -- ^ lower bound
-> Double -- ^ upper bound
-> Double -- ^ desired relative error
-> Int -- ^ number of subdivisions
-> IO IntegralResult -- ^ value, error estimate, error code
integration f lower upper relError subdiv = do
errorEstimatePtr <- mallocBytes (sizeOf (0 :: Double))
errorCodePtr <- mallocBytes (sizeOf (0 :: Int))
fPtr <- funPtr f
result <-
c_integration fPtr lower upper relError subdiv errorEstimatePtr errorCodePtr
errorEstimate <- peek errorEstimatePtr
errorCode <- peek errorCodePtr
let out = IntegralResult {_value = result, _error = errorEstimate, _code = errorCode}
free errorEstimatePtr
free errorCodePtr
freeHaskellFunPtr fPtr
return out
This works but there's a problem regarding the error code of the integration. When the integration is ok, the error code should be 0. Sometimes it is 0, as expected. But sometimes it is a huge integer number, nonsensical, though the integration is fine.
Would you have an idea about this issue? Why this strange error code? Is there something bad in my code? I'm not fluent in C++ (nor in Haskell). But apart this strange error code, the library seems to work very well.
On the x86_64
platform, a Haskell Int
is 64 bits while a C int
is 32 bits. (A C long
is 64 bits.) In your code, you are picking up garbage in the upper bytes and getting an absurdly large 64-bit integer whose lowest 32-bits are zero, no doubt.
Anyway, there's a module Foreign.C
that contains newtypes CInt
, CDouble
, etc. that are intended to match the corresponding C types on the target platform, and I think it's considered good practice to always use those:
import Foreign
import Foreign.C
foreign import ccall safe "wrapper" funPtr
:: (CDouble -> CDouble) -> IO (FunPtr (CDouble -> CDouble))
foreign import ccall safe "integration" c_integration
:: FunPtr (CDouble -> CDouble) -> CDouble -> CDouble -> CDouble -> CInt
-> Ptr CDouble -> Ptr CInt -> IO CDouble
Of course, since these are newtypes, there is the bother of wrapping and unwrapping, though usually fromIntegral
takes care of that automatically, while also converting across bit sizes:
errorCode <- fromIntegral <$> peek errorCodePtr
But, as a less portable alternative, you could stick with Double
and Int32
in your foreign import
declarations, provided you are only targeting platforms where C int
s are 32 bits.
Also note that, if you get the types right, then malloc
is a type-safe alternative to mallocBytes
:
errorCodePtr <- malloc
The type of the Ptr
determines the correct number of bytes to allocate here.