Merge branch 'remove-gsl'

.github/workflows/build-and-test.yaml

@@ -20,7 +20,6 @@ jobs:
       - name: Install dependencies
         run: |
           sudo apt-get update
-          sudo apt-get install -y libgsl-dev libgslcblas0
           pip install --upgrade pip setuptools wheel auditwheel patchelf
       - name: Build the wheel
         run: python setup.py bdist_wheel
@@ -46,7 +45,6 @@ jobs:
       - name: Install dependencies
         shell: powershell
         run: |
-          vcpkg install gsl:x64-windows-static
           pip install --upgrade pip setuptools wheel
       - name: Build the wheel
         run: python setup.py bdist_wheel

Makefile

@@ -2,7 +2,7 @@
 .PHONY: clean bdist
 
 traceon/backend/traceon_backend.so: traceon/backend/traceon-backend.c
-	clang -O3 -march=native -shared -fPIC -ffast-math -Wno-extern-initializer ./traceon/backend/traceon-backend.c -o traceon/backend/traceon_backend.so -lgsl -lm 
+	clang -O3 -march=native -shared -fPIC -ffast-math -Wno-extern-initializer ./traceon/backend/traceon-backend.c -o traceon/backend/traceon_backend.so -lm 
 
 clean:
 	rm -rf ./build ./dist

setup.py

@@ -3,22 +3,20 @@
 from setuptools import setup, Extension
 import platform
 
-VCPKG_DIR = 'C:\\vcpkg'
-
 if platform.system() == 'Linux':
     os.environ['CC'] = 'clang' # Clang is known to produce faster binaries.
-    compiler_kwargs = dict(extra_compile_args=['-O3', '-mavx', '-ffast-math', '-DNDEBUG'], extra_link_args=['-lm', '-l:libgsl.a', '-l:libgslcblas.a'])
+    compiler_kwargs = dict(extra_compile_args=['-O3', '-mavx', '-ffast-math', '-DNDEBUG'], extra_link_args=['-lm'])
     extra_objects = []
     include_dirs = []
 elif platform.system() == 'Darwin':
     os.environ['CC'] = 'clang' # Clang is known to produce faster binaries.
-    compiler_kwargs = dict(extra_compile_args=['-O3', '-mavx', '-ffast-math', '-DNDEBUG'], extra_link_args=['-lm', '-L/opt/homebrew/lib', '-lgsl', '-lgslcblas'])
+    compiler_kwargs = dict(extra_compile_args=['-O3', '-mavx', '-ffast-math', '-DNDEBUG'], extra_link_args=['-lm', '-L/opt/homebrew/lib'])
     extra_objects = []
     include_dirs = ['/opt/homebrew/include/']
 elif platform.system() == 'Windows':
     compiler_kwargs = dict(extra_compile_args=['/fp:fast', '/Ox', '/Ob3', '/Oi', '/GL', '/arch:AVX', '-I .\\traceon\\backend\\'])
-    extra_objects = [VCPKG_DIR + '\\packages\\gsl_x64-windows-static\\lib\\gsl.lib']
-    include_dirs = VCPKG_DIR + '\\packages\\gsl_x64-windows-static\\include'
+    extra_objects = []
+    include_dirs = []
 
 
 backend_extension = Extension(

traceon/backend/__init__.py

@@ -8,6 +8,8 @@
 
 from numpy.ctypeslib import ndpointer
 import numpy as np
+from scipy import LowLevelCallable
+from scipy.integrate import quad
 
 from .. import logging
 
@@ -182,6 +184,8 @@ def __init__(self, eff, *args, **kwargs):
     'current_field': (None, v3, v3, currents_2d, jac_buffer_3d, pos_buffer_3d, sz),
     'current_axial_derivatives_radial_ring': (None, arr(ndim=2), currents_2d, jac_buffer_3d, pos_buffer_3d, sz, z_values, sz),
     'fill_jacobian_buffer_radial': (None, jac_buffer_2d, pos_buffer_2d, vertices, sz),
+    'self_potential_radial': (dbl, dbl, vp),
+    'self_field_dot_normal_radial': (dbl, dbl, vp),
     'fill_matrix_radial': (None, arr(ndim=2), lines, arr(dtype=C.c_uint8, ndim=1), arr(ndim=1), jac_buffer_2d, pos_buffer_2d, sz, sz, C.c_int, C.c_int),
     'fill_jacobian_buffer_3d': (None, jac_buffer_3d, pos_buffer_3d, vertices, sz),
     'fill_matrix_3d': (None, arr(ndim=2), vertices, arr(dtype=C.c_uint8, ndim=1), arr(ndim=1), jac_buffer_3d, pos_buffer_3d, sz, sz, C.c_int, C.c_int),
@@ -196,14 +200,19 @@ def ensure_contiguous_aligned(arr):
     assert not DEBUG or (new_arr is arr), "Made copy while ensuring contiguous array"
     return new_arr
 
+# These are passed directly to scipy.LowLevelCallable, so we shouldn't wrap them in a function
+# that checks the numpy arrays
+numpy_wrapper_exceptions = ['self_potential_radial', 'self_field_dot_normal_radial']
+
 for (fun, (res, *args)) in backend_functions.items():
     libfun = getattr(backend_lib, fun)
-
-    def backend_check_numpy_requirements_wrapper(*args, _cfun_reference=libfun, _cfun_name=fun):
-        new_args = [ (ensure_contiguous_aligned(a) if isinstance(a, np.ndarray) else a) for a in args ]
-        return _cfun_reference(*new_args)
-
-    setattr(backend_lib, fun, backend_check_numpy_requirements_wrapper)
+
+    if fun not in numpy_wrapper_exceptions:
+        def backend_check_numpy_requirements_wrapper(*args, _cfun_reference=libfun, _cfun_name=fun):
+            new_args = [ (ensure_contiguous_aligned(a) if isinstance(a, np.ndarray) else a) for a in args ]
+            return _cfun_reference(*new_args)
+
+        setattr(backend_lib, fun, backend_check_numpy_requirements_wrapper)
 
     libfun.restype = res
     libfun.argtypes = args
@@ -552,6 +561,30 @@ def fill_jacobian_buffer_radial(vertices):
 
     return jac_buffer, pos_buffer
 
+def self_potential_radial(vertices):
+    assert vertices.shape == (4,3) and vertices.dtype == np.double
+    user_data = vertices.ctypes.data_as(C.c_void_p)
+    low_level = LowLevelCallable(backend_lib.self_potential_radial, user_data=user_data)
+    return quad(low_level, -1, 1, points=(0,), epsabs=1e-9, epsrel=1e-9)[0]
+
+class SelfFieldDotNormalRadialArgs(C.Structure):
+    _fields_ = [("line_points", C.POINTER(C.c_double)), ("K", C.c_double)]
+
+def self_field_dot_normal_radial(vertices, K):
+    assert vertices.shape == (4,3) and vertices.dtype == np.double
+    user_data = vertices.ctypes.data_as(C.c_void_p)
+
+    args = SelfFieldDotNormalRadialArgs()
+    args.K = float(K)
+    args.line_points = vertices.ctypes.data_as(dbl_p)
+
+    user_data = C.cast(C.pointer(args), vp)
+
+    low_level = LowLevelCallable(backend_lib.self_field_dot_normal_radial, user_data=user_data)
+    return quad(low_level, -1, 1, points=(0,), epsabs=1e-9, epsrel=1e-9)[0]
+
+
+
 def fill_matrix_radial(matrix, lines, excitation_types, excitation_values, jac_buffer, pos_buffer, start_index, end_index):
     N = len(lines)
     assert matrix.shape[0] == N and matrix.shape[1] == N and matrix.shape[0] == matrix.shape[1]

traceon/backend/kronrod.c

@@ -0,0 +1,77 @@
+#include <stdio.h>
+#include <math.h>
+
+#define G7_WEIGHTS_SIZE 4
+#define K15_WEIGHTS_SIZE 8
+
+const double G7_WEIGHTS[G7_WEIGHTS_SIZE] = {
+    0.417959183673469,
+    0.381830050505119,
+    0.279705391489277,
+    0.129484966168870,
+};
+
+const double G7_NODES[G7_WEIGHTS_SIZE] = {
+    0.000000000000000,
+    0.405845151377397,
+    0.741531185599394,
+    0.949107912342759,
+};
+
+const double K15_WEIGHTS[K15_WEIGHTS_SIZE] = {
+    0.209482141084728,
+    0.204432940075299,
+    0.190350578064785,
+    0.169004726639268,
+    0.140653259715526,
+    0.104790010322250,
+    0.063092092629979,
+    0.022935322010529,
+};
+
+const double K15_NODES[K15_WEIGHTS_SIZE] = {
+    0.000000000000000,
+    0.207784955007898,
+    0.405845151377397,
+    0.586087235467691,
+    0.741531185599394,
+    0.864864423359769,
+    0.949107912342759,
+    0.991455371120813,
+};
+
+double kronrod_adaptive(double (*f)(double, void*), double a, double b, void* args, double abs_tol, double rel_tol) {
+    double result = 0.0;
+    double current_start = a;
+    double current_end = b;
+
+    while (current_start < b) {
+        double c = (current_start + current_end) / 2.0;
+        double h = (current_end - current_start) / 2.0;
+
+        double g7 = h * G7_WEIGHTS[0] * f(c, args);
+        double k15 = h * K15_WEIGHTS[0] * f(c, args);
+
+        for (int i = 1; i < G7_WEIGHTS_SIZE; i++) {
+            double x = h * G7_NODES[i];
+            g7 += h * G7_WEIGHTS[i] * (f(c - x, args) + f(c + x, args));
+        }
+
+        for (int i = 1; i < K15_WEIGHTS_SIZE; i++) {
+            double x = h * K15_NODES[i];
+            k15 += h * K15_WEIGHTS[i] * (f(c - x, args) + f(c + x, args));
+        }
+
+        double error = fabs(k15 - g7);
+
+        if (error <= abs_tol || error <= rel_tol * fabs(k15)) {
+            result += k15;
+            current_start = current_end;
+            current_end = b;
+        } else {
+            current_end = c;
+        }
+    }
+
+    return result;
+}

traceon/backend/traceon-backend.c

@@ -6,9 +6,7 @@
 #include <stdbool.h>
 #include <time.h>
 
-#include <gsl/gsl_integration.h>
-#include <gsl/gsl_errno.h>
-
+#include "kronrod.c"
 #include "defs.c"
 #include "utilities_3d.c"
 #include "triangle_contribution.c"
@@ -1284,6 +1282,50 @@ enum ExcitationType{
 	MAGNETOSTATIC_POT=5,
 	MAGNETIZABLE=6};
 
+double self_potential_radial(double alpha, double line_points[4][3]) {
+
+	double *v1 = line_points[0];
+	double *v2 = line_points[2];
+	double *v3 = line_points[3];
+	double *v4 = line_points[1];
+
+	double pos[2], jac;
+	position_and_jacobian_radial(alpha, v1, v2, v3, v4, pos, &jac);
+
+	double target[2], jac2;
+	position_and_jacobian_radial(0, v1, v2, v3, v4, target, &jac2);
+
+	return jac*potential_radial_ring(target[0], target[1], pos[0], pos[1], NULL);
+}
+
+struct self_field_dot_normal_radial_args {
+	double (*line_points)[3];
+	double K;
+};
+
+double self_field_dot_normal_radial(double alpha, struct self_field_dot_normal_radial_args* args) {
+
+	double *v1 = args->line_points[0];
+	double *v2 = args->line_points[2];
+	double *v3 = args->line_points[3];
+	double *v4 = args->line_points[1];
+
+	double pos[2], jac;
+	position_and_jacobian_radial(alpha, v1, v2, v3, v4, pos, &jac);
+
+	double target[2], jac2;
+	position_and_jacobian_radial(0, v1, v2, v3, v4, target, &jac2);
+
+	double normal[2];
+	higher_order_normal_radial(0.0, v1, v2, v3, v4, normal);
+
+	struct {double *normal; double K;} cb_args = {normal, args->K};
+
+	return jac*field_dot_normal_radial(target[0], target[1], pos[0], pos[1], (void*) &cb_args);
+}
+
+
+
 
 struct self_voltage_radial_args {
 	double (*line_points)[3];
@@ -1310,59 +1352,6 @@ double self_voltage_radial(double alpha, void *args_p) {
 	return jac * args->cb_fun(args->target[0], args->target[1], pos[0], pos[1], &cb_args);
 }
 
-void fill_self_voltages_radial(double *matrix, 
-                        vertices_2d line_points,
-						uint8_t *excitation_types,
-						double *excitation_values,
-						size_t N_lines,
-						size_t N_matrix,
-                        int lines_range_start, 
-                        int lines_range_end) {
-
-	gsl_integration_workspace * w = gsl_integration_workspace_alloc(ADAPTIVE_MAX_ITERATION);
-
-	for(int i = lines_range_start; i <= lines_range_end; i++) {
-		double *v1 = &line_points[i][0][0];
-		double *v2 = &line_points[i][2][0];
-		double *v3 = &line_points[i][3][0];
-		double *v4 = &line_points[i][1][0];
-
-		double target[2], jac;
-		position_and_jacobian_radial(0.0, v1, v2, v3, v4, target, &jac);
-
-		double normal[2];
-		higher_order_normal_radial(0.0, v1, v2, v3, v4, normal);
-		//normal_2d(v1, v2, normal);
-
-		enum ExcitationType type_ = excitation_types[i];
-
-		struct self_voltage_radial_args integration_args = {
-			.target = target,
-			.line_points = &line_points[i][0],
-			.normal = normal,
-			.K = excitation_values[i],
-			.cb_fun = (type_ != DIELECTRIC && type_ != MAGNETIZABLE) ? potential_radial_ring : field_dot_normal_radial
-		};
-
-		gsl_function F;
-		F.function = &self_voltage_radial;
-		F.params = &integration_args;
-
-		double result, error;
-		double singular_points[3] = {-1, 0, 1};
-		gsl_integration_qagp(&F, singular_points, 3, 1e-9, 1e-9, ADAPTIVE_MAX_ITERATION, w, &result, &error);
-
-		if(type_ == DIELECTRIC || type_ == MAGNETIZABLE) {
-			matrix[N_matrix*i + i] = result - 1;
-		}
-		else {
-			matrix[N_matrix*i + i] = result;
-		}
-	}
-
-	gsl_integration_workspace_free(w);
-}
-
 EXPORT void fill_jacobian_buffer_radial(
 	jacobian_buffer_2d jacobian_buffer,
 	position_buffer_2d pos_buffer,
@@ -1399,7 +1388,6 @@ EXPORT void fill_matrix_radial(double *matrix,
                         int lines_range_start, 
                         int lines_range_end) {
 
-	gsl_set_error_handler_off();
 	assert(lines_range_start < N_lines && lines_range_end < N_lines);
 	assert(N_matrix >= N_lines);
 
@@ -1450,8 +1438,6 @@ EXPORT void fill_matrix_radial(double *matrix,
             exit(1);
 		}
 	}
-
-	fill_self_voltages_radial(matrix, line_points, excitation_types, excitation_values, N_lines, N_matrix, lines_range_start, lines_range_end);
 }
 
 EXPORT void fill_jacobian_buffer_3d(
@@ -1497,7 +1483,6 @@ EXPORT void fill_matrix_3d(double *restrict matrix,
                     int lines_range_start, 
                     int lines_range_end) {
 
-	gsl_set_error_handler_off();
 	assert(lines_range_start < N_lines && lines_range_end < N_lines);
 
     for (int i = lines_range_start; i <= lines_range_end; i++) {

traceon/backend/triangle_contribution.c

@@ -2,9 +2,6 @@
 #include <time.h>
 #include <math.h>
 
-#include <gsl/gsl_integration.h>
-#include <gsl/gsl_errno.h>
-
 // To efficiently compute the double integrals we define
 // a coordinate system as follows.
 // Let v0, v1, v2 be the vertices of the source triangle
@@ -86,16 +83,7 @@ potential_triangle(double v0[3], double v1[3], double v2[3], double target[3]) {
 
 	struct _normalized_triangle tri = {x0, y0, a,b,c,z};
 
-	gsl_function F;
-	F.function = _potential_integrand;
-	F.params = (void*) &tri;
-
-	double result, error;
-	gsl_integration_workspace *w = gsl_integration_workspace_alloc(ADAPTIVE_MAX_ITERATION);
-    gsl_integration_qag(&F, 0, c, 1e-9, 1e-9, ADAPTIVE_MAX_ITERATION, GSL_INTEG_GAUSS31, w, &result, &error);
-	gsl_integration_workspace_free(w);
-
-	return result;
+	return kronrod_adaptive(_potential_integrand, 0, c, (void*) &tri, 1e-9, 1e-9);
 }
 
 EXPORT double self_potential_triangle_v0(double v0[3], double v1[3], double v2[3]) {
@@ -200,16 +188,7 @@ flux_triangle(double v0[3], double v1[3], double v2[3], double target[3], double
 
 	struct _normalized_triangle tri = {x0, y0, a, b, c, z, new_normal};
 
-	gsl_function F;
-	F.function = _flux_integrand;
-	F.params = (void*) &tri;
-
-	double result, error;
-	gsl_integration_workspace *w = gsl_integration_workspace_alloc(ADAPTIVE_MAX_ITERATION);
-    gsl_integration_qag(&F, 0, c, 1e-9, 1e-9, ADAPTIVE_MAX_ITERATION, GSL_INTEG_GAUSS31, w, &result, &error);
-	gsl_integration_workspace_free(w);
-
-	return result;
+	return kronrod_adaptive(_flux_integrand, 0, c, (void*) &tri, 1e-9, 1e-9);
 }