fix some coding style isses in new class XYTabulated2DFunction and its test
This commit is contained in:
Andreas Lauser
parent
20b6068a87
commit
7a4c6546b4
@ -42,91 +42,69 @@
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \brief Implements a scalar function that depends on two variables, which are sampled
|
||||
* in X and Y directions with sampling points, respectively.
|
||||
* The table might be able to be extrapolated in certian directions.
|
||||
* \brief Implements a function that depends on two variables.
|
||||
*
|
||||
* The function is sampled in regular intervals in both directions, i.e., the
|
||||
* interpolation cells are rectangles. The table can be specified to be extrapolated in
|
||||
* either direction.
|
||||
*/
|
||||
|
||||
template <class Scalar>
|
||||
class XYTabulated2DFunction
|
||||
{
|
||||
|
||||
public:
|
||||
XYTabulated2DFunction()
|
||||
{ }
|
||||
|
||||
XYTabulated2DFunction(const std::vector<Scalar>& x_pos,
|
||||
const std::vector<Scalar>& y_pos,
|
||||
const std::vector<std::vector<Scalar> >& data,
|
||||
const bool x_extrapolate = false,
|
||||
const bool y_extrapolate = false)
|
||||
: xPos_(x_pos)
|
||||
, yPos_(y_pos)
|
||||
template <class DataContainer>
|
||||
XYTabulated2DFunction(const std::vector<Scalar>& xPos,
|
||||
const std::vector<Scalar>& yPos,
|
||||
const DataContainer& data,
|
||||
const bool xExtrapolate = false,
|
||||
const bool yExtrapolate = false)
|
||||
: xPos_(xPos)
|
||||
, yPos_(yPos)
|
||||
, samples_(data)
|
||||
, xExtrapolate_(x_extrapolate)
|
||||
, yExtrapolate_(y_extrapolate)
|
||||
, xExtrapolate_(xExtrapolate)
|
||||
, yExtrapolate_(yExtrapolate)
|
||||
{
|
||||
// make sure the size is correct
|
||||
if (numX() != samples_.size()) {
|
||||
throw std::runtime_error("numX() is not equal to the number of rows of the sample points");
|
||||
#ifndef NDEBUG
|
||||
// in debug mode, ensure that the x and y positions arrays are strictly
|
||||
// mononically increasing.
|
||||
for (unsigned i = 0; i < xPos.size() - 1; ++ i) {
|
||||
if (xPos[i + 1] <= xPos[i])
|
||||
throw std::runtime_error("The array for the x-positions is not strictly increasing!");
|
||||
}
|
||||
|
||||
for (unsigned ix = 0; ix < numX(); ++ix) {
|
||||
if (samples_[ix].size() != numY()) {
|
||||
for (unsigned i = 0; i < yPos.size() - 1; ++ i) {
|
||||
if (yPos[i + 1] <= yPos[i])
|
||||
throw std::runtime_error("The array for the y-positions is not strictly increasing!");
|
||||
}
|
||||
#endif
|
||||
|
||||
// make sure the size is correct
|
||||
if (numX() != samples_.size())
|
||||
throw std::runtime_error("numX() is not equal to the number of rows of the sampling points");
|
||||
|
||||
for (unsigned xIdx = 0; xIdx < numX(); ++xIdx) {
|
||||
if (samples_[xIdx].size() != numY()) {
|
||||
std::ostringstream oss;
|
||||
oss << "the " << ix << "th row of the sample points has different number of data from numY() ";
|
||||
oss << "The " << xIdx << "-th row of the sampling points has different size than numY() ";
|
||||
throw std::runtime_error(oss.str());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Evaluate the function at a given (x,y) position.
|
||||
*
|
||||
* If this method is called for a value outside of the tabulated
|
||||
* range, and extrapolation is not allowed in the corresponding direction,
|
||||
* a \c Opm::NumericalIssue exception is thrown.
|
||||
* \brief Returns the number of sampling points in X direction.
|
||||
*/
|
||||
template <typename DataType>
|
||||
void eval(const DataType& x, const DataType& y, DataType& result) const
|
||||
{
|
||||
if ( (!xExtrapolate_ && !appliesX(x)) ||
|
||||
(!yExtrapolate_ && !appliesY(y)) ) {
|
||||
std::ostringstream oss;
|
||||
oss << "Attempt to get undefined table value (" << x << ", " << y << ")";
|
||||
throw NumericalIssue(oss.str());
|
||||
};
|
||||
size_t numX() const
|
||||
{ return xPos_.size(); }
|
||||
|
||||
// bi-linear interpolation: first, calculate the x and y indices in the lookup
|
||||
// table ...
|
||||
const unsigned i = xSegmentIndex(x);
|
||||
const unsigned j = ySegmentIndex(y);
|
||||
|
||||
// bi-linear interpolation / extrapolation
|
||||
const DataType alpha = xToAlpha(x, i);
|
||||
const DataType beta = yToBeta(y, j);
|
||||
|
||||
const DataType s1 = valueAt(i, j) * (1.0 - beta) + valueAt(i, j + 1) * beta;
|
||||
const DataType s2 = valueAt(i + 1, j) * (1.0 - beta) + valueAt(i + 1, j + 1) * beta;
|
||||
|
||||
Valgrind::CheckDefined(s1);
|
||||
Valgrind::CheckDefined(s2);
|
||||
|
||||
// ... and combine them using the x position
|
||||
result = s1 * (1.0 - alpha) + s2 * alpha;
|
||||
Valgrind::CheckDefined(result);
|
||||
}
|
||||
|
||||
private:
|
||||
// the sampling points in the x-drection
|
||||
std::vector<Scalar> xPos_;
|
||||
// the sampling points in the y-drection
|
||||
std::vector<Scalar> yPos_;
|
||||
// data at the sampling points
|
||||
std::vector<std::vector<Scalar> > samples_;
|
||||
|
||||
bool xExtrapolate_ = false;
|
||||
bool yExtrapolate_ = false;
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in Y direction.
|
||||
*/
|
||||
size_t numY() const
|
||||
{ return yPos_.size(); }
|
||||
|
||||
/*!
|
||||
* \brief Returns the minimum of the X coordinate of the sampling points.
|
||||
@ -158,83 +136,126 @@ private:
|
||||
Scalar valueAt(size_t i, size_t j) const
|
||||
{ return samples_[i][j]; }
|
||||
|
||||
/*!
|
||||
* \brief Returns true if a coordinate lies in the tabulated range
|
||||
*/
|
||||
template <class Evaluation>
|
||||
bool applies(const Evaluation& x, const Evaluation& y) const
|
||||
{ return appliesX(x) && appliesY(y); }
|
||||
|
||||
/*!
|
||||
* \brief Returns true if a coordinate lies in the tabulated range on the x direction
|
||||
*/
|
||||
template <class Evaluation>
|
||||
bool appliesX(const Evaluation& x) const
|
||||
{
|
||||
if (x < xMin() || xMax() < x) {
|
||||
return false;
|
||||
} else {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
{ return xMin() <= x && x <= xMax(); }
|
||||
|
||||
/*!
|
||||
* \brief Returns true if a coordinate lies in the tabulated range on the y direction
|
||||
*/
|
||||
template <class Evaluation>
|
||||
bool appliesY(const Evaluation& y) const
|
||||
{ return yMin() <= y && y <= yMax(); }
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Evaluate the function at a given (x,y) position.
|
||||
*
|
||||
* If this method is called for a value outside of the tabulated
|
||||
* range, and extrapolation is not allowed in the corresponding direction,
|
||||
* a \c Opm::NumericalIssue exception is thrown.
|
||||
*/
|
||||
template <typename Evaluation>
|
||||
void eval(const Evaluation& x, const Evaluation& y, Evaluation& result) const
|
||||
{
|
||||
if (y < yMin() || yMax() < y) {
|
||||
return false;
|
||||
} else {
|
||||
return true;
|
||||
}
|
||||
if ((!xExtrapolate_ && !appliesX(x)) || (!yExtrapolate_ && !appliesY(y))) {
|
||||
std::ostringstream oss;
|
||||
oss << "Attempt to get undefined table value (" << x << ", " << y << ")";
|
||||
throw NumericalIssue(oss.str());
|
||||
};
|
||||
|
||||
// bi-linear interpolation: first, calculate the x and y indices in the lookup
|
||||
// table ...
|
||||
const unsigned i = xSegmentIndex_(x);
|
||||
const unsigned j = ySegmentIndex_(y);
|
||||
|
||||
// bi-linear interpolation / extrapolation
|
||||
const Evaluation alpha = xToAlpha(x, i);
|
||||
const Evaluation beta = yToBeta(y, j);
|
||||
|
||||
const Evaluation s1 = valueAt(i, j) * (1.0 - beta) + valueAt(i, j + 1) * beta;
|
||||
const Evaluation s2 = valueAt(i + 1, j) * (1.0 - beta) + valueAt(i + 1, j + 1) * beta;
|
||||
|
||||
Valgrind::CheckDefined(s1);
|
||||
Valgrind::CheckDefined(s2);
|
||||
|
||||
// ... and combine them using the x position
|
||||
result = s1 * (1.0 - alpha) + s2 * alpha;
|
||||
Valgrind::CheckDefined(result);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in X direction.
|
||||
*/
|
||||
size_t numX() const
|
||||
{ return xPos_.size(); }
|
||||
private:
|
||||
// the sampling points in the x-drection
|
||||
std::vector<Scalar> xPos_;
|
||||
// the sampling points in the y-drection
|
||||
std::vector<Scalar> yPos_;
|
||||
// data at the sampling points
|
||||
std::vector<std::vector<Scalar> > samples_;
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in Y direction.
|
||||
*/
|
||||
size_t numY() const
|
||||
{ return yPos_.size(); }
|
||||
bool xExtrapolate_ = false;
|
||||
bool yExtrapolate_ = false;
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the x-axis.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
unsigned xSegmentIndex(const Evaluation& x) const
|
||||
unsigned xSegmentIndex_(const Evaluation& x) const
|
||||
{
|
||||
assert(xExtrapolate_ || appliesX(x) );
|
||||
|
||||
return segmentIndex(x, xPos_);
|
||||
return segmentIndex_(x, xPos_);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the y-axis.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
unsigned ySegmentIndex(const Evaluation& y) const
|
||||
unsigned ySegmentIndex_(const Evaluation& y) const
|
||||
{
|
||||
assert(yExtrapolate_ || appliesY(y) );
|
||||
|
||||
return segmentIndex(y, yPos_);
|
||||
return segmentIndex_(y, yPos_);
|
||||
}
|
||||
|
||||
|
||||
template <class Evaluation>
|
||||
static unsigned segmentIndex(const Evaluation& y, const std::vector<Scalar>& pos)
|
||||
static unsigned segmentIndex_(const Evaluation& v, const std::vector<Scalar>& vPos)
|
||||
{
|
||||
const unsigned n = vPos.size();
|
||||
assert(n >= 2);
|
||||
|
||||
const unsigned num_pos = pos.size();
|
||||
assert(num_pos >= 2);
|
||||
|
||||
if (y <= pos.front() || num_pos == 2) {
|
||||
if (v <= vPos.front() || n == 2)
|
||||
return 0;
|
||||
} else if (y >= pos.back()) {
|
||||
return num_pos - 2;
|
||||
} else {
|
||||
assert(num_pos >= 3);
|
||||
else if (v >= vPos.back())
|
||||
return n - 2;
|
||||
|
||||
return --std::lower_bound(pos.begin(), pos.end(), y) - pos.begin();
|
||||
assert(n > 2 && v > vPos.front() && v < vPos.back());
|
||||
|
||||
// bisection. this assumes that the vPos array is strictly mononically
|
||||
// increasing.
|
||||
size_t lowerIdx = 0;
|
||||
size_t upperIdx = vPos.size() - 1;
|
||||
while (lowerIdx + 1 < upperIdx) {
|
||||
size_t pivotIdx = (lowerIdx + upperIdx) / 2;
|
||||
if (v < vPos[pivotIdx])
|
||||
upperIdx = pivotIdx;
|
||||
else
|
||||
lowerIdx = pivotIdx;
|
||||
}
|
||||
|
||||
assert(vPos[lowerIdx] <= v);
|
||||
assert(v <= vPos[lowerIdx + 1]);
|
||||
return lowerIdx;
|
||||
}
|
||||
|
||||
/*!
|
||||
|
||||
@ -39,363 +39,362 @@
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
template <class ScalarT>
|
||||
struct Test
|
||||
{
|
||||
typedef ScalarT Scalar;
|
||||
typedef ScalarT Scalar;
|
||||
|
||||
static Scalar testFn1(Scalar x, Scalar /* y */)
|
||||
{ return x; }
|
||||
static Scalar testFn1(Scalar x, Scalar /* y */)
|
||||
{ return x; }
|
||||
|
||||
static Scalar testFn2(Scalar /* x */, Scalar y)
|
||||
{ return y; }
|
||||
static Scalar testFn2(Scalar /* x */, Scalar y)
|
||||
{ return y; }
|
||||
|
||||
static Scalar testFn3(Scalar x, Scalar y)
|
||||
{ return x*y; }
|
||||
static Scalar testFn3(Scalar x, Scalar y)
|
||||
{ return x*y; }
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformTabulated2DFunction<Scalar> >
|
||||
createUniformTabulatedFunction(Fn& f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
unsigned m = 50;
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformTabulated2DFunction<Scalar> >
|
||||
createUniformTabulatedFunction(Fn& f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
unsigned m = 50;
|
||||
|
||||
Scalar yMin = -1/2.0;
|
||||
Scalar yMax = 1/3.0;
|
||||
unsigned n = 40;
|
||||
Scalar yMin = -1/2.0;
|
||||
Scalar yMax = 1/3.0;
|
||||
unsigned n = 40;
|
||||
|
||||
auto tab = std::make_shared<Opm::UniformTabulated2DFunction<Scalar>>(
|
||||
xMin, xMax, m,
|
||||
yMin, yMax, n);
|
||||
for (unsigned i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
for (unsigned j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n - 1) * (yMax - yMin);
|
||||
tab->setSamplePoint(i, j, f(x, y));
|
||||
auto tab = std::make_shared<Opm::UniformTabulated2DFunction<Scalar>>(
|
||||
xMin, xMax, m,
|
||||
yMin, yMax, n);
|
||||
for (unsigned i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
for (unsigned j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n - 1) * (yMax - yMin);
|
||||
tab->setSamplePoint(i, j, f(x, y));
|
||||
}
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformXTabulated2DFunction<Scalar> >
|
||||
createUniformXTabulatedFunction(Fn& f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
unsigned m = 50;
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformXTabulated2DFunction<Scalar> >
|
||||
createUniformXTabulatedFunction(Fn& f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
unsigned m = 50;
|
||||
Scalar yMin = -1/2.0;
|
||||
Scalar yMax = 1/3.0;
|
||||
unsigned n = 40;
|
||||
|
||||
Scalar yMin = -1/2.0;
|
||||
Scalar yMax = 1/3.0;
|
||||
unsigned n = 40;
|
||||
|
||||
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>();
|
||||
for (unsigned i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
tab->appendXPos(x);
|
||||
for (unsigned j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin);
|
||||
tab->appendSamplePoint(i, y, f(x, y));
|
||||
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>();
|
||||
for (unsigned i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
tab->appendXPos(x);
|
||||
for (unsigned j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin);
|
||||
tab->appendSamplePoint(i, y, f(x, y));
|
||||
}
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformXTabulated2DFunction<Scalar> >
|
||||
createUniformXTabulatedFunction2(Fn& f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
Scalar m = 50;
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformXTabulated2DFunction<Scalar> >
|
||||
createUniformXTabulatedFunction2(Fn& f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
Scalar m = 50;
|
||||
Scalar yMin = - 4.0;
|
||||
Scalar yMax = 5.0;
|
||||
|
||||
Scalar yMin = - 4.0;
|
||||
Scalar yMax = 5.0;
|
||||
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>();
|
||||
for (unsigned i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
tab->appendXPos(x);
|
||||
|
||||
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>();
|
||||
for (unsigned i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
tab->appendXPos(x);
|
||||
Scalar n = i + 10;
|
||||
|
||||
Scalar n = i + 10;
|
||||
|
||||
for (unsigned j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin);
|
||||
tab->appendSamplePoint(i, y, f(x, y));
|
||||
for (unsigned j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin);
|
||||
tab->appendSamplePoint(i, y, f(x, y));
|
||||
}
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::XYTabulated2DFunction<Scalar> >
|
||||
createXYTabulated2DFunction(Fn& f)
|
||||
{
|
||||
const Scalar xMin = -2.0;
|
||||
const Scalar xMax = 3.0;
|
||||
const unsigned m = 50;
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::XYTabulated2DFunction<Scalar> >
|
||||
createXYTabulated2DFunction(Fn& f)
|
||||
{
|
||||
const Scalar xMin = -2.0;
|
||||
const Scalar xMax = 3.0;
|
||||
const unsigned m = 50;
|
||||
const Scalar yMin = -1/2.0;
|
||||
const Scalar yMax = 1/3.0;
|
||||
const unsigned n = 40;
|
||||
|
||||
const Scalar yMin = -1/2.0;
|
||||
const Scalar yMax = 1/3.0;
|
||||
const unsigned n = 40;
|
||||
std::vector<Scalar> x_samples(m);
|
||||
std::vector<Scalar> y_samples(n);
|
||||
std::vector<std::vector<Scalar>> data(m, std::vector<Scalar>(n));
|
||||
|
||||
std::vector<Scalar> x_samples(m);
|
||||
std::vector<Scalar> y_samples(n);
|
||||
std::vector<std::vector<Scalar>> data(m, std::vector<Scalar>(n));
|
||||
for (unsigned i = 0; i < m; ++i) {
|
||||
x_samples[i] = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
|
||||
for (unsigned i = 0; i < m; ++i) {
|
||||
x_samples[i] = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
|
||||
for (unsigned j = 0; j < n; ++j) {
|
||||
y_samples[j] = yMin + Scalar(j)/(n -1) * (yMax - yMin);
|
||||
data[i][j] = f(x_samples[i], y_samples[j]);
|
||||
for (unsigned j = 0; j < n; ++j) {
|
||||
y_samples[j] = yMin + Scalar(j)/(n -1) * (yMax - yMin);
|
||||
data[i][j] = f(x_samples[i], y_samples[j]);
|
||||
}
|
||||
}
|
||||
|
||||
auto tab = std::make_shared<Opm::XYTabulated2DFunction<Scalar>>(x_samples, y_samples, data, true, true);
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
auto tab = std::make_shared<Opm::XYTabulated2DFunction<Scalar>>(x_samples, y_samples, data, true, true);
|
||||
template <class Fn, class Table>
|
||||
bool compareTableWithAnalyticFn(const Table& table,
|
||||
Scalar xMin,
|
||||
Scalar xMax,
|
||||
unsigned numX,
|
||||
|
||||
return tab;
|
||||
}
|
||||
Scalar yMin,
|
||||
Scalar yMax,
|
||||
unsigned numY,
|
||||
|
||||
template <class Fn, class Table>
|
||||
bool compareTableWithAnalyticFn(const Table& table,
|
||||
Scalar xMin,
|
||||
Scalar xMax,
|
||||
unsigned numX,
|
||||
Fn& f,
|
||||
Scalar tolerance = 1e-8)
|
||||
{
|
||||
// make sure that the tabulated function evaluates to the same thing as the analytic
|
||||
// one (modulo tolerance)
|
||||
for (unsigned i = 1; i <= numX; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/numX*(xMax - xMin);
|
||||
|
||||
Scalar yMin,
|
||||
Scalar yMax,
|
||||
unsigned numY,
|
||||
for (unsigned j = 0; j < numY; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/numY*(yMax - yMin);
|
||||
if (std::abs(table->eval(x, y) - f(x, y)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<<x<<","<<y<<") != f("<<x<<","<<y<<"): " << table->eval(x,y) << " != " << f(x,y) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Fn& f,
|
||||
Scalar tolerance = 1e-8)
|
||||
{
|
||||
// make sure that the tabulated function evaluates to the same thing as the analytic
|
||||
// one (modulo tolerance)
|
||||
for (unsigned i = 1; i <= numX; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/numX*(xMax - xMin);
|
||||
return true;
|
||||
}
|
||||
|
||||
for (unsigned j = 0; j < numY; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/numY*(yMax - yMin);
|
||||
if (std::abs(table->eval(x, y) - f(x, y)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<<x<<","<<y<<") != f("<<x<<","<<y<<"): " << table->eval(x,y) << " != " << f(x,y) << "\n";
|
||||
template <class Fn, class Table>
|
||||
bool compareTableWithAnalyticFn2(const Table& table,
|
||||
const Scalar xMin,
|
||||
const Scalar xMax,
|
||||
unsigned numX,
|
||||
const Scalar yMin,
|
||||
const Scalar yMax,
|
||||
unsigned numY,
|
||||
|
||||
Fn& f,
|
||||
Scalar tolerance = 1e-8)
|
||||
{
|
||||
// make sure that the tabulated function evaluates to the same thing as the analytic
|
||||
// one (modulo tolerance)
|
||||
for (unsigned i = 1; i <= numX; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/numX*(xMax - xMin);
|
||||
|
||||
for (unsigned j = 0; j < numY; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/numY*(yMax - yMin);
|
||||
Scalar result;
|
||||
table->eval(x, y, result);
|
||||
if (std::abs(result - f(x, y)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<<x<<","<<y<<") != f("<<x<<","<<y<<"): " << result << " != " << f(x,y) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class UniformTablePtr, class UniformXTablePtr, class Fn>
|
||||
bool compareTables(const UniformTablePtr uTable,
|
||||
const UniformXTablePtr uXTable,
|
||||
Fn& f,
|
||||
Scalar tolerance = 1e-8)
|
||||
{
|
||||
// make sure the uniform and the non-uniform tables exhibit the same dimensions
|
||||
if (std::abs(uTable->xMin() - uXTable->xMin()) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMin() != uXTable->xMin(): " << uTable->xMin() << " != " << uXTable->xMin() << "\n";
|
||||
return false;
|
||||
}
|
||||
if (std::abs(uTable->xMax() - uXTable->xMax()) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMax() != uXTable->xMax(): " << uTable->xMax() << " != " << uXTable->xMax() << "\n";
|
||||
return false;
|
||||
}
|
||||
if (uTable->numX() != uXTable->numX()) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numX() != uXTable->numX(): " << uTable->numX() << " != " << uXTable->numX() << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
for (unsigned i = 0; i < uTable->numX(); ++i) {
|
||||
if (std::abs(uTable->yMin() - uXTable->yMin(i)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMin() != uXTable->yMin("<<i<<"): " << uTable->yMin() << " != " << uXTable->yMin(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
if (std::abs(uTable->yMax() - uXTable->yMax(i)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMax() != uXTable->yMax("<<i<<"): " << uTable->yMax() << " != " << uXTable->yMax(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
if (uTable->numY() != uXTable->numY(i)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numY() != uXTable->numY("<<i<<"): " << uTable->numY() << " != " << uXTable->numY(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class Fn, class Table>
|
||||
bool compareTableWithAnalyticFn2(const Table& table,
|
||||
const Scalar xMin,
|
||||
const Scalar xMax,
|
||||
unsigned numX,
|
||||
const Scalar yMin,
|
||||
const Scalar yMax,
|
||||
unsigned numY,
|
||||
|
||||
Fn& f,
|
||||
Scalar tolerance = 1e-8)
|
||||
{
|
||||
// make sure that the tabulated function evaluates to the same thing as the analytic
|
||||
// one (modulo tolerance)
|
||||
for (unsigned i = 1; i <= numX; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/numX*(xMax - xMin);
|
||||
|
||||
for (unsigned j = 0; j < numY; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/numY*(yMax - yMin);
|
||||
Scalar result;
|
||||
table->eval(x, y, result);
|
||||
if (std::abs(result - f(x, y)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<<x<<","<<y<<") != f("<<x<<","<<y<<"): " << result << " != " << f(x,y) << "\n";
|
||||
// make sure that the x and y values are identical
|
||||
for (unsigned i = 0; i < uTable->numX(); ++i) {
|
||||
if (std::abs(uTable->iToX(i) - uXTable->iToX(i)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->iToX("<<i<<") != uXTable->iToX("<<i<<"): " << uTable->iToX(i) << " != " << uXTable->iToX(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class UniformTablePtr, class UniformXTablePtr, class Fn>
|
||||
bool compareTables(const UniformTablePtr uTable,
|
||||
const UniformXTablePtr uXTable,
|
||||
Fn& f,
|
||||
Scalar tolerance = 1e-8)
|
||||
{
|
||||
// make sure the uniform and the non-uniform tables exhibit the same dimensions
|
||||
if (std::abs(uTable->xMin() - uXTable->xMin()) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMin() != uXTable->xMin(): " << uTable->xMin() << " != " << uXTable->xMin() << "\n";
|
||||
return false;
|
||||
}
|
||||
if (std::abs(uTable->xMax() - uXTable->xMax()) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMax() != uXTable->xMax(): " << uTable->xMax() << " != " << uXTable->xMax() << "\n";
|
||||
return false;
|
||||
}
|
||||
if (uTable->numX() != uXTable->numX()) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numX() != uXTable->numX(): " << uTable->numX() << " != " << uXTable->numX() << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
for (unsigned i = 0; i < uTable->numX(); ++i) {
|
||||
if (std::abs(uTable->yMin() - uXTable->yMin(i)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMin() != uXTable->yMin("<<i<<"): " << uTable->yMin() << " != " << uXTable->yMin(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
if (std::abs(uTable->yMax() - uXTable->yMax(i)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMax() != uXTable->yMax("<<i<<"): " << uTable->yMax() << " != " << uXTable->yMax(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
if (uTable->numY() != uXTable->numY(i)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numY() != uXTable->numY("<<i<<"): " << uTable->numY() << " != " << uXTable->numY(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
// make sure that the x and y values are identical
|
||||
for (unsigned i = 0; i < uTable->numX(); ++i) {
|
||||
if (std::abs(uTable->iToX(i) - uXTable->iToX(i)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->iToX("<<i<<") != uXTable->iToX("<<i<<"): " << uTable->iToX(i) << " != " << uXTable->iToX(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
for (unsigned j = 0; j < uTable->numY(); ++j) {
|
||||
if (std::abs(uTable->jToY(j) - uXTable->jToY(i, j)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->jToY("<<j<<") != uXTable->jToY("<<i<<","<<j<<"): " << uTable->jToY(i) << " != " << uXTable->jToY(i, j) << "\n";
|
||||
return false;
|
||||
for (unsigned j = 0; j < uTable->numY(); ++j) {
|
||||
if (std::abs(uTable->jToY(j) - uXTable->jToY(i, j)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->jToY("<<j<<") != uXTable->jToY("<<i<<","<<j<<"): " << uTable->jToY(i) << " != " << uXTable->jToY(i, j) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// check that the appicable range is correct. Note that due to rounding errors it is
|
||||
// undefined whether the table applies to the boundary of the tabulated domain or not
|
||||
Scalar xMin = uTable->xMin();
|
||||
Scalar yMin = uTable->yMin();
|
||||
Scalar xMax = uTable->xMax();
|
||||
Scalar yMax = uTable->yMax();
|
||||
// check that the appicable range is correct. Note that due to rounding errors it is
|
||||
// undefined whether the table applies to the boundary of the tabulated domain or not
|
||||
Scalar xMin = uTable->xMin();
|
||||
Scalar yMin = uTable->yMin();
|
||||
Scalar xMax = uTable->xMax();
|
||||
Scalar yMax = uTable->yMax();
|
||||
|
||||
Scalar x = xMin - tolerance;
|
||||
Scalar y = yMin - tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
Scalar x = xMin - tolerance;
|
||||
Scalar y = yMin - tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMin - tolerance;
|
||||
y = yMin + tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
x = xMin - tolerance;
|
||||
y = yMin + tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMin + tolerance;
|
||||
y = yMin - tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
x = xMin + tolerance;
|
||||
y = yMin - tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMin + tolerance;
|
||||
y = yMin + tolerance;
|
||||
if (!uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (!uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
x = xMin + tolerance;
|
||||
y = yMin + tolerance;
|
||||
if (!uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (!uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMax + tolerance;
|
||||
y = yMax + tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
x = xMax + tolerance;
|
||||
y = yMax + tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMax - tolerance;
|
||||
y = yMax + tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
x = xMax - tolerance;
|
||||
y = yMax + tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMax + tolerance;
|
||||
y = yMax - tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
x = xMax + tolerance;
|
||||
y = yMax - tolerance;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMax - tolerance;
|
||||
y = yMax - tolerance;
|
||||
if (!uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (!uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
x = xMax - tolerance;
|
||||
y = yMax - tolerance;
|
||||
if (!uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (!uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
// make sure that the function values at the sampling points are identical and that
|
||||
// they correspond to the analytic function
|
||||
unsigned m2 = uTable->numX()*5;
|
||||
unsigned n2 = uTable->numY()*5;
|
||||
if (!compareTableWithAnalyticFn(uTable,
|
||||
xMin, xMax, m2,
|
||||
yMin, yMax, n2,
|
||||
f,
|
||||
tolerance))
|
||||
return false;
|
||||
if (!compareTableWithAnalyticFn(uXTable,
|
||||
xMin, xMax, m2,
|
||||
yMin, yMax, n2,
|
||||
f,
|
||||
tolerance))
|
||||
return false;
|
||||
// make sure that the function values at the sampling points are identical and that
|
||||
// they correspond to the analytic function
|
||||
unsigned m2 = uTable->numX()*5;
|
||||
unsigned n2 = uTable->numY()*5;
|
||||
if (!compareTableWithAnalyticFn(uTable,
|
||||
xMin, xMax, m2,
|
||||
yMin, yMax, n2,
|
||||
f,
|
||||
tolerance))
|
||||
return false;
|
||||
if (!compareTableWithAnalyticFn(uXTable,
|
||||
xMin, xMax, m2,
|
||||
yMin, yMax, n2,
|
||||
f,
|
||||
tolerance))
|
||||
return false;
|
||||
|
||||
return true;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
template <class TestType>
|
||||
inline int testAll( const typename TestType::Scalar tolerance = 1e-6 )
|
||||
inline int testAll(const typename TestType::Scalar tolerance = 1e-6)
|
||||
{
|
||||
TestType test;
|
||||
auto uniformTab = test.createUniformTabulatedFunction(TestType::testFn1);
|
||||
@ -415,10 +414,10 @@ inline int testAll( const typename TestType::Scalar tolerance = 1e-6 )
|
||||
|
||||
uniformXTab = test.createUniformXTabulatedFunction2(TestType::testFn3);
|
||||
if (!test.compareTableWithAnalyticFn(uniformXTab,
|
||||
-2.0, 3.0, 100,
|
||||
-4.0, 5.0, 100,
|
||||
TestType::testFn3,
|
||||
/*tolerance=*/1e-2))
|
||||
-2.0, 3.0, 100,
|
||||
-4.0, 5.0, 100,
|
||||
TestType::testFn3,
|
||||
/*tolerance=*/1e-2))
|
||||
return 1;
|
||||
|
||||
{
|
||||
@ -434,19 +433,19 @@ inline int testAll( const typename TestType::Scalar tolerance = 1e-6 )
|
||||
const unsigned n = 170;
|
||||
|
||||
// extrapolation and interpolation involved, the tolerance needs to be bigger
|
||||
const ScalarType temp_tolerance = 1000. * tolerance;
|
||||
const ScalarType tmpTolerance = 1000. * tolerance;
|
||||
|
||||
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn1, temp_tolerance))
|
||||
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn1, tmpTolerance))
|
||||
return 1;
|
||||
|
||||
xytab = test.createXYTabulated2DFunction(TestType::testFn2);
|
||||
|
||||
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn2, temp_tolerance))
|
||||
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn2, tmpTolerance))
|
||||
return 1;
|
||||
|
||||
xytab = test.createXYTabulated2DFunction(TestType::testFn3);
|
||||
|
||||
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn3, temp_tolerance))
|
||||
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn3, tmpTolerance))
|
||||
return 1;
|
||||
}
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user