Skip to content

Commit

Permalink
Merge pull request #37189 from robertapplin/36502-teixeira-water-in-i…
Browse files Browse the repository at this point in the history 
…qtfit

Add TeixeiraWaterIqt to IqtFit
  • Loading branch information
@SilkeSchomann
SilkeSchomann committed Apr 26, 2024
2 parents f552f6d + 214f363 commit 08c54b8
Show file tree
Hide file tree
Showing 15 changed files with 119 additions and 82 deletions.
7 changes: 4 additions & 3 deletions Framework/PythonInterface/mantid/api/src/FitFunctions/IFunctionAdapter.cpp
View file
Expand Up @@ -203,11 +203,12 @@ void IFunctionAdapter::setAttributePythonValue(IFunction &self, const std::strin
* @param attr An attribute object
*/
void IFunctionAdapter::setAttribute(const std::string &attName, const Attribute &attr) {
try {
auto self = getSelf();
if (typeHasAttribute(self, "setAttributeValue")) {
object value = object(handle<>(getAttributeValue(*this, attr)));
callMethod<void, std::string, object>(getSelf(), "setAttributeValue", attName, value);
callMethod<void, std::string, object>(self, "setAttributeValue", attName, value);
storeAttributeValue(attName, attr);
} catch (UndefinedAttributeError &) {
} else {
IFunction::setAttribute(attName, attr);
}
}
Expand Down
41 changes: 7 additions & 34 deletions Framework/PythonInterface/plugins/functions/TeixeiraWaterIqt.py
View file
Expand Up @@ -5,16 +5,10 @@
# & Institut Laue - Langevin
# SPDX - License - Identifier: GPL - 3.0 +
# pylint: disable=no-init,invalid-name

"""
@author Spencer Howells, ISIS
@date December 05, 2013
"""

import numpy as np

from mantid.api import IFunction1D, FunctionFactory
from scipy.special import spherical_jn # bessel
from scipy.special import spherical_jn


class TeixeiraWaterIqt(IFunction1D):
Expand All @@ -32,46 +26,25 @@ def init(self):
self.declareAttribute("Q", 0.4)
self.declareAttribute("a", 0.98)

def setAttributeValue(self, name, value):
if name == "Q":
self.Q = value
if name == "a":
self.radius = value

def function1D(self, xvals):
amp = self.getParameterValue("Amp")
tau1 = self.getParameterValue("Tau1")
gamma = self.getParameterValue("Gamma")

q_value = self.getAttributeValue("Q")
radius = self.getAttributeValue("a")

xvals = np.array(xvals)

qr = np.array(self.Q * self.radius)
qr = np.array(q_value * radius)
j0 = spherical_jn(0, qr)
j1 = spherical_jn(1, qr)
j2 = spherical_jn(2, qr)
with np.errstate(divide="ignore"):
rotational = j0 * j0 + 3 * j1 * j1 * np.exp(-xvals / (3 * tau1)) + 5 * j2 * j2 * np.exp(-xvals / tau1)
rotational = np.square(j0) + 3 * np.square(j1) * np.exp(-xvals / (3 * tau1)) + 5 * np.square(j2) * np.exp(-xvals / tau1)
translational = np.exp(-gamma * xvals)
iqt = amp * rotational * translational
return iqt

def functionDeriv1D(self, xvals, jacobian):
amp = self.getParameterValue("Amp")
tau1 = self.getParameterValue("Tau1")
gamma = self.getParameterValue("Gamma")

qr = np.array(self.Q * self.radius)
j0 = spherical_jn(0, qr)
j1 = spherical_jn(1, qr)
j2 = spherical_jn(2, qr)
with np.errstate(divide="ignore"):
for i, x in enumerate(xvals, start=0):
rotational = j0 * j0 + 3 * j1 * j1 * np.exp(-x / (3 * tau1)) + 5 * j2 * j2 * np.exp(-x / tau1)
translational = np.exp(-gamma * x)
partial_tau = (x / np.square(tau1)) * (j1 * j1 * np.exp(-x / (3 * tau1)) + 5 * j2 * j2 * np.exp(-x / tau1))
jacobian.set(i, 0, rotational * translational)
jacobian.set(i, 1, amp * rotational * partial_tau)
jacobian.set(i, 2, -x * amp * rotational * translational)


# Required to have Mantid recognise the new function
FunctionFactory.subscribe(TeixeiraWaterIqt)
2 changes: 1 addition & 1 deletion docs/source/release/v6.10.0/Inelastic/New_features/36502.rst
View file
@@ -1 +1 @@
- Added new `TeixeiraWaterIqt` fitting function, to fit linewidth and molecular residence time in intermediate scattering functions with Teixeira's model for water.
- Added new `TeixeiraWaterIqt` fitting function, to fit linewidth and molecular residence time in intermediate scattering functions with Teixeira's model for water. This is now available in the IqtFit tab of the :ref:`QENS Fitting <interface-inelastic-qens-fitting>` interface.
2 changes: 1 addition & 1 deletion qt/python/mantidqt/mantidqt/_common.sip
View file
Expand Up @@ -351,7 +351,7 @@ public:

UserSubWindow *createSubWindow(
const QString &interface_name,
QWidget *parent = nullptr);
QWidget *parent = nullptr) throw(std::exception);
QStringList getUserSubWindowKeys() const;

void showHelpPage(const QString &url=QString());
Expand Down
1 change: 1 addition & 0 deletions qt/scientific_interfaces/Inelastic/QENSFitting/FitTabConstants.h
View file
Expand Up @@ -23,6 +23,7 @@ static const auto FUNCTION_STRINGS =
{"Lorentzian", "L"},
{"StretchedExpFT", "SFT"},
{"TeixeiraWater", "TxWater"},
{"TeixeiraWaterIqt", "TxWater"},
{"TeixeiraWaterSQE", "TxWater"},
{"FickDiffusionSQE", "FickDiff"},
{"ChudleyElliotSQE", "ChudElliot"},
Expand Down
2 changes: 0 additions & 2 deletions qt/scientific_interfaces/Inelastic/QENSFitting/FunctionBrowser/ConvFunctionTemplateModel.cpp
View file
Expand Up @@ -159,8 +159,6 @@ void ConvFunctionTemplateModel::checkSingleFunction(const IFunction_sptr &fun, b
}
}

