From 22217f05447fdc4818c0d2d1d1d992a86ef619d6 Mon Sep 17 00:00:00 2001 From: dominik Date: Wed, 23 Mar 2022 15:32:56 +0100 Subject: [PATCH] fixed kww calculation with array of betas --- nmreval/gui_qt/main/mainwindow.py | 4 ++-- nmreval/gui_qt/nmr/t1_from_tau.py | 8 +++---- nmreval/math/kww.py | 36 ++++++++++++++++++++++++++++ nmreval/nmr/relaxation.py | 39 +++++++++++++++++++++++++++++++ 4 files changed, 81 insertions(+), 6 deletions(-) diff --git a/nmreval/gui_qt/main/mainwindow.py b/nmreval/gui_qt/main/mainwindow.py index b0316a7..d5694be 100644 --- a/nmreval/gui_qt/main/mainwindow.py +++ b/nmreval/gui_qt/main/mainwindow.py @@ -514,7 +514,7 @@ class NMRMainWindow(QtWidgets.QMainWindow, Ui_BaseWindow): def _select_t1tauwidget(self, onoff: bool, pick_required: bool, block_window: bool): if onoff: # tau from t1 if self.current_graph_widget is None: - return + return pick_required, block_window idx = self.tabWidget.indexOf(self.t1tauwidget) if len(self.current_graph_widget) != 1: @@ -522,7 +522,7 @@ class NMRMainWindow(QtWidgets.QMainWindow, Ui_BaseWindow): 'Only one T1 curve can be evaluated at once') self.tabWidget.removeTab(idx) self.tabWidget.setCurrentIndex(0) - return + return pick_required, block_window self.tabWidget.setTabText(idx, 'T1 min (%s)' % {self.current_graph_widget.title}) diff --git a/nmreval/gui_qt/nmr/t1_from_tau.py b/nmreval/gui_qt/nmr/t1_from_tau.py index 51355e7..71fed17 100644 --- a/nmreval/gui_qt/nmr/t1_from_tau.py +++ b/nmreval/gui_qt/nmr/t1_from_tau.py @@ -1,7 +1,7 @@ from typing import List, Tuple from ...nmr.coupling import * -from ...distributions import ColeCole, ColeDavidson, HavriliakNegami +from ...distributions import ColeCole, ColeDavidson, HavriliakNegami, KWW from ...utils import pi from ...utils.text import convert @@ -19,7 +19,7 @@ class QRelaxCalc(QtWidgets.QDialog, Ui_Dialog): self.graphs = {} - self.specdens = [ColeCole, ColeDavidson, HavriliakNegami] + self.specdens = [ColeCole, ColeDavidson, HavriliakNegami, KWW] self.coupling = [Quadrupolar, Czjzek, HomoDipolar] self.tau_parameter = [] @@ -162,9 +162,9 @@ class QRelaxCalc(QtWidgets.QDialog, Ui_Dialog): if isinstance(p, Widget): parameter.append(p.value) elif isinstance(p, SelectionWidget): - kwargs[p.argname] = p.value + kwargs.update(p.value) else: - raise TypeError('WTF?: Unknown widget', p) + raise TypeError('WTF: Unknown widget', p) dic['cp_param'] = (parameter, kwargs) diff --git a/nmreval/math/kww.py b/nmreval/math/kww.py index dd3e115..48bf395 100644 --- a/nmreval/math/kww.py +++ b/nmreval/math/kww.py @@ -204,6 +204,24 @@ def _kwws_lim_hig(b): def kww_cos(w, tau, beta): + beta = np.asanyarray(beta) + if beta.ndim == 0: + return _kww_cos_single(w, tau, beta) + + else: + tau = np.asanyarray(tau) + if tau.size != beta.size: + raise ValueError('KWW mismatched number of beta and tau') + else: + beta.reshape(tau.shape) + ret_val = np.zeros((w.size, tau.size)) + for i, (tau_i, beta_i) in enumerate(zip(tau.ravel(), beta.ravel())): + ret_val[:, i] = _kww_cos_single(w, tau_i, beta_i) + + return ret_val.reshape(-1, *tau.shape) + + +def _kww_cos_single(w, tau, beta): """ Args: @@ -246,6 +264,24 @@ def kww_cos(w, tau, beta): # \int_0^\infty dt sin(w*t) exp(-(t/tau)^beta) def kww_sin(w, tau, beta): + beta = np.asanyarray(beta) + if beta.ndim == 0: + return _kww_sin_single(w, tau, beta) + + else: + tau = np.asanyarray(tau) + if tau.size != beta.size: + raise ValueError('KWW mismatched number of beta and tau') + else: + beta.reshape(tau.shape) + ret_val = np.zeros((w.size, tau.size)) + for i, (tau_i, beta_i) in enumerate(zip(tau.ravel(), beta.ravel())): + ret_val[:, i] = _kww_sin_single(w, tau_i, beta_i) + + return ret_val.reshape(-1, *tau.shape) + + +def _kww_sin_single(w, tau, beta): # check input data if not (0.1 <= beta <= 2.): raise ValueError('KWW beta is not inside range 0.1-2') diff --git a/nmreval/nmr/relaxation.py b/nmreval/nmr/relaxation.py index c7e403d..cdf6f82 100755 --- a/nmreval/nmr/relaxation.py +++ b/nmreval/nmr/relaxation.py @@ -104,6 +104,8 @@ class Relaxation: return self.t1_bpp(omega, tau, *specdens_args, **kwargs) elif mode == 'dipolar': return self.t1_dipolar(omega, tau, *specdens_args, **kwargs) + else: + return self.t1_csa(omega, tau, *specdens_args, **kwargs) def t1_dipolar(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any, inverse: bool = True, prefactor: float = None, omega_coup: ArrayLike = None, @@ -302,6 +304,43 @@ class Relaxation: else: return rate + def t2(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any, + mode: str = 'bpp', **kwargs) -> Union[np.ndarray, float]: + r""" + Convenience function + + Args: + omega (array-like): Frequency in 1/s (not Hz) + tau (array-like): Correlation times in s + *specdens_args: Additional parameter for spectral density. + If given this will replace previously set arguments during calculation + mode (str, {`bpp`, `dipolar`, `csa`}): Type of relaxation + + - `bpp` : Homonuclear dipolar/quadrupolar interaction :func:`(see here) ` + - `dipolar` : Heteronuclear dipolar interaction :func:`(see here) ` + - `csa` : Chemical shift interaction :func:`(see here) ` + + Default is `bpp` + + **kwargs: + + Returns: + float or ndarray: + A k by n array where k is length of ``omega`` and n is length of ``tau``. + If both are of length 1 a number is returned instead of array. + + """ + if mode not in ['bpp', 'dipolar', 'csa']: + raise ValueError(f'Unknown mode {mode} not `bpp`, `dipolar`, `csa`.') + + # defaults to BPP + if mode == 'bpp': + return self.t2_bpp(omega, tau, *specdens_args, **kwargs) + elif mode == 'dipolar': + return self.t2_dipolar(omega, tau, *specdens_args, **kwargs) + else: + return self.t2_csa(omega, tau, *specdens_args, **kwargs) + def t2_bpp(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any, inverse: bool = True, prefactor: float = None) -> Union[np.ndarray, float]: r"""Calculate T2 under homonuclear dipolar coupling or quadrupolar coupling.