diff --git a/Framework/CurveFitting/src/ComplexMatrix.cpp b/Framework/CurveFitting/src/ComplexMatrix.cpp
index 0f0759a1111c89c0da8805d82dd5cfc4ddebcea7..f3cd9ea3eb9aa2f1c9338e8bd4e07d680d9ad264 100644
--- a/Framework/CurveFitting/src/ComplexMatrix.cpp
+++ b/Framework/CurveFitting/src/ComplexMatrix.cpp
@@ -301,8 +301,6 @@ void ComplexMatrix::eigenSystemHermitian(GSLVector &eigenValues,
eigenValues.resize(n);
eigenVectors.resize(n, n);
auto workspace = gsl_eigen_hermv_alloc(n);
- // RT: I think there is a bug in this function. It returns garbage in
- // eigenvectors in a repeated call.
gsl_eigen_hermv(gsl(), eigenValues.gsl(), eigenVectors.gsl(), workspace);
gsl_eigen_hermv_free(workspace);
}
diff --git a/Framework/CurveFitting/src/GSLVector.cpp b/Framework/CurveFitting/src/GSLVector.cpp
index afe73a1b750506c2ead052320e6c8b4752db2ee1..62213159dbc6eefe3fa52366583b0f565e1f1e98 100644
--- a/Framework/CurveFitting/src/GSLVector.cpp
+++ b/Framework/CurveFitting/src/GSLVector.cpp
@@ -272,6 +272,7 @@ void GSLVector::sort(const std::vector<size_t> &indices) {
data[i] = m_data[indices[i]];
}
std::swap(m_data, data);
+ m_view = gsl_vector_view_array(m_data.data(), m_data.size());
}
/// Create a new GSLVector and move all data to it. Destroys this vector.
diff --git a/Framework/CurveFitting/src/IMWDomainCreator.cpp b/Framework/CurveFitting/src/IMWDomainCreator.cpp
index d2081c2ded0e494603103293e059eb8c472183e7..9a2c956805e59e0722cad8fd30986383613b16cd 100644
--- a/Framework/CurveFitting/src/IMWDomainCreator.cpp
+++ b/Framework/CurveFitting/src/IMWDomainCreator.cpp
@@ -437,58 +437,75 @@ void IMWDomainCreator::addFunctionValuesToWS(
function->function(*domain, *resultValues);
size_t nParams = function->nParams();
- // and errors
- SimpleJacobian J(nData, nParams);
- try {
- function->functionDeriv(*domain, J);
- } catch (...) {
- function->calNumericalDeriv(*domain, J);
- }
// the function should contain the parameter's covariance matrix
auto covar = function->getCovarianceMatrix();
-
- if (covar) {
- // if the function has a covariance matrix attached - use it for the errors
- const Kernel::Matrix<double> &C = *covar;
- // The formula is E = J * C * J^T
- // We don't do full 3-matrix multiplication because we only need the
- // diagonals of E
- std::vector<double> E(nData);
- for (size_t k = 0; k < nData; ++k) {
- double s = 0.0;
- for (size_t i = 0; i < nParams; ++i) {
- double tmp = J.get(k, i);
- s += C[i][i] * tmp * tmp;
- for (size_t j = i + 1; j < nParams; ++j) {
- s += J.get(k, i) * C[i][j] * J.get(k, j) * 2;
- }
+ bool hasErrors = false;
+ if (!covar) {
+ for (size_t j = 0; j < nParams; ++j) {
+ if (function->getError(j) != 0.0) {
+ hasErrors = true;
+ break;
}
- E[k] = s;
}
+ }
- double chi2 = function->getChiSquared();
- auto &yValues = ws->mutableY(wsIndex);
- auto &eValues = ws->mutableE(wsIndex);
- for (size_t i = 0; i < nData; i++) {
- yValues[i] = resultValues->getCalculated(i);
- eValues[i] = std::sqrt(E[i] * chi2);
+ if (covar || hasErrors) {
+ // and errors
+ SimpleJacobian J(nData, nParams);
+ try {
+ function->functionDeriv(*domain, J);
+ } catch (...) {
+ function->calNumericalDeriv(*domain, J);
}
+ if (covar) {
+ // if the function has a covariance matrix attached - use it for the errors
+ const Kernel::Matrix<double> &C = *covar;
+ // The formula is E = J * C * J^T
+ // We don't do full 3-matrix multiplication because we only need the
+ // diagonals of E
+ std::vector<double> E(nData);
+ for (size_t k = 0; k < nData; ++k) {
+ double s = 0.0;
+ for (size_t i = 0; i < nParams; ++i) {
+ double tmp = J.get(k, i);
+ s += C[i][i] * tmp * tmp;
+ for (size_t j = i + 1; j < nParams; ++j) {
+ s += J.get(k, i) * C[i][j] * J.get(k, j) * 2;
+ }
+ }
+ E[k] = s;
+ }
- } else {
- // otherwise use the parameter errors which is OK for uncorrelated
- // parameters
- auto &yValues = ws->mutableY(wsIndex);
- auto &eValues = ws->mutableE(wsIndex);
- for (size_t i = 0; i < nData; i++) {
- yValues[i] = resultValues->getCalculated(i);
- double err = 0.0;
- for (size_t j = 0; j < nParams; ++j) {
- double d = J.get(i, j) * function->getError(j);
- err += d * d;
+ double chi2 = function->getChiSquared();
+ auto &yValues = ws->mutableY(wsIndex);
+ auto &eValues = ws->mutableE(wsIndex);
+ for (size_t i = 0; i < nData; i++) {
+ yValues[i] = resultValues->getCalculated(i);
+ eValues[i] = std::sqrt(E[i] * chi2);
+ }
+
+ } else {
+ // otherwise use the parameter errors which is OK for uncorrelated
+ // parameters
+ auto &yValues = ws->mutableY(wsIndex);
+ auto &eValues = ws->mutableE(wsIndex);
+ for (size_t i = 0; i < nData; i++) {
+ yValues[i] = resultValues->getCalculated(i);
+ double err = 0.0;
+ for (size_t j = 0; j < nParams; ++j) {
+ double d = J.get(i, j) * function->getError(j);
+ err += d * d;
+ }
+ eValues[i] = std::sqrt(err);
}
- eValues[i] = std::sqrt(err);
}
+ } else {
+ // No errors
+ auto &yValues = ws->mutableY(wsIndex);
+ for (size_t i = 0; i < nData; i++) {
+ yValues[i] = resultValues->getCalculated(i);
+ }
}
}
diff --git a/Framework/CurveFitting/test/Functions/CrystalFieldHeatCapacityTest.h b/Framework/CurveFitting/test/Functions/CrystalFieldHeatCapacityTest.h
index d78f5edce30e2f601d2f3f984266361c3be2f7cd..379d12b6c59e5c081642fe02370b7b33ca3d6c56 100644
--- a/Framework/CurveFitting/test/Functions/CrystalFieldHeatCapacityTest.h
+++ b/Framework/CurveFitting/test/Functions/CrystalFieldHeatCapacityTest.h
@@ -3,11 +3,14 @@
#include <cxxtest/TestSuite.h>
+#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/ParameterTie.h"
+#include "MantidCurveFitting/Algorithms/EvaluateFunction.h"
#include "MantidCurveFitting/Functions/CrystalFieldHeatCapacity.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
using namespace Mantid;
using namespace Mantid::API;
@@ -75,6 +78,32 @@ public:
}
TS_ASSERT_EQUALS(nTies, 1); // Fixed values not ties.
}
+
+ void test_evaluate_function() {
+ auto ws = WorkspaceCreationHelper::create2DWorkspaceBinned(1, 100, 1.0, 3.0);
+ std::string fun = "name=CrystalFieldHeatCapacity,Ion=Ce,Symmetry=C2v,"
+ "B20=0.37,B22=3.9, B40=-0.03,B42=-0.1,B44=-0.12, "
+ "ties=(BmolX=0,BmolY=0,BmolZ=0,BextX=0,BextY=0,BextZ=BextX)";
+
+ Algorithms::EvaluateFunction eval;
+ eval.initialize();
+ eval.setProperty("Function", fun);
+ eval.setProperty("InputWorkspace", ws);
+ eval.setPropertyValue("OutputWorkspace", "out");
+ eval.execute();
+
+ auto out =
+ API::AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>("out");
+
+ auto &y = out->histogram(1).counts();
+ TS_ASSERT_DELTA(y[10], 0.0305, 1e-4);
+ TS_ASSERT_DELTA(y[30], 3.7753, 1e-4);
+ TS_ASSERT_DELTA(y[70], 5.1547, 1e-4);
+ TS_ASSERT_DELTA(y[99], 3.4470, 1e-4);
+
+ API::AnalysisDataService::Instance().clear();
+ }
+
};
#endif /*CRYSTALFIELDHEATCAPACITYTEST_H_*/
diff --git a/scripts/test/CrystalFieldTest.py b/scripts/test/CrystalFieldTest.py
index 4335307448c4bebc32059000d0dfefc11d706a0a..2662c3366b32908827cfa767c1c80cec3876a6fc 100644
--- a/scripts/test/CrystalFieldTest.py
+++ b/scripts/test/CrystalFieldTest.py
@@ -1331,8 +1331,8 @@ class CrystalFieldFitTest(unittest.TestCase):
Temperature=44.0, FWHM=1.0, ResolutionModel=rm, FWHMVariation=0.3)
sp = cf.getSpectrum()
self.assertAlmostEqual(cf.peaks.param[0]['FWHM'], 1.0, 8)
- self.assertAlmostEqual(cf.peaks.param[1]['FWHM'], 1.58101468, 8)
- self.assertAlmostEqual(cf.peaks.param[2]['FWHM'], 1.884945866, 8)
+ self.assertAlmostEqual(cf.peaks.param[1]['FWHM'], 1.58958, 1)
+ self.assertAlmostEqual(cf.peaks.param[2]['FWHM'], 1.85644, 1)
def test_ResolutionModel_set_multi(self):
from CrystalField import ResolutionModel, CrystalField, CrystalFieldFit