void ConvFunctionTemplateModel::setQValues(const std::vector<double> &qValues) { m_qValues = qValues; }

void ConvFunctionTemplateModel::addFunction(std::string const &prefix, std::string const &funStr) {
if (!prefix.empty())
throw std::runtime_error("Function doesn't have member function with prefix " + prefix);
Expand Down
5 changes: 2 additions & 3 deletions qt/scientific_interfaces/Inelastic/QENSFitting/FunctionBrowser/ConvFunctionTemplateModel.h
View file
Expand Up @@ -41,7 +41,6 @@ class MANTIDQT_INELASTIC_DLL ConvFunctionTemplateModel : public MultiFunctionTem
std::map<std::size_t, int> getSubTypes() const override;
std::string setBackgroundA0(double value) override;
void setResolution(const std::vector<std::pair<std::string, size_t>> &fitResolutions) override;
void setQValues(const std::vector<double> &qValues) override;

bool hasDeltaFunction() const;
bool hasTempCorrection() const;
Expand All @@ -55,14 +54,15 @@ class MANTIDQT_INELASTIC_DLL ConvFunctionTemplateModel : public MultiFunctionTem
void applyParameterFunction(const std::function<void(ParamID)> &paramFun) const override;

void clearData();
void setModel();
void setModel() override;

std::optional<std::string> getLor1Prefix() const;
std::optional<std::string> getLor2Prefix() const;
std::optional<std::string> getFitTypePrefix() const;
std::optional<std::string> getDeltaPrefix() const;
std::optional<std::string> getBackgroundPrefix() const;

std::string buildFunctionString(int const domainIndex) const override { return ""; }
std::string buildLorentzianFunctionString() const;
std::string buildTeixeiraFunctionString() const;
std::string buildFickFunctionString() const;
Expand Down Expand Up @@ -96,7 +96,6 @@ class MANTIDQT_INELASTIC_DLL ConvFunctionTemplateModel : public MultiFunctionTem
BackgroundSubType m_backgroundSubtype;

std::vector<std::pair<std::string, size_t>> m_fitResolutions;
std::vector<double> m_qValues;
bool m_isQDependentFunction = false;
};

Expand Down
2 changes: 2 additions & 0 deletions qt/scientific_interfaces/Inelastic/QENSFitting/FunctionBrowser/FitTypes.cpp
View file
Expand Up @@ -25,6 +25,8 @@ std::map<IqtTypes::FitType, TemplateSubTypeDescriptor> TemplateSubTypeImpl<IqtTy
{IqtTypes::FitType::None, {"None", "", {ParamID::NONE, ParamID::NONE}}},
{IqtTypes::FitType::StretchExponential,
{"Stretch Exponential", "StretchExp", {ParamID::STRETCH_HEIGHT, ParamID::STRETCH_STRETCHING}}},
{IqtTypes::FitType::TeixeiraWaterIqt,
{"Teixeira Water Iqt", "TeixeiraWaterIqt", {ParamID::TWI_AMPLITUDE, ParamID::TWI_GAMMA}}},
};

template <>
Expand Down
2 changes: 1 addition & 1 deletion qt/scientific_interfaces/Inelastic/QENSFitting/FunctionBrowser/FitTypes.h
View file
Expand Up @@ -23,7 +23,7 @@ enum class ExponentialType {
TwoExponentials,
};

enum class FitType { None, StretchExponential };
enum class FitType { None, StretchExponential, TeixeiraWaterIqt };

enum class BackgroundType { None, Flat };

Expand Down
74 changes: 49 additions & 25 deletions qt/scientific_interfaces/Inelastic/QENSFitting/FunctionBrowser/IqtFunctionTemplateModel.cpp
View file
Expand Up @@ -64,11 +64,8 @@ void IqtFunctionTemplateModel::clearData() {
}

void IqtFunctionTemplateModel::setModel() {
m_model->setFunctionString(buildFunctionString());
m_model->setGlobalParameters(makeGlobalList());

MultiFunctionTemplateModel::setModel();
tieIntensities();

estimateFunctionParameters();
}

Expand All @@ -82,6 +79,8 @@ void IqtFunctionTemplateModel::setFunction(IFunction_sptr fun) {
m_exponentialType = ExponentialType::OneExponential;
} else if (name == "StretchExp") {
m_fitType = FitType::StretchExponential;
} else if (name == "TeixeiraWaterIqt") {
m_fitType = FitType::TeixeiraWaterIqt;
} else if (name == "FlatBackground") {
m_backgroundType = BackgroundType::Flat;
} else {
Expand All @@ -91,7 +90,7 @@ void IqtFunctionTemplateModel::setFunction(IFunction_sptr fun) {
return;
}
bool areExponentialsSet = false;
bool isStretchSet = false;
bool isFitTypeSet = false;
bool isBackgroundSet = false;
for (size_t i = 0; i < fun->nFunctions(); ++i) {
auto f = fun->getFunction(i);
Expand All @@ -107,19 +106,23 @@ void IqtFunctionTemplateModel::setFunction(IFunction_sptr fun) {
areExponentialsSet = true;
}
} else if (name == "StretchExp") {
if (isStretchSet) {
if (isFitTypeSet) {
throw std::runtime_error("Function has wrong structure.");
}
m_fitType = FitType::StretchExponential;
areExponentialsSet = true;
isStretchSet = true;
isFitTypeSet = true;
} else if (name == "TeixeiraWaterIqt") {
m_fitType = FitType::TeixeiraWaterIqt;
areExponentialsSet = true;
isFitTypeSet = true;
} else if (name == "FlatBackground") {
if (isBackgroundSet) {
throw std::runtime_error("Function has wrong structure.");
}
m_backgroundType = BackgroundType::Flat;
areExponentialsSet = true;
isStretchSet = true;
isFitTypeSet = true;
isBackgroundSet = true;
} else {
clear();
Expand All @@ -146,10 +149,15 @@ void IqtFunctionTemplateModel::addFunction(std::string const &prefix, std::strin
newPrefix = *getExp1Prefix();
}
} else if (name == "StretchExp") {
if (hasStretchExponential())
if (hasFitType(FitType::StretchExponential))
throw std::runtime_error("Cannot add more stretched exponentials.");
m_fitType = FitType::StretchExponential;
newPrefix = *getStretchPrefix();
newPrefix = *getFitTypePrefix(m_fitType);
} else if (name == "TeixeiraWaterIqt") {
if (hasFitType(FitType::TeixeiraWaterIqt))
throw std::runtime_error("Cannot add another TeixeiraWaterIqt function.");
m_fitType = FitType::TeixeiraWaterIqt;
newPrefix = *getFitTypePrefix(m_fitType);
} else if (name == "FlatBackground") {
if (hasBackground())
throw std::runtime_error("Cannot add more backgrounds.");
Expand Down Expand Up @@ -213,7 +221,7 @@ void IqtFunctionTemplateModel::removeFunction(std::string const &prefix) {
m_exponentialType = ExponentialType::OneExponential;
return;
}
prefix1 = getStretchPrefix();
prefix1 = getFitTypePrefix(m_fitType);
if (prefix1 && *prefix1 == prefix) {
m_fitType = FitType::None;
return;
Expand All @@ -230,13 +238,14 @@ int IqtFunctionTemplateModel::numberOfExponentials() const { return static_cast<

bool IqtFunctionTemplateModel::hasExponential() const { return m_exponentialType != ExponentialType::None; }

bool IqtFunctionTemplateModel::hasStretchExponential() const { return m_fitType == FitType::StretchExponential; }
bool IqtFunctionTemplateModel::hasFitType() const { return m_fitType != FitType::None; }

bool IqtFunctionTemplateModel::hasFitType(FitType fitType) const { return m_fitType == fitType; }

void IqtFunctionTemplateModel::removeBackground() {
auto oldValues = getCurrentValues();
m_backgroundType = BackgroundType::None;
m_model->setFunctionString(buildFunctionString());
m_model->setGlobalParameters(makeGlobalList());
setModel();
setCurrentValues(oldValues);
}

Expand Down Expand Up @@ -278,8 +287,6 @@ void IqtFunctionTemplateModel::setResolution(const std::vector<std::pair<std::st
(void)fitResolutions;
}

void IqtFunctionTemplateModel::setQValues(const std::vector<double> &qValues) { (void)qValues; }

void IqtFunctionTemplateModel::tieIntensities() {
auto heightName = getParameterName(ParamID::STRETCH_HEIGHT);
if (!heightName)
Expand All @@ -299,7 +306,9 @@ std::optional<std::string> IqtFunctionTemplateModel::getPrefix(ParamID name) con
} else if (name <= ParamID::EXP2_LIFETIME) {
return getExp2Prefix();
} else if (name <= ParamID::STRETCH_STRETCHING) {
return getStretchPrefix();
return getFitTypePrefix(FitType::StretchExponential);
} else if (name <= ParamID::TWI_GAMMA) {
return getFitTypePrefix(FitType::TeixeiraWaterIqt);
} else {
return getBackgroundPrefix();
}
Expand All @@ -314,11 +323,16 @@ void IqtFunctionTemplateModel::applyParameterFunction(const std::function<void(P
paramFun(ParamID::EXP2_HEIGHT);
paramFun(ParamID::EXP2_LIFETIME);
}
if (hasStretchExponential()) {
if (hasFitType(FitType::StretchExponential)) {
paramFun(ParamID::STRETCH_HEIGHT);
paramFun(ParamID::STRETCH_LIFETIME);
paramFun(ParamID::STRETCH_STRETCHING);
}
if (hasFitType(FitType::TeixeiraWaterIqt)) {
paramFun(ParamID::TWI_AMPLITUDE);
paramFun(ParamID::TWI_TAU);
paramFun(ParamID::TWI_GAMMA);
}
if (hasBackground()) {
paramFun(ParamID::FLAT_BG_A0);
}
Expand All @@ -333,21 +347,31 @@ std::string IqtFunctionTemplateModel::buildStretchExpFunctionString() const {
"0,Lifetime>0,0<Stretching<1.001)";
}

std::string IqtFunctionTemplateModel::buildTeixeiraWaterIqtFunctionString(int const domainIndex) const {
auto const qValue = domainIndex < static_cast<int>(m_qValues.size()) ? m_qValues[domainIndex] : 0.4;
return "name=TeixeiraWaterIqt,Q=" + std::to_string(qValue) +
",Amp=1,Tau1=0.05,Gamma=1.2,constraints=(Amp>"
"0,Tau1>0,Gamma>0)";
}

std::string IqtFunctionTemplateModel::buildBackgroundFunctionString() const {
return "name=FlatBackground,A0=0,constraints=(A0>0)";
}

std::string IqtFunctionTemplateModel::buildFunctionString() const {
std::string IqtFunctionTemplateModel::buildFunctionString(int const domainIndex) const {
QStringList functions;
if (hasExponential()) {
functions << QString::fromStdString(buildExpDecayFunctionString());
}
if (numberOfExponentials() > 1) {
functions << QString::fromStdString(buildExpDecayFunctionString());
}
if (hasStretchExponential()) {
if (hasFitType(FitType::StretchExponential)) {
functions << QString::fromStdString(buildStretchExpFunctionString());
}
if (hasFitType(FitType::TeixeiraWaterIqt)) {
functions << QString::fromStdString(buildTeixeiraWaterIqtFunctionString(domainIndex));
}
if (hasBackground()) {
functions << QString::fromStdString(buildBackgroundFunctionString());
}
Expand All @@ -357,7 +381,7 @@ std::string IqtFunctionTemplateModel::buildFunctionString() const {
std::optional<std::string> IqtFunctionTemplateModel::getExp1Prefix() const {
if (!hasExponential())
return std::optional<std::string>();
if (numberOfExponentials() == 1 && !hasStretchExponential() && !hasBackground())
if (numberOfExponentials() == 1 && !hasFitType() && !hasBackground())
return std::string("");
return std::string("f0.");
}
Expand All @@ -368,8 +392,8 @@ std::optional<std::string> IqtFunctionTemplateModel::getExp2Prefix() const {
return std::string("f1.");
}

std::optional<std::string> IqtFunctionTemplateModel::getStretchPrefix() const {
if (!hasStretchExponential())
std::optional<std::string> IqtFunctionTemplateModel::getFitTypePrefix(FitType fitType) const {
if (!hasFitType(fitType))
return std::optional<std::string>();
if (!hasExponential() && !hasBackground())
return std::string("");
Expand All @@ -379,9 +403,9 @@ std::optional<std::string> IqtFunctionTemplateModel::getStretchPrefix() const {
std::optional<std::string> IqtFunctionTemplateModel::getBackgroundPrefix() const {
if (!hasBackground())
return std::optional<std::string>();
if (!hasExponential() && !hasStretchExponential())
if (!hasExponential() && !hasFitType())
return std::string("");
return "f" + std::to_string(numberOfExponentials() + (hasStretchExponential() ? 1 : 0)) + ".";
return "f" + std::to_string(numberOfExponentials() + (hasFitType() ? 1 : 0)) + ".";
}

} // namespace MantidQt::CustomInterfaces::Inelastic

0 comments on commit 08c54b8

Please sign in to comment.