Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects
Commit 3e798115 authored by ext-rogers_c's avatar ext-rogers_c
Browse files

Tests now all pass

parent d140879f
No related branches found
No related tags found
1 merge request!459Resolve "Unit Test TestSolvePolynomialQuadraticSmoothed disabled"
Hide whitespace changes
Inline Side-by-side
src/Classic/Fields/Interpolation/PPSolveFactory.cpp
+ 30
− 9
View file @ 3e798115
......@@ -226,7 +226,7 @@ void PPSolveFactory::getDerivPoints() {
isEqual = polyIndex[k] == derivIndices_m.back()[k];
}
if (isEqual) {
derivIndexByPower_m.push_back(index);
// derivIndexByPower_m.push_back(index);
break;
}
++index;
......@@ -256,25 +256,35 @@ void PPSolveFactory::getDerivPoints() {
}
void PPSolveFactory::getDerivs(const Mesh::Iterator& it) {
// calculate the derivatives near to polynomial at position indexed by it
#ifdef DEBUG
bool verbose = false;
if (verbose) {
std::cerr << "PPSolveFactory::GetDerivs at position ";
for (auto p: it.getPosition()) {std::cerr << p << " ";}
std::cerr << std::endl;
}
#endif
// derivatives go in derivValues_m
derivValues_m = std::vector< std::vector<double> >(derivPoints_m.size());
std::vector<double> itPos = it.getPosition();
int posDim = it.getState().size();
// get the outer layer of points
Mesh::Iterator end = it.getMesh()->end()-1;
Mesh::Iterator last = it.getMesh()->end()-1;
int deltaOrder = smoothingOrder_m - polyPatchOrder_m;
if (deltaOrder <= 0)
if (deltaOrder <= 0) { // no smoothing needed
return;
}
for (size_t i = 0; i < derivPoints_m.size(); ++i) {
Mesh::Iterator nearest = it;
for (int j = 0; j < posDim; ++j) {
bool inTheMesh = true;
for (int j = 0; j < posDim && inTheMesh; ++j) {
nearest[j] += derivOrigins_m[i][j];
if (nearest[j] > end[j]) {
nearest = it;
break;
}
inTheMesh = nearest[j] <= last[j];
}
if (polynomials_m[nearest.toInteger()] == NULL) {
if (!inTheMesh) {
derivValues_m[i] = std::vector<double>(polyDim_m, 0.);
} else {
MMatrix<double> coeffs =
......@@ -285,6 +295,17 @@ void PPSolveFactory::getDerivs(const Mesh::Iterator& it) {
derivValues_m[i][j] = values(j+1);
}
}
#ifdef DEBUG
if (verbose) {
std::cerr << "Point ";
for (auto p: derivPoints_m[i]) {std::cerr << p << " ";}
std::cerr << " values ";
for (auto v: derivValues_m[i]) {std::cerr << v << " ";}
std::cerr << " derivIndices ";
for (auto in: derivIndices_m[i]) {std::cerr << in << " ";}
std::cerr << std::endl;
}
#endif
}
}
......
tests/classic_src/Fields/Interpolation/PPSolveFactoryTest.cpp
+ 100
− 74
View file @ 3e798115
......@@ -99,7 +99,7 @@ TEST(PPSolveFactoryTest, TestNearbyPointsSquares) {
class PPSolveFactoryTestFixture : public ::testing::Test {
protected:
std::vector<std::vector<double> > values;
int np;
int np1, np2, np3;
// 3D
SquarePolynomialVector ref;
ThreeDGrid* grid;
......@@ -118,10 +118,12 @@ class PPSolveFactoryTestFixture : public ::testing::Test {
refCoeffs = MMatrix<double>(2, 27, &data[0]);
ref = SquarePolynomialVector(3, refCoeffs);
// we get some data
np = 6;
np1 = 6;
np2 = 7;
np3 = 8;
// choose the grid so that a midpoint falls at 0/0/0; we want to compare it
// later on...
grid = new ThreeDGrid(1., 2., 3., -2.5, -5.0, -7.5, np, np, np);
grid = new ThreeDGrid(1., 2., 3., -2.5, -5.0, -7.5, np1, np2, np3);
for (Mesh::Iterator it = grid->begin(); it < grid->end(); ++it) {
std::vector<double> value(2);
ref.F(&it.getPosition()[0], &value[0]);
......@@ -137,14 +139,11 @@ class PPSolveFactoryTestFixture : public ::testing::Test {
ref2D = SquarePolynomialVector(3, refCoeffs2D);
std::vector<double> start = {-1., -2.};
std::vector<double> step = {2., 4.};
std::vector<int> nsteps = {np, np};
std::vector<int> nsteps = {np1, np2};
grid2D = new NDGrid(2, &nsteps[0], &step[0], &start[0]);
for (Mesh::Iterator it = grid2D->begin(); it < grid2D->end(); ++it) {
std::vector<double> value2D(2);
ref2D.F(&it.getPosition()[0], &value2D[0]);
std::cerr << "values2D " << it.getPosition()[0] << " "
<< it.getPosition()[1] << " ** "
<< value2D[0] << " " << value2D[1] << std::endl;
values2D.push_back(value2D);
}
}
......@@ -198,7 +197,7 @@ TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadratic) {
}
// now check that the values in each polynomial are correct
ThreeDGrid testGrid(1./4., 2./4., 3./4., -1., -2., -3., np*4-1, np*4-1, np*4-1);
ThreeDGrid testGrid(1./4., 2./4., 3./4., -1., -2., -3., np1*4-1, np2*4-1, np3*4-1);
for (Mesh::Iterator it = testGrid.begin(); it < testGrid.end(); ++it) {
std::vector<double> refValue(2);
std::vector<double> testValue(2);
......@@ -216,7 +215,9 @@ double f(double x) {
}
TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadratic1DSmoothed) {
std::vector<std::vector<double> > gridPoints = {{0., 1., 2., 3., 4., 5., 6.}};
OpalTestUtilities::SilenceTest silencer;
std::vector<std::vector<double> > gridPoints = {{-1., 1., 3., 5., 7., 9., 11.}};
std::vector<std::vector<double> > values1D(gridPoints[0].size(), std::vector<double>(1, 0.));
for (size_t i = 0; i < gridPoints[0].size(); ++i) {
values1D[i][0] = f(gridPoints[0][i]);
......@@ -230,7 +231,7 @@ TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadratic1DSmoothed) {
patch = fac->solve();
delete fac;
} catch (GeneralClassicException& exc) {
std::cerr << exc.what() << " " << exc.where() << std::endl;
std::cout << exc.what() << " " << exc.where() << std::endl;
throw;
}
// Check that we correctly generate polynomials with values that match f(x)
......@@ -260,9 +261,9 @@ TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadratic1DSmoothed) {
}
poly->F(&end[0], &values[2][0]);
derivPoly.F(&end[0], &values[3][0]);
std::cerr << "At pos: " << position[0] << " ";
std::cerr << "Values: " << values[0][0] << " " << values[2][0] << " ";
std::cerr << "Derivs: " << values[1][0] << " " << values[3][0] << std::endl;
std::cout << "At pos: " << position[0] << " ";
std::cout << "Values: " << values[0][0] << " " << values[2][0] << " ";
std::cout << "Derivs: " << values[1][0] << " " << values[3][0] << std::endl;
}
EXPECT_NEAR(values[3][0], 0., 1e-6); // derivative 0.0 at boundary
}
......@@ -270,6 +271,7 @@ TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadratic1DSmoothed) {
TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadratic2DSmoothed) {
//OpalTestUtilities::SilenceTest silencer;
Mesh* mesh = grid2D->clone();
//Mesh* dual = grid2D->dual();
PPSolveFactory* fac = NULL;
PolynomialPatch* patch = NULL;
try {
......@@ -280,69 +282,93 @@ TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadratic2DSmoothed) {
std::cerr << exc.what() << " " << exc.where() << std::endl;
throw;
}
// first check that the polynomial at 0, 0, 0 is the same as ref
std::vector<double> zero({0., 0.});
SquarePolynomialVector* test = patch->getPolynomialVector(&zero[0]);
MMatrix<double> testCoeffs = test->GetCoefficientsAsMatrix();
// check values match at grid points
std::vector<double> point({0, 0}), testValue(2), refValue(2);
test->F(&point[0], &testValue[0]);
ref2D.F(&point[0], &refValue[0]);
EXPECT_NEAR(testValue[0], refValue[0], 1e-12);
EXPECT_NEAR(testValue[1], refValue[1], 1e-12);
std::cout << "Testing at " << point[0] << " " << point[1] << std::endl;
std::cout << "Gives " << testValue[0] << " " << testValue[1] << std::endl;
std::cout << "Compared to " << refValue[0] << " " << refValue[1] << std::endl;
std::cout << "Ref a" << std::endl;
std::cout << refCoeffs2D << std::endl;
std::cout << "Test a" << std::endl;
std::cout << testCoeffs << std::endl;
return;
// check derivatives match at grid points
std::vector<double> d11RefPosition({2., 4.});
std::vector<int> d11Index({1, 1});
// get derivatives of neighbouring polynomials; test is at 0, 0 and ref is
// at 2, 4; so diagonal, the 11 derivatives should match (and also 01, 10
// but I don't bother testing those)
SquarePolynomialVector d11Test = test->Deriv(&d11Index[0]);
SquarePolynomialVector d11Ref =
patch->getPolynomialVector(&d11RefPosition[0])->Deriv(&d11Index[0]);
// position in local coordinate system of the neighbouring polynomials
std::vector<double> d11RefPoint({-1., -2.}), d11TestPoint({1., 2.});
std::vector<double> d11RefValue(2, 0.), d11TestValue(2, 0);
d11Ref.F(&d11RefPoint[0], &d11RefValue[0]);
d11Test.F(&d11TestPoint[0], &d11TestValue[0]);
EXPECT_NEAR(d11RefValue[0], d11TestValue[0], 1e-12);
EXPECT_NEAR(d11RefValue[1], d11TestValue[1], 1e-12);
std::cout << "Ref" << std::endl;
std::cout << refCoeffs2D << std::endl;
std::cout << "Test" << std::endl;
std::cout << testCoeffs << std::endl;
return;
ASSERT_EQ(testCoeffs.num_row(), refCoeffs2D.num_row());
ASSERT_EQ(testCoeffs.num_col(), refCoeffs2D.num_col());
for (size_t i = 0; i < testCoeffs.num_row(); ++i) {
for (size_t j = 0; j < testCoeffs.num_row(); ++j) {
EXPECT_NEAR(testCoeffs(i+1, j+1), refCoeffs2D(i+1, j+1), 1e-6)
<< "col " << i << " row " << j;
}
// first check that the fitted values (0th derivative) match input values on
// grid points
for (Mesh::Iterator it = mesh->begin(); it != mesh->end(); it++) {
std::vector<double> pos = it.getPosition();
std::vector<double> value1(2);
patch->function(&pos[0], &value1[0]);
std::vector<double> value2 = values2D[it.toInteger()];
EXPECT_NEAR(value1[0], value2[0], 1e-6);
EXPECT_NEAR(value1[1], value2[1], 1e-6);
}
return;
// now check that the values returned by each polynomial are correct
ThreeDGrid testGrid(1./4., 2./4., 3./4., -1., -2., -3., np*4-1, np*4-1, np*4-1);
for (Mesh::Iterator it = testGrid.begin(); it < testGrid.end(); ++it) {
std::vector<double> refValue(2);
std::vector<double> testValue(2);
ref.F(&it.getPosition()[0], &refValue[0]);
patch->function(&it.getPosition()[0], &testValue[0]);
for (size_t i = 0; i < 2; ++i) {
EXPECT_NEAR(refValue[i], testValue[i], 1e-6) << std::endl << it;
// now check that appropriate derivatives match neighbours on grid points:
NDGrid* dual = dynamic_cast<NDGrid*>(mesh->dual());
// grid size:
std::vector<double> step = {dual->coord(2, 0)-dual->coord(1, 0), dual->coord(2, 1)-dual->coord(1, 1)};
// indexes the derivative at each point:
std::vector< std::vector< std::vector<int> > > derivIndexVector = {
{{1, 0}},
{{0, 1}},
{{1, 0}, {0, 1}, {1, 1}}
};
// Pointer to the neighbouring polynomial in the grid
std::vector< std::vector<int> > deltaIndex = {{1, 0}, {0, 1}, {1, 1}};
// vector from the current polynomial to the point at which the derivatives
// are supposed to match with neighbour
std::vector< std::vector<double> > forwardDeltaVector = {
{+step[0]/2., -step[1]/2.},
{-step[0]/2., +step[1]/2.},
{+step[0]/2., +step[1]/2.},
};
// vector from the neighbouring polynomial to the point at which the
// derivatives are supposed to match with current polynomial
std::vector< std::vector<double> > backwardDeltaVector = {
{-step[0]/2., -step[1]/2.},
{-step[0]/2., -step[1]/2.},
{-step[0]/2., -step[1]/2.},
};
// iterate over all mesh points and calculate derivatives; check that
// derivatives match, as they are supposed to, with neighbour
for (Mesh::Iterator it00 = dual->begin(); it00 != dual->end(); it00++) {
bool verbose = false;
if (verbose) {std::cout << "\n";}
for (size_t deltaI = 0; deltaI < deltaIndex.size() && deltaI < derivIndexVector.size(); deltaI++) {
for (size_t derivI = 0; derivI < derivIndexVector[deltaI].size(); derivI++) {
std::vector<int> thisDeltaVector = deltaIndex[deltaI];
std::vector<int> thisDerivVector = derivIndexVector[deltaI][derivI];
std::vector<double> pos00 = it00.getPosition();
Mesh::Iterator itDelta(it00);
Mesh::Iterator boundCheck(it00);
boundCheck[0] += 1;
boundCheck[1] += 1;
itDelta[0] += thisDeltaVector[0];
itDelta[1] += thisDeltaVector[1];
std::vector<double> pos1(2, 0.0), pos2(2, 0.0), posDelta(2, 0.0);
std::vector<double> derivCalc1(2, 0.0);
if (!itDelta.isOutOfBounds()) {
posDelta = itDelta.getPosition();
SquarePolynomialVector* polyDelta = patch->getPolynomials()[itDelta.toInteger()];
SquarePolynomialVector derivDelta = polyDelta->Deriv(&thisDerivVector[0]);
pos1 = backwardDeltaVector[deltaI];
derivDelta.F(&pos1[0], &derivCalc1[0]);
}
SquarePolynomialVector* poly00 = patch->getPolynomials()[it00.toInteger()];
SquarePolynomialVector deriv = poly00->Deriv(&thisDerivVector[0]);
std::vector<double> derivCalc2(2, 0.0);
pos2 = forwardDeltaVector[deltaI];
deriv.F(&pos2[0], &derivCalc2[0]);
if (verbose) {
std::cout << "Delta " << thisDeltaVector[0] << " " << thisDeltaVector[1];
std::cout << " deriv " << thisDerivVector[0] << " " << thisDerivVector[1] << std::endl;
std::cout << " OOB: " << boundCheck.isOutOfBounds() << " " << itDelta.isOutOfBounds();
std::cout << " here: " << it00[0] << " " << it00[1];
std::cout << " chck: " << boundCheck[0] << " " << boundCheck[1];
std::cout << " dlt: " << itDelta[0] << " " << itDelta[1];
std::cout << " last: " << (dual->end()-1)[0] << " " << (dual->end()-1)[1] << std::endl;
std::cout << " Offset from this " << pos00[0] << " " << pos00[1]
<< " by "<< pos2[0] << " " << pos2[1] << std::endl;
std::cout << " Offset from neighbour " << posDelta[0] << " " << posDelta[1]
<< " by "<< pos1[0] << " " << pos1[1] << std::endl;
std::cout << " deriv this " << derivCalc2[0] << " " << derivCalc2[1] << std::endl;
std::cout << " deriv neighbour " << derivCalc1[0] << " " << derivCalc1[1] << std::endl;
}
EXPECT_NEAR(derivCalc1[0], derivCalc2[0], 1e-6);
EXPECT_NEAR(derivCalc1[1], derivCalc2[1], 1e-6);
}
}
}
}
......
tests/classic_src/Fields/Interpolation/PPSolveFactoryTestPlot.py 0 → 100644
+ 150
− 0
View file @ 3e798115
import glob
import matplotlib.pyplot
class PlotInterpolation(object):
"""
Make plots to support plotting of the PPSolveFactoryTest output
"""
def __init__(self, file_name_list):
"""
Initialise the plotter - file_name_list is the list of files that will
be plotted
"""
self.file_name_list = file_name_list
self.data = {}
self.plots = {}
self.setup_plots()
def setup_plots(self):
"""
Set up empty plots
"""
self.plots["fit_fig"] = matplotlib.pyplot.figure(0, (12.8, 9.6))
self.plots["fit_fig"].xlabel = "s"
self.plots["fit_fig"].ylabel = "value"
kwd = {"ylim":[-1, 1], "xlim":[-1.0, 20.0]}
ax = self.plots["fit_fig"].add_subplot(1, 1, 1, **kwd)
self.plots["residual_fig"] = matplotlib.pyplot.figure(1, (12.8, 9.6))
self.plots["residual_fig"].xlabel = "s"
self.plots["residual_fig"].ylabel = "residual"
kwd = {"yscale":"log", "ylim":[1e-6, 1000], "xlim":[-1.0, 25.0]}
ax = self.plots["residual_fig"].add_subplot(1, 1, 1, **kwd)
def get_fit_fig(self):
matplotlib.pyplot.figure(0)
return self.plots["fit_fig"]
def get_residual_fig(self):
matplotlib.pyplot.figure(1)
return self.plots["residual_fig"]
def finish_plots(self):
self.get_fit_fig()
matplotlib.pyplot.legend()
self.get_residual_fig()
matplotlib.pyplot.legend()
matplotlib.pyplot.show(block=False)
self.plots["fit_fig"].savefig("ppsolvefactorytest_fit.png")
self.plots["residual_fig"].savefig("ppsolvefactorytest_residual.png")
def load_file(self, file_name):
fin = open(file_name)
self.data = {}
self.data["title"] = fin.readline().rstrip()
smooth = self.data["title"].split("smooth")[1]
smooth = smooth.split("_")[1]
poly = self.data["title"].split("poly")[1]
poly = poly.split("_")[1]
self.data["smooth"] = int(smooth)
self.data["poly"] = int(poly)
columns = fin.readline()
columns.rstrip()
self.data["columns"] = columns.split()
for col in self.data["columns"]:
self.data[col] = []
for line in fin.readlines():
line.rstrip()
for i, a_word in enumerate(line.split()):
col = self.data["columns"][i]
self.data[col].append(float(a_word))
self.data["length"] = len(self.data[col])
self.get_s_data()
self.get_fit_errors()
def get_s_data(self):
s_list = [0.]
for i in range(self.data["length"]):
x, y, z = self.data["pos_x"][i], self.data["pos_y"][i], self.data["pos_z"][i]
if i == 0:
print("Stepping from", x, y, z, end=' ')
if i > 0:
s = ((x-x_old)**2+(y-y_old)**2+(z-z_old)**2)**0.5
s_list.append(s+s_list[-1])
x_old, y_old, z_old = x, y, z
print("to", x, y, z)
self.data["s_pos"] = s_list
def get_fit_errors(self):
for key in "x", "y", "z":
true = self.data["true_"+key]
fit = self.data["fit_"+key]
error = [abs(fit[i]-true[i]) for i in range(self.data["length"])]
self.data["residual_"+key] = error
def get_title_suffix(self):
return "t: "+ str(self.data["smooth"])+ " f: "+str(self.data["poly"])
def do_fit_plots(self):
self.get_fit_fig()
title = "Fitted - "+self.get_title_suffix()
abs_list = self.data["s_pos"]
ord_list = self.data["fit_x"]
matplotlib.pyplot.plot(abs_list, ord_list,
color=self.color(), marker=self.marker(),
fillstyle='none', label=title)
def color(self):
i = self.data['smooth']
c_list = ['red', 'blue', 'green', 'orange', 'magenta']
color = c_list[i%len(c_list)]
return color
def marker(self):
i = self.data['poly']
m_list = ['o', 'v', '^', 's', 'x']
marker = m_list[i%len(m_list)]
return marker
def do_truth_plot(self):
self.get_fit_fig()
abs_list = self.data["s_pos"]
ord_list = self.data["true_x"]
color = self.color()
matplotlib.pyplot.plot(abs_list, ord_list, label="True x")
def do_residual_plot(self):
self.get_residual_fig()
title = "Residual - "+self.get_title_suffix()
abs_list = self.data["s_pos"]
ord_list = self.data["residual_x"]
matplotlib.pyplot.semilogy(abs_list, ord_list,
color=self.color(), marker=self.marker(),
fillstyle='none', label=title)
def main(self):
for i, file_name in enumerate(sorted(self.file_name_list)):
self.load_file(file_name)
self.do_fit_plots()
if i == len(self.file_name_list)-1:
self.do_truth_plot()
self.do_residual_plot()
self.finish_plots()
def main():
plotter = PlotInterpolation(glob.glob("PPSolveFactoryTest_*"))
plotter.main()
input("Press <CR> to end")
if __name__ == "__main__":
main()
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment