From d3bec142f13552d4362d1fd5266f3e8d4021174b Mon Sep 17 00:00:00 2001
From: Roman Tolchenov <roman.tolchenov@stfc.ac.uk>
Date: Mon, 21 Feb 2011 15:28:24 +0000
Subject: [PATCH] Added algorithm FFTDerivative. re #2517
---
.../Framework/Algorithms/CMakeLists.txt | 7 +-
.../inc/MantidAlgorithms/FFTDerivative.h | 61 ++++++++
.../Algorithms/src/FFTDerivative.cpp | 132 ++++++++++++++++++
.../Algorithms/test/FFTDerivativeTest.h | 89 ++++++++++++
4 files changed, 287 insertions(+), 2 deletions(-)
create mode 100644 Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FFTDerivative.h
create mode 100644 Code/Mantid/Framework/Algorithms/src/FFTDerivative.cpp
create mode 100644 Code/Mantid/Framework/Algorithms/test/FFTDerivativeTest.h
diff --git a/Code/Mantid/Framework/Algorithms/CMakeLists.txt b/Code/Mantid/Framework/Algorithms/CMakeLists.txt
index ba00fe0db53..5344bd916a6 100644
--- a/Code/Mantid/Framework/Algorithms/CMakeLists.txt
+++ b/Code/Mantid/Framework/Algorithms/CMakeLists.txt
@@ -49,6 +49,7 @@ set ( SRC_FILES src/AbsorptionCorrection.cpp
src/ExtractMasking.cpp
src/ExtractSingleSpectrum.cpp
src/FFT.cpp
+ src/FFTDerivative.cpp
src/FFTSmooth2.cpp
src/FFTSmooth.cpp
src/FilterBadPulses.cpp
@@ -184,6 +185,7 @@ set ( INC_FILES inc/MantidAlgorithms/AbsorptionCorrection.h
inc/MantidAlgorithms/ExtractMasking.h
inc/MantidAlgorithms/ExtractSingleSpectrum.h
inc/MantidAlgorithms/FFT.h
+ inc/MantidAlgorithms/FFTDerivative.h
inc/MantidAlgorithms/FFTSmooth2.h
inc/MantidAlgorithms/FFTSmooth.h
inc/MantidAlgorithms/FilterBadPulses.h
@@ -317,6 +319,7 @@ set ( TEST_FILES test/AddSampleLogTest.h
test/ExtractSingleSpectrumTest.h
test/FFTSmooth2Test.h
test/FFTTest.h
+ test/FFTDerivativeTest.h
test/FilterByLogValueTest.h
test/FilterByTimeTest.h
test/FindCenterOfMassPositionTest.h
@@ -335,7 +338,7 @@ set ( TEST_FILES test/AddSampleLogTest.h
test/IdentifyNoisyDetectorsTest.h
test/IntegrationTest.h
test/InterpolatingRebinTest.h
- test/IQTransformTest.h
+ test/IQTransformTest.h
test/LogarithmTest.h
test/MaskBinsTest.h
test/MedianDetectorTestTest.h
@@ -389,7 +392,7 @@ set ( TEST_FILES test/AddSampleLogTest.h
test/UnaryOperationTest.h
test/UnGroupWorkspacesTest.h
test/UnwrapTest.h
- test/WeightedMeanTest.h
+ test/WeightedMeanTest.h
test/WorkspaceGroupTest.h
)
diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FFTDerivative.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FFTDerivative.h
new file mode 100644
index 00000000000..fb9c5c23952
--- /dev/null
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FFTDerivative.h
@@ -0,0 +1,61 @@
+#ifndef FFTDERIVATIVE_H_
+#define FFTDERIVATIVE_H_
+
+#include "MantidAPI/Algorithm.h"
+
+namespace Mantid
+{
+namespace Algorithms
+{
+
+/** Calculates derivatives of the spectra in a MatrixWorkspace using a Fast Fourier Transform.
+
+ @author Roman Tolchenov
+ @date 21/02/2011
+
+ Copyright © 2008 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
+
+ This file is part of Mantid.
+
+ Mantid is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ Mantid is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ File change history is stored at: <https://svn.mantidproject.org/mantid/trunk/Code/Mantid>
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+ */
+class FFTDerivative : public Mantid::API::Algorithm
+{
+public:
+ /// (Empty) Constructor
+ FFTDerivative() : Mantid::API::Algorithm() {}
+ /// Virtual destructor
+ virtual ~FFTDerivative() {}
+ /// Algorithm's name
+ virtual const std::string name() const { return "FFTDerivative"; }
+ /// Algorithm's version
+ virtual int version() const { return (1); }
+ /// Algorithm's category for identification
+ virtual const std::string category() const { return "General"; }
+
+private:
+ /// Initialisation code
+ void init();
+ ///Execution code
+ void exec();
+
+};
+
+} // FFT
+} // Mandid
+
+#endif /*DERIVATIVE_H_*/
diff --git a/Code/Mantid/Framework/Algorithms/src/FFTDerivative.cpp b/Code/Mantid/Framework/Algorithms/src/FFTDerivative.cpp
new file mode 100644
index 00000000000..61d2070ebb0
--- /dev/null
+++ b/Code/Mantid/Framework/Algorithms/src/FFTDerivative.cpp
@@ -0,0 +1,132 @@
+#include "MantidAlgorithms/FFTDerivative.h"
+
+#include <algorithm>
+#include <functional>
+
+namespace Mantid
+{
+namespace Algorithms
+{
+
+// Register the algorithm into the AlgorithmFactory
+DECLARE_ALGORITHM(FFTDerivative)
+
+using namespace Mantid::Kernel;
+using namespace Mantid::API;
+
+void FFTDerivative::init()
+{
+ declareProperty(new WorkspaceProperty<>("InputWorkspace","",Direction::Input,"Input workspace for differentiation"));
+ declareProperty(new WorkspaceProperty<>("OutputWorkspace","",Direction::Output,"Workspace with result derivatives"));
+}
+
+void FFTDerivative::exec()
+{
+
+ MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
+ MatrixWorkspace_sptr outWS;
+
+ int n = inWS->getNumberHistograms();
+ API::Progress progress(this,0,1,n);
+
+ int ny = inWS->readY(0).size();
+ int nx = inWS->readX(0).size();
+
+ // Workspace for holding a copy of a spectrum. Each spectrum is symmetrized to minimize
+ // possible edge effects.
+ MatrixWorkspace_sptr copyWS = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+ (Mantid::API::WorkspaceFactory::Instance().create(inWS,1,nx+ny,ny+ny));
+
+ bool isHist = (nx != ny);
+
+ for(int spec = 0; spec < n; ++spec)
+ {
+
+ const Mantid::MantidVec& x0 = inWS->readX(spec);
+ const Mantid::MantidVec& y0 = inWS->readY(spec);
+
+ Mantid::MantidVec& x1 = copyWS->dataX(0);
+ Mantid::MantidVec& y1 = copyWS->dataY(0);
+
+ double xx = 2*x0[0];
+
+ x1[ny] = x0[0];
+ y1[ny] = y0[0];
+
+ for(int i = 1; i < ny; ++i)
+ {
+ int j1 = ny - i;
+ int j2 = ny + i;
+ x1[j1] = xx - x0[i];
+ x1[j2] = x0[i];
+ y1[j1] = y1[j2] = y0[i];
+ }
+
+ x1[0] = 2*x1[1] - x1[2];
+ y1[0] = y0.back();
+
+ if (isHist)
+ {
+ x1[y1.size()] = x0[ny];
+ }
+
+ // Transform symmetrized spectrum
+ IAlgorithm_sptr fft = createSubAlgorithm("RealFFT");
+ fft->setProperty("InputWorkspace",copyWS);
+ fft->setProperty("WorkspaceIndex",0);
+ fft->setProperty("Transform","Forward");
+ fft->execute();
+
+ MatrixWorkspace_sptr transWS = fft->getProperty("OutputWorkspace");
+
+ Mantid::MantidVec& nu = transWS->dataX(0);
+ Mantid::MantidVec& re = transWS->dataY(0);
+ Mantid::MantidVec& im = transWS->dataY(1);
+
+ // Multiply the transform by 2*pi*i*w
+ for(int j=0; j < re.size(); ++j)
+ {
+ double w = 2 * M_PI * nu[j];
+ double a = re[j]*w;
+ double b = -im[j]*w;
+ re[j] = b;
+ im[j] = a;
+ }
+
+ // Inverse transform
+ fft = createSubAlgorithm("RealFFT");
+ fft->setProperty("InputWorkspace",transWS);
+ fft->setProperty("Transform","Backward");
+ fft->execute();
+
+ transWS = fft->getProperty("OutputWorkspace");
+
+ int m2 = transWS->readY(0).size() / 2;
+ int my = m2 + (transWS->readY(0).size() % 2 ? 1 : 0);
+ int mx = my + (transWS->isHistogramData() ? 1 : 0);
+
+ if (!outWS)
+ {
+ outWS = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+ (Mantid::API::WorkspaceFactory::Instance().create(inWS,n,mx,my));
+ }
+
+ // Save the upper half of the inverse transform for output
+ Mantid::MantidVec& x = outWS->dataX(spec);
+ Mantid::MantidVec& y = outWS->dataY(spec);
+ double dx = x1[0];
+ std::copy(transWS->dataX(0).begin() + m2,transWS->dataX(0).end(),x.begin());
+ std::transform(x.begin(),x.end(),x.begin(),std::bind2nd(std::plus<double>(),dx));
+ std::copy(transWS->dataY(0).begin() + m2,transWS->dataY(0).end(),y.begin());
+
+ progress.report();
+
+ }
+
+ setProperty("OutputWorkspace",outWS);
+
+}
+
+} // Algorithms
+} // Mandid
+
diff --git a/Code/Mantid/Framework/Algorithms/test/FFTDerivativeTest.h b/Code/Mantid/Framework/Algorithms/test/FFTDerivativeTest.h
new file mode 100644
index 00000000000..540f044faea
--- /dev/null
+++ b/Code/Mantid/Framework/Algorithms/test/FFTDerivativeTest.h
@@ -0,0 +1,89 @@
+#ifndef FFT_DERIVATIVE_TEST_H_
+#define FFT_DERIVATIVE_TEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidAlgorithms/FFTDerivative.h"
+#include "MantidAPI/AnalysisDataService.h"
+#include "MantidAPI/FrameworkManager.h"
+#include "MantidAPI/MatrixWorkspace.h"
+#include "MantidDataObjects/Workspace2D.h"
+
+using namespace Mantid;
+using namespace Mantid::API;
+
+class FFTDerivativeTest : public CxxTest::TestSuite
+{
+public:
+ FFTDerivativeTest()
+ {
+ }
+ ~FFTDerivativeTest()
+ {
+ }
+
+ void testGaussianDerivative()
+ {
+ const int N = 100;
+
+ createWS(N,0,"exp");
+
+ IAlgorithm* fft = Mantid::API::FrameworkManager::Instance().createAlgorithm("FFTDerivative");
+ fft->initialize();
+ fft->setPropertyValue("InputWorkspace","FFTDerivative_WS_exp");
+ fft->setPropertyValue("OutputWorkspace","FFTDerivative_out");
+ fft->execute();
+
+ MatrixWorkspace_sptr fWS = boost::dynamic_pointer_cast<MatrixWorkspace>
+ (AnalysisDataService::Instance().retrieve("FFTDerivative_out"));
+
+ const MantidVec& X = fWS->readX(0);
+ const MantidVec& Y = fWS->readY(0);
+
+ TS_ASSERT_EQUALS(Y.size(),101);
+
+ for(int i=0; i < Y.size(); ++i)
+ {
+ double xx = X[i] - 5.1;
+ TS_ASSERT_DELTA( Y[i],(-2 * xx * exp(-(xx*xx)*2)*2),0.000001);
+ }
+
+ FrameworkManager::Instance().deleteWorkspace("FFTDerivative_WS_exp");
+ FrameworkManager::Instance().deleteWorkspace("FFTDerivative_out");
+ }
+
+
+private:
+
+ MatrixWorkspace_sptr createWS(int n,int dn,const std::string& name)
+ {
+ FrameworkManager::Instance();
+ Mantid::DataObjects::Workspace2D_sptr ws = boost::dynamic_pointer_cast<Mantid::DataObjects::Workspace2D>
+ (WorkspaceFactory::Instance().create("Workspace2D",1,n+dn,n));
+
+ const double dX = 10.0 / (n-1);
+ const double x0 = 0.;
+ const double c = 5.1;
+ Mantid::MantidVec& X = ws->dataX(0);
+ Mantid::MantidVec& Y = ws->dataY(0);
+ Mantid::MantidVec& E = ws->dataE(0);
+
+ for(int i=0;i<n;i++)
+ {
+ double x = x0 + dX*(i);
+ X[i] = x;
+ Y[i] = exp(-(x-c)*(x-c)*2);
+ E[i] = 1.;
+ }
+
+ if (dn > 0)
+ X[n] = X[n-1] + dX;
+
+ AnalysisDataService::Instance().add("FFTDerivative_WS_"+name,ws);
+
+ return ws;
+ }
+
+};
+
+#endif /*FFT_DERIVATIVE_TEST_H_*/
--
GitLab