calculate t1 for generalized gamma

2024-01-09 14:20:20 +01:00 · 2024-01-09 14:20:20 +01:00 · 50a811b7ec
commit 50a811b7ec
parent 73bdc71a83
5 changed files with 65 additions and 36 deletions
--- a/src/gui_qt/nmr/t1widget.py
+++ b/src/gui_qt/nmr/t1widget.py
@ -31,7 +31,7 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):
        self.t1calculator = RelaxationEvaluation()

        self.sd_parameter = []
-        self.sdmodels = [Debye, ColeCole, ColeDavidson, KWW, HavriliakNegami, LogGaussian]
+        self.sdmodels = [Debye, ColeCole, ColeDavidson, KWW, HavriliakNegami, LogGaussian, GGAlpha]
        for i in self.sdmodels:
            self.specdens_combobox.addItem(i.name)
        self.specdens_combobox.currentIndexChanged.connect(self.update_specdens)
@ -51,8 +51,14 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):
        self.conv_y = QT1Widget.time_conversion[self.t1_combobox.currentIndex()]

        self.minimum = (1, np.inf)
-        self.min_pos = PlotItem(x=np.array([]), y=np.array([]),
-                                symbol='+', symbolBrush=mkBrush(color='r'), symbolPen=mkPen(color='r'), symbolSize=14)
+        self.min_pos = PlotItem(
+            x=np.array([]),
+            y=np.array([]),
+            symbol='+',
+            symbolBrush=mkBrush(color='r'),
+            symbolPen=mkPen(color='r'),
+            symbolSize=14,
+        )
        self.parabola = PlotItem(x=np.array([]), y=np.array([]))

        self.lineEdit_2.setValidator(QtGui.QDoubleValidator())
@ -83,10 +89,10 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):
        right_b = min(np.argmin(y)+3, len(x)-1)

        self.lineEdit_2.blockSignals(True)
-        self.lineEdit_2.setText('{:.2f}'.format(x[left_b]))
+        self.lineEdit_2.setText(f'{x[left_b]:.2f}')
        self.lineEdit_2.blockSignals(False)
        self.lineEdit_3.blockSignals(True)
-        self.lineEdit_3.setText('{:.2f}'.format(x[right_b]))
+        self.lineEdit_3.setText(f'{x[right_b]:.2f}')
        self.lineEdit_3.blockSignals(False)

        self.t1calculator.set_data(x, y)
@ -110,6 +116,7 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):

        if self.sdmodels[idx].parameter is not None:
            for name in self.sdmodels[idx].parameter:
+                print(name)
                _temp = FormWidget(parent=self, name=name, fixable=True)
                _temp.value = 1
                _temp.setChecked(True)
@ -133,7 +140,7 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):
        try:
            for i, v, in enumerate(values):
                self.sd_parameter[i].blockSignals(True)
-                self.sd_parameter[i].value = '{:.3g}'.format(round(v, 3))
+                self.sd_parameter[i].value = f'{v:.3g}'
                self.sd_parameter[i].blockSignals(False)
        except IndexError:
            pass
@ -219,7 +226,7 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):
        self.update_model()

    @QtCore.pyqtSlot(int, name='on_interpol_combobox_currentIndexChanged')
-    def determine_minimum(self, idx):
+    def determine_minimum(self, idx: int):
        if idx == 0:
            self.checkBox_interpol.setChecked(False)
            self.checkBox_interpol.hide()
@ -229,9 +236,10 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):
            self.checkBox_interpol.show()
            self.frame.show()
            try:
-                m, i_func = self.t1calculator.calculate_t1_min(interpolate=idx,
-                                                               trange=(float(self.lineEdit_2.text()),
-                                                                       float(self.lineEdit_3.text())))
+                m, i_func = self.t1calculator.calculate_t1_min(
+                    interpolate=idx,
+                    trange=(float(self.lineEdit_2.text()), float(self.lineEdit_3.text())),
+                )
            except ValueError:
                m, i_func = self.t1calculator.calculate_t1_min(interpolate=None)

@ -273,11 +281,13 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):
            return

        with busy_cursor():
-            calc_stretching, mini = self.t1calculator.get_increase(height=self.minimum[1],
-                                                                   idx=var_idx, mode=notfix,
-                                                                   omega=2*np.pi*self.frequency,
-                                                                   dist_parameter=sd_args, prefactor=cp_args,
-                                                                   coupling_kwargs=cp_kwargs)
+            calc_stretching, mini = self.t1calculator.get_increase(
+                height=self.minimum[1],
+                idx=var_idx, mode=notfix,
+                omega=2*np.pi*self.frequency,
+                dist_parameter=sd_args, prefactor=cp_args,
+                coupling_kwargs=cp_kwargs
+            )

        self.label_t1min.setText(f'{mini:.4g} s')

@ -292,9 +302,13 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):
        sd_args, _ = self.get_sd_values()
        cp_args, cp_kwargs, _ = self.get_cp_values()
        tau_mode = ['fit', 'peak', 'mean', 'logmean'][self.tau_combox.currentIndex()]
-        corr, opts = self.t1calculator.correlation_from_t1(omega=2*np.pi*self.frequency, dist_parameter=sd_args,
-                                                           coupling_param=cp_args, coupling_kwargs=cp_kwargs,
-                                                           mode=tau_mode, interpolate=self.checkBox_interpol.isChecked())
+        corr, opts = self.t1calculator.correlation_from_t1(
+            omega=2*np.pi*self.frequency,
+            dist_parameter=sd_args,
+            coupling_param=cp_args, coupling_kwargs=cp_kwargs,
+            mode=tau_mode,
+            interpolate=self.checkBox_interpol.isChecked()
+        )

        name = self.name + '-' + str(self.t1calculator) + '('
        name += ','.join([f'{a:.3g}' for a in sd_args])
@ -332,4 +346,4 @@ class QT1Widget(QtWidgets.QDialog, Ui_t1dialog):

    @QtCore.pyqtSlot(int, name='on_graph_checkbox_stateChanged')
    def changed_state(self, checked):
-        self.graph_combobox.setEnabled(checked != QtCore.Qt.Checked)
+        self.graph_combobox.setEnabled(checked != QtCore.Qt.CheckState.Checked)
--- a/src/nmreval/distributions/init.py
+++ b/src/nmreval/distributions/init.py
@ -26,3 +26,4 @@ from .coledavidson import ColeDavidson
 from .debye import Debye
 from .kww import KWW
 from .loggaussian import LogGaussian
+from .gengamma import GGAlpha
--- a/src/nmreval/distributions/gengamma.py
+++ b/src/nmreval/distributions/gengamma.py
@ -20,10 +20,10 @@ class AbstractGG(Distribution, ABC):

    @classmethod
    def correlation(cls, t, tau0, *args):
-        tt = np.asanyarray(t)
+        tt = np.atleast_1d(t)
        taus, ln_tau = AbstractGG._prepare_integration(tau0)
        g_tau = cls.distribution(taus, tau0, *args)
-        ret_val = np.array([simpson(np.exp(-t_i/taus) * g_tau, ln_tau) for t_i in tt])
+        ret_val = np.array([simpson(np.exp(-t_i/taus) * g_tau, ln_tau) for t_i in tt]).squeeze()

        return ret_val

@ -32,11 +32,11 @@ class AbstractGG(Distribution, ABC):
        r"""
        Calculate spectral density \int G(ln(tau) tau/(1+(w*tau)^2) dln(tau)
        """
-        w = np.asanyarray(omega)
+        w = np.atleast_1d(omega)
        taus, ln_tau = AbstractGG._prepare_integration(tau0)
        g_tau = cls.distribution(taus, tau0, *args)

-        ret_val = np.array([simpson(g_tau / (1 - 1j*w_i*taus), ln_tau) for w_i in w])
+        ret_val = np.array([simpson(g_tau / (1 - 1j*w_i*taus), ln_tau) for w_i in w]).squeeze()

        return ret_val

@ -45,17 +45,23 @@ class AbstractGG(Distribution, ABC):
        r"""
        Calculate spectral density \int G(ln(tau) tau/(1+(w*tau)^2) dln(tau)
        """
-        w = np.asanyarray(omega)
-        taus, ln_tau = AbstractGG._prepare_integration(tau0)
-        g_tau = cls.distribution(taus, tau0, *args)
+        w = np.atleast_1d(omega)
+        _t = np.atleast_1d(tau0)
+        ret_val = np.zeros((w.size, _t.size))

-        ret_val = np.array([simpson(g_tau * taus / (1 + (w_i*taus)**2), ln_tau) for w_i in w])
+        for i, tau_i in enumerate(_t):
+            taus, ln_tau = AbstractGG._prepare_integration(tau_i)
+            g_tau = cls.distribution(taus, tau_i, *args)

-        return ret_val
+            ret_val[:, i] = np.array([simpson(g_tau * taus / (1 + (w_i*taus)**2), ln_tau) for w_i in w]).squeeze()
+
+        return ret_val.squeeze()

    @staticmethod
    def _prepare_integration(
-            tau0: float, limits: tuple[int, int] = (20, 20), num_steps: int = 4001
+            tau0: float,
+            limits: tuple[int, int] = (20, 20),
+            num_steps: int = 4001,
    ) -> tuple[np.ndarray, np.ndarray]:
        """
        Create array of correlation times for integration over ln(tau)
@ -77,8 +83,8 @@ class AbstractGG(Distribution, ABC):

 # noinspection PyMethodOverriding
 class GGAlpha(AbstractGG):
-    name = r'General \Gamma (\alpha)'
-    parameter = [r'\tau', r'\alpha', r'\beta']
+    name = r'General Gamma (alpha)'
+    parameter = [r'\alpha', r'\beta']

    @staticmethod
    def distribution(taus: float | np.ndarray, tau: float, alpha: float, beta: float) -> float | np.ndarray:
@ -92,8 +98,8 @@ class GGAlpha(AbstractGG):

 # noinspection PyMethodOverriding
 class GGAlphaEW(AbstractGG):
-    name = r'General \Gamma (\alpha + EW)'
-    parameter = [r'\tau', r'\alpha', r'\beta', r'\sigma', r'\gamma']
+    name = r'General Gamma (alpha + EW)'
+    parameter = [r'\alpha', r'\beta', r'\sigma', r'\gamma']

    @staticmethod
    def distribution(tau: float | np.ndarray, tau0: float,
@ -117,8 +123,8 @@ class GGAlphaEW(AbstractGG):

 # noinspection PyMethodOverriding
 class GGBeta(AbstractGG):
-    name = r'General \Gamma (\beta)'
-    parameter = [r'\tau', 'a', 'b']
+    name = r'General Gamma (beta)'
+    parameter = ['a', 'b']

    @staticmethod
    def distribution(tau: float | np.ndarray, tau0: float, a: float, b: float) -> float | np.ndarray:
--- a/src/nmreval/models/fieldcycling.py
+++ b/src/nmreval/models/fieldcycling.py
@ -2,6 +2,7 @@ import numpy as np

 from ..distributions import *
 from ..distributions.energy import EnergyBarriers
+from ..distributions.gengamma import GGAlpha
 from ..distributions.intermolecular import FFHS
 from ..nmr.relaxation import Relaxation
 from ..utils.constants import gamma
@ -82,6 +83,13 @@ class FFHSFC(_AbstractFC):
    relax = Relaxation(distribution=FFHS)


+class GGAFC(_AbstractFC):
+    name = 'GG(alpha)'
+    params = _AbstractFC.params + [r'\alpha', r'\beta']
+    bounds = _AbstractFC.bounds + [(None, None), (None, None)]
+    relax = Relaxation(distribution=GGAlpha)
+
+
 class EnergyFC(_AbstractFC):
    name = 'Energy distribution'
    params = ['C', 'T'] + EnergyBarriers.parameter
--- a/src/nmreval/nmr/relaxation.py
+++ b/src/nmreval/nmr/relaxation.py
@ -525,7 +525,7 @@ class RelaxationEvaluation(Relaxation):
            dist_parameter: tuple | list = None,
            prefactor: tuple | list | float = None,
            coupling_kwargs: dict = None,
-    ) -> None:
+    ) -> tuple[float, float] :
        """
        Determine a single parameter from a T1 minimum.
        It replaces the previously set value.