Skip to content

Commit

Permalink
Re #11131. ChebfunBase unit tests. Formatting.
Browse files Browse the repository at this point in the history
  • Loading branch information
@mantid-roman
mantid-roman committed Mar 4, 2015
1 parent cf8dd2a commit a664f02
Show file tree
Hide file tree
Showing 11 changed files with 1,049 additions and 998 deletions.
2 changes: 2 additions & 0 deletions Code/Mantid/Framework/CurveFitting/CMakeLists.txt
View file
Expand Up @@ -162,6 +162,7 @@ set ( INC_FILES
inc/MantidCurveFitting/Gaussian.h
inc/MantidCurveFitting/GaussianComptonProfile.h
inc/MantidCurveFitting/GramCharlierComptonProfile.h
inc/MantidCurveFitting/HalfComplex.h
inc/MantidCurveFitting/IkedaCarpenterPV.h
inc/MantidCurveFitting/Jacobian.h
inc/MantidCurveFitting/LeBailFit.h
Expand Down Expand Up @@ -230,6 +231,7 @@ set ( TEST_FILES
BoundaryConstraintTest.h
CalculateGammaBackgroundTest.h
CalculateMSVesuvioTest.h
ChebfunBaseTest.h
ChebyshevTest.h
CompositeFunctionTest.h
ComptonPeakProfileTest.h
Expand Down
281 changes: 162 additions & 119 deletions Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ChebfunBase.h
View file
Expand Up @@ -10,141 +10,184 @@

namespace Mantid {

namespace API{
class IFunction;
namespace API {
class IFunction;
}

namespace CurveFitting {

typedef std::vector<double> ChebfunVec;
/// Type of the approximated function
typedef std::function<double(double)> ChebfunFunctionType;


/**
* @brief The ChebfunBase class provides a base for a single chebfun.
*
* It keeps the x-points (knots? or just points?) and various weights.
* A single chebfun base can be shared by multiple chebfuns.
*/
class MANTID_CURVEFITTING_DLL ChebfunBase
{
The ChebfunBase class provides a base for function approximation
with Chebyshev polynomials.
A smooth function on a finite interval [a,b] can be approximated
by a Chebyshev expansion of order n. Finding an approximation is
very easy: the function needs to be evaluated at n+1 specific x-
points. These n+1 values can be used to interpolate the function
at any x-point in interval [a,b]. This is done by calling the fit(...)
method.
Different functions require different polynomial orders to reach
the same accuracy of approximation. Static method bestFit(...) tries
to find the smallest value of n that provides the required accuracy.
If it fails to find an n smaller than some maximum number it returns
an empty shared pointer.
Knowing the vector of the function values (P) at the n+1 base x-points and the
related vector of the Chebyshev expansion coefficients (A) (claculated
by calcA(...) method) allows one to perform various manipulations on
the approximation:
- algebraic operations: +,-,*,/
- applying a function
- root finding
- differentiation
- integration
- convolution
- solving of (integro-)differential equations
- etc
This calss doesn't represent a function approximation itself but keeps
proerties that can be shared by multiple approximations.
This class is based on the ideas from the Chebfun matlab package
(http://www.chebfun.org/).
Copyright &copy; 2007-8 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
National Laboratory & European Spallation Source
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://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class MANTID_CURVEFITTING_DLL ChebfunBase {
public:
ChebfunBase(size_t n, double start, double end, double tolerance = 0.0);
/// Copy constructor
ChebfunBase( const ChebfunBase& other );
/// Get the polynomial order of the chebfun based on this base.
size_t order() const {return m_n;}
/// Get the size of the base which is the number of x-points.
size_t size() const {return m_x.size();}
/// Start of the interval
double startX() const {return m_x.front();}
/// End of the interval
double endX() const {return m_x.back();}
/// Get the width of the interval
double width() const {return endX() - startX();}
/// Get a reference to the x-points
const std::vector<double>& xPoints() const {return m_x;}
/// Get a reference to the integration weights
const std::vector<double>& integrationWeights() const;
/// Calculate an integral
double integrate(const ChebfunVec& p) const;
/// Calculate expansion coefficients
ChebfunVec calcA(const ChebfunVec& p)const;
/// Calculate function values
ChebfunVec calcP(const ChebfunVec& a)const;
/// Calculate function values at chebfun x-points
ChebfunVec fit(ChebfunFunctionType f ) const;
/// Calculate function values at chebfun x-points
ChebfunVec fit(const API::IFunction& f ) const;
/// Test an array of chebfun coefficients for convergence
bool isConverged(const std::vector<double>& a, double maxA = 0.0);

/// Evaluate a function
double eval(double x, const std::vector<double> &p) const;
/// Evaluate a function
void evalVector(const std::vector<double> &x, const std::vector<double> &p, std::vector<double> &res) const;
/// Evaluate a function
std::vector<double> evalVector(const std::vector<double> &x, const std::vector<double> &p) const;
/// Evaluate a function for a range of x-values.
template<class XIter, class ResIter>
void evalIter(XIter xbegin, XIter xend, const std::vector<double> &p, ResIter res) const;
/// Calculate the derivative
void derivative(const std::vector<double>& a, std::vector<double>& aout) const;
/// Calculate the integral
boost::shared_ptr<ChebfunBase> integral(const std::vector<double>& a, std::vector<double>& aout) const;
std::vector<double> ChebfunBase::roots(const std::vector<double>& p) const;

/// Fit a function until full convergence
static boost::shared_ptr<ChebfunBase> bestFit(double start, double end, ChebfunFunctionType, ChebfunVec& p, ChebfunVec& a, double maxA = 0.0, double tolerance = 0.0, size_t maxSize = 0 );
/// Fit a function until full convergence
static boost::shared_ptr<ChebfunBase> bestFit(double start, double end,const API::IFunction&, ChebfunVec& p, ChebfunVec& a, double maxA = 0.0, double tolerance = 0.0, size_t maxSize = 0 );
/// Tolerance for comparing doubles
double tolerance() {return m_tolerance;}

std::vector<double> linspace(size_t n) const;
/// Get an interpolating matrix
GSLMatrix createInterpolatingMatrix(const std::vector<double> &x, bool isZeroOutside = false) const;
GSLMatrix createConvolutionMatrix(ChebfunFunctionType fun) const;
std::vector<double> smooth(const std::vector<double> &xvalues, const std::vector<double> &yvalues) const;
ChebfunBase(size_t n, double start, double end, double tolerance = 0.0);
/// Copy constructor
ChebfunBase(const ChebfunBase &other);
/// Get the polynomial order of this base.
size_t order() const { return m_n; }
/// Get the size of the base which is the number of x-points.
size_t size() const { return m_x.size(); }
/// Start of the interval
double startX() const { return m_x.front(); }
/// End of the interval
double endX() const { return m_x.back(); }
/// Get the width of the interval
double width() const { return endX() - startX(); }
/// Get a reference to the x-points
const std::vector<double> &xPoints() const { return m_x; }
/// Get a reference to the integration weights
const std::vector<double> &integrationWeights() const;
/// Calculate an integral
double integrate(const std::vector<double> &p) const;
/// Calculate expansion coefficients
std::vector<double> calcA(const std::vector<double> &p) const;
/// Calculate function values
std::vector<double> calcP(const std::vector<double> &a) const;
/// Calculate function values at chebfun x-points
std::vector<double> fit(ChebfunFunctionType f) const;
/// Calculate function values at chebfun x-points
std::vector<double> fit(const API::IFunction &f) const;

/// Evaluate a function
double eval(double x, const std::vector<double> &p) const;
/// Evaluate a function
void evalVector(const std::vector<double> &x, const std::vector<double> &p,
std::vector<double> &res) const;
/// Evaluate a function
std::vector<double> evalVector(const std::vector<double> &x,
const std::vector<double> &p) const;
/// Calculate the derivative
void derivative(const std::vector<double> &a,
std::vector<double> &aout) const;
/// Calculate the integral
boost::shared_ptr<ChebfunBase> integral(const std::vector<double> &a,
std::vector<double> &aout) const;
/// Find all roots of a function on this interval
std::vector<double> ChebfunBase::roots(const std::vector<double> &a) const;

/// Fit a function until full convergence
static boost::shared_ptr<ChebfunBase>
bestFit(double start, double end, ChebfunFunctionType, std::vector<double> &p,
std::vector<double> &a, double maxA = 0.0, double tolerance = 0.0,
size_t maxSize = 0);
/// Fit a function until full convergence
static boost::shared_ptr<ChebfunBase>
bestFit(double start, double end, const API::IFunction &,
std::vector<double> &p, std::vector<double> &a, double maxA = 0.0,
double tolerance = 0.0, size_t maxSize = 0);
/// Tolerance for comparing doubles
double tolerance() { return m_tolerance; }

/// Create a vector of x values linearly spaced on the approximation interval
std::vector<double> linspace(size_t n) const;

private:
/// Private assingment operator to stress the immutability of ChebfunBase.
ChebfunBase& operator=( const ChebfunBase& other );
/// Calculate the x-values based on the (start,end) interval.
void calcX();
/// Calculate the integration weights
void calcIntegrationWeights() const;

/// Calculate function values at odd-valued indices of chebfun x-points
ChebfunVec fitOdd(ChebfunFunctionType f, ChebfunVec& p) const;
/// Calculate function values at odd-valued indices of chebfun x-points
ChebfunVec fitOdd(const API::IFunction& f, ChebfunVec& p) const;
/// Test an array of chebfun coefficients for convergence
static bool hasConverged(const std::vector<double>& a, double maxA, double tolerance);
template<class FunctionType>
static boost::shared_ptr<ChebfunBase> bestFitTempl(double start, double end, FunctionType f, ChebfunVec& p, ChebfunVec& a , double maxA, double tolerance, size_t maxSize);

/// Actual tolerance in comparing doubles
const double m_tolerance;
/// Number of points on the x-axis.
size_t m_n;
/// Start of the interval
double m_start;
/// End of the interval
double m_end;
/// The x-points
std::vector<double> m_x;
/// The barycentric weights.
std::vector<double> m_bw;
/// Integration weights
mutable std::vector<double> m_integrationWeights;
/// Maximum tolerance in comparing doubles
static const double g_tolerance;
/// Maximum number of (x) points in a base.
static const size_t g_maxNumberPoints;

/// Private assingment operator to stress the immutability of ChebfunBase.
ChebfunBase &operator=(const ChebfunBase &other);
/// Calculate the x-values based on the (start,end) interval.
void calcX();
/// Calculate the integration weights
void calcIntegrationWeights() const;

/// Calculate function values at odd-valued indices of the base x-points
std::vector<double> fitOdd(ChebfunFunctionType f,
std::vector<double> &p) const;
/// Calculate function values at odd-valued indices of the base x-points
std::vector<double> fitOdd(const API::IFunction &f,
std::vector<double> &p) const;
/// Test an array of Chebyshev coefficients for convergence
static bool hasConverged(const std::vector<double> &a, double maxA,
double tolerance, size_t shift = 0);
/// Templated implementation of bestFit method
template <class FunctionType>
static boost::shared_ptr<ChebfunBase>
bestFitTempl(double start, double end, FunctionType f, std::vector<double> &p,
std::vector<double> &a, double maxA, double tolerance,
size_t maxSize);

/// Actual tolerance in comparing doubles
const double m_tolerance;
/// Polynomial order.
size_t m_n;
/// Start of the interval
double m_start;
/// End of the interval
double m_end;
/// The x-points
std::vector<double> m_x;
/// The barycentric weights.
std::vector<double> m_bw;
/// The integration weights
mutable std::vector<double> m_integrationWeights;
/// Maximum tolerance in comparing doubles
static const double g_tolerance;
/// Maximum number of (x) points in a base.
static const size_t g_maxNumberPoints;
};

/**
* Evaluate a function for a range of x-values.
* @param xbegin :: Iterator of the start of a range of x-values
* @param xend :: Iterator of the end of a range of x-values
* @param p :: The function parameters.
* @param res :: Iterator to the start of the results container. The size of the container must
* not be smaller than distance(xbegin,xend).
*/
template<class XIter, class ResIter>
void ChebfunBase::evalIter(XIter xbegin, XIter xend, const std::vector<double> &p, ResIter res) const
{
using namespace std::placeholders;
std::transform( xbegin, xend, res, std::bind(&ChebfunBase::eval, this, _1, p) );
}

typedef boost::shared_ptr<ChebfunBase> ChebfunBase_sptr;

} // CurveFitting
} // Mantid


#endif // MANTID_CURVEFITTING_CHEBFUNBASE_H
70 changes: 0 additions & 70 deletions Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/Convolution.h
View file
Expand Up @@ -41,76 +41,6 @@ Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class DLLExport Convolution : public API::CompositeFunction {
public:
/**
* Class for helping to read the transformed data. It represent an output of
* the
* GSL real fast fourier transform routine. The routine transforms an array of
* n
* real numbers into an array of about n/2 complex numbers which are the
* amplitudes of
* the positive frequencies of the full complex fourier transform.
*/
class HalfComplex {
size_t m_size; ///< size of the transformed data
double *m_data; ///< pointer to the transformed data
bool m_even; ///< true if the size of the original data is even
public:
/**
* Constructor.
* @param data :: A pointer to the transformed complex data
* @param n :: The size of untransformed real data
*/
HalfComplex(double *data, const size_t &n)
: m_size(n / 2 + 1), m_data(data), m_even(n / 2 * 2 == n) {}
/// Returns the size of the transform
size_t size() const { return m_size; }
/**
* The real part of i-th transform coefficient
* @param i :: The index of the complex transform coefficient
* @return The real part
*/
double real(size_t i) const {
if (i >= m_size)
return 0.;
if (i == 0)
return m_data[0];
return m_data[2 * i - 1];
}
/**
* The imaginary part of i-th transform coefficient
* @param i :: The index of the complex transform coefficient
* @return The imaginary part
*/
double imag(size_t i) const {
if (i >= m_size)
return 0.;
if (i == 0)
return 0;
if (m_even && i == m_size - 1)
return 0;
return m_data[2 * i];
}
/**
* Set a new value for i-th complex coefficient
* @param i :: The index of the coefficient
* @param re :: The real part of the new value
* @param im :: The imaginary part of the new value
*/
void set(size_t i, const double &re, const double &im) {
if (i >= m_size)
return;
if (i == 0) // this is purely real
{
m_data[0] = re;
} else if (m_even && i == m_size - 1) // this is also purely real
{
m_data[2 * i - 1] = re;
} else {
m_data[2 * i - 1] = re;
m_data[2 * i] = im;
}
}
};

/// Constructor
Convolution();
Expand Down

0 comments on commit a664f02

Please sign in to comment.