From ead30d127a92968bb3efba51b29471526b19172e Mon Sep 17 00:00:00 2001
From: Dominik Demuth <dominik.demuth@physik.tu-darmstadt.de>
Date: Mon, 10 Apr 2023 19:33:23 +0200
Subject: [PATCH] add C function for FFHS integration

---
 Makefile                                    | 22 +++++++-------
 src/nmreval/clib/integrate.c                | 15 ++++++++++
 src/nmreval/distributions/intermolecular.py | 32 ++++++++++++++++++++-
 src/nmreval/fit/minimizer.py                |  1 -
 4 files changed, 57 insertions(+), 13 deletions(-)
 create mode 100644 src/nmreval/clib/integrate.c

diff --git a/Makefile b/Makefile
index 8c9b2d5..5f03e4f 100755
--- a/Makefile
+++ b/Makefile
@@ -2,37 +2,37 @@
 
 #binaries
 PYUIC = pyuic5
-PYRCC = pyrcc5
+# PYRCC = pyrcc5
+
+CC = gcc
+CFLAGS = -O2 -fPIC
+LDFLAGs = -shared
+
 
 #Directory with ui files
-RESOURCE_DIR = resources/_ui
+RESOURCE_DIR = src/resources/_ui
 
 #Directory for compiled resources
-COMPILED_DIR = nmreval/gui_qt/_py
+PYQT_DIR = src/gui_qt/_py
 
 #UI files to compile, uses every *.ui found in RESOURCE_DIR
 UI_FILES = $(foreach dir, $(RESOURCE_DIR), $(notdir $(wildcard $(dir)/*.ui)))
-COMPILED_UI = $(UI_FILES:%.ui=$(COMPILED_DIR)/%.py)
+PYQT_UI = $(UI_FILES:%.ui=$(PYQT_DIR)/%.py)
 
 SVG_FILES = $(foreach dir, $(RCC_DIR), $(notdir $(wildcard $(dir)/*.svg)))
 PNG_FILES = $(SVG_FILES:%.svg=$(RCC_DIR)/%.png)
 
 all : ui
  
-ui : $(COMPILED_UI)
+ui : $(PYQT_UI)
 
 rcc:  $(PNG_FILES)
 
- 
 $(COMPILED_DIR)/%.py : $(RESOURCE_DIR)/%.ui
 	$(PYUIC) $< -o $@
-# 	replace import of ressource to correct path
-#	@sed -i s/images_rc/nmrevalqt.$(COMPILED_DIR).images_rc/g $@
-#	@sed -i /images_rc/d $@
-	
+
 $(RCC_DIR)/%.png : $(RCC_DIR)/%.svg
 	convert -background none $< $@
-	$(PYRCC) $(RCC_DIR)/images.qrc -o $(COMPILED_DIR)/images_rc.py
 
 clean:
 	find . -name '*.pyc' -exec rm -f {} +
diff --git a/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c
new file mode 100644
index 0000000..3cec66e
--- /dev/null
+++ b/src/nmreval/clib/integrate.c
@@ -0,0 +1,15 @@
+/* integrands used in quadrature integration with scipy's LowLevelCallables */
+
+/* integrand for distributions/intermolecular/FFHS */
+double ffhs_sd(double x, void *user_data) {
+	double *c = (double *)user_data;
+	double tau = c[1];
+	double omega = c[0];
+
+	double u = x*x;
+	double res = u*u * c[1];
+
+	res /= 81. + 9.*u - 2.*u*u + u*u*u;
+	res /= u*u + omega*omega * tau*tau;
+	return res;
+}
diff --git a/src/nmreval/distributions/intermolecular.py b/src/nmreval/distributions/intermolecular.py
index 1607361..c1b6354 100644
--- a/src/nmreval/distributions/intermolecular.py
+++ b/src/nmreval/distributions/intermolecular.py
@@ -1,8 +1,24 @@
+import ctypes
+import os.path
+from pathlib import Path
+
 import numpy as np
+from scipy import LowLevelCallable
 from scipy.integrate import quad
 
 from .base import Distribution
 
+try:
+    print(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
+    lib = ctypes.CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
+    lib.ffhs_sd.restype = ctypes.c_double
+    lib.ffhs_sd.argtypes = (ctypes.c_double, ctypes.c_void_p)
+    HAS_C_FUNCS = True
+    print('Use C functions')
+except OSError:
+    HAS_C_FUNCS = False
+    print('Use python functions')
+
 
 class FFHS(Distribution):
     name = 'Intermolecular (FFHS)'
@@ -24,7 +40,7 @@ class FFHS(Distribution):
         return ret_val
 
     @staticmethod
-    def specdens(omega, tau0, *args):
+    def specdens_py(omega, tau0, *args):
         def integrand(u, o, tau0):
             return u**4 * tau0 / (81 + 9*u**2 - 2*u**4 + u**6) / (u**4 + (o*tau0)**2)
             # return FFHS.distribution(u, tau0) * u / (1+o**2 * u**2)
@@ -33,6 +49,17 @@ class FFHS(Distribution):
 
         return ret_val * 54 / np.pi
 
+    @staticmethod
+    def specdens_c(omega, tau0, *args):
+        res = []
+        for o in omega:
+            c = (ctypes.c_double * 2)(o, tau0)
+            user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
+            func = LowLevelCallable(lib.ffhs_sd, user_data)
+            res.append(quad(func, 0, np.infty)[0])
+
+        return np.array(res) * 54 / np.pi
+
     @staticmethod
     def susceptibility(omega, tau0, *args):
         def integrand_real(u, o, tau0):
@@ -48,6 +75,9 @@ class FFHS(Distribution):
         return ret_val
 
 
+FFHS.specdens = FFHS.specdens_c if HAS_C_FUNCS else FFHS.specdens_py
+
+
 # class Bessel(Distribution):
 #     name = 'Intermolecular (Bessel)'
 #     parameter = None
diff --git a/src/nmreval/fit/minimizer.py b/src/nmreval/fit/minimizer.py
index 2655e4a..57e1b56 100644
--- a/src/nmreval/fit/minimizer.py
+++ b/src/nmreval/fit/minimizer.py
@@ -1,4 +1,3 @@
-import time
 import warnings
 from itertools import product