Speed up clr

ivmodels/tests/conditional_likelihood_ratio.py

@@ -1,3 +1,6 @@
+import ctypes
+from pathlib import Path
+
 import numpy as np
 import scipy
 
@@ -82,21 +85,35 @@ def conditional_likelihood_ratio_critical_value_function(
         return 1 - scipy.stats.chi2(q).cdf(z)
 
     if method in ["numerical_integration"]:
-        beta = scipy.stats.beta((q - p) / 2, p / 2)
-        chi2 = scipy.stats.chi2(q)
+        # Load the shared library
+        lib = ctypes.CDLL(Path(__file__).parent / "integrand.dylib")
+        integrand = lib.integrand
+        integrand.restype = ctypes.c_double
+        integrand.argtypes = (ctypes.c_double, ctypes.c_void_p)
+
         a = s_min / (z + s_min)
+        alpha = (q - p) / 2
+        beta = p / 2
+        # Create a NumPy array with the parameter values
+        params = np.array([a, z, alpha, beta, q], dtype=np.double)
+
+        # Convert the parameter array to a ctypes pointer
+        params_ctypes = params.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
+
+        # Cast the ctypes pointer to a void pointer
+        user_data = ctypes.cast(params_ctypes, ctypes.c_void_p)
 
-        def integrand(b):
-            return beta.pdf(b) * chi2.cdf(z / (1 - a * b))
+        # Create the LowLevelCallable with the user_data
+        integrand_callable = scipy.LowLevelCallable(integrand, user_data)
 
         res = scipy.integrate.quad(
-            integrand,
+            integrand_callable,
             0,
             1,
             epsabs=tol,
         )
 
-        return 1 - res[0]
+        return 1 - res[0] / scipy.special.beta(alpha, beta)
 
     elif method == "power_series":
         a = s_min / (z + s_min)

ivmodels/tests/integrand.c

@@ -0,0 +1,28 @@
+#include <math.h>
+#include <gsl/gsl_cdf.h>
+#include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_randist.h>
+
+// Parameters struct to hold the values of a and z
+typedef struct {
+    double a;
+    double z;
+    double alpha;
+    double beta;
+} Params;
+
+// C function to compute the integrand
+double integrand(double x, void *user_data) {
+    double *p = (double *)user_data;
+
+    double a = p[0];
+    double z = p[1];
+    double alpha = p[2];
+    double beta = p[3];
+    int k = (int) p[4];
+
+    double beta_pdf = exp(log(x) * (alpha - 1) + log(1 - x) * (beta - 1));
+    double chi2_cdf = gsl_cdf_chisq_P(z / (1.0 - a * x), k);
+
+    return beta_pdf * chi2_cdf;
+}

setup.py

@@ -0,0 +1,19 @@
+from setuptools import setup, Extension
+from Cython.Build import cythonize
+import numpy as np
+
+ext_modules = [
+    Extension(
+        "integrand",
+        sources=["ivmodels/tests/integrand.pyx"],
+        libraries=["gsl", "gslcblas"],
+        library_dirs=["/opt/homebrew/lib"],
+        include_dirs=["/opt/homebrew/include", np.get_include()],
+        extra_compile_args=["-fPIC"]
+    )
+]
+
+setup(
+    name="integrand",
+    ext_modules=cythonize(ext_modules),
+)