Commit f8c716f8 authored by JasonPries's avatar JasonPries
Browse files

Complete instrumentation of Sketch::solve()

Uniform testing of Sketch::build() results
parent 8496e143
......@@ -21,4 +21,4 @@ include_directories(./src/Mesh/include/)
add_subdirectory(./src/)
add_subdirectory(./test/)
add_subdirectory(./test/)
\ No newline at end of file
......@@ -43,7 +43,7 @@ void Angle::update(Eigen::MatrixXd &J, Eigen::VectorXd &r) const {
r(EquationIndex) = scale * (dot - cos(M_PI * Dim / 180.0));
f = scale * (vx1 - dot * vx0) / d0;
J(EquationIndex, Line0->start()->x_index()) -= f; //TODO: Change get_index() and set_index() to size_t index() and size_t index(size_t) (overloaded)
J(EquationIndex, Line0->start()->x_index()) -= f;
J(EquationIndex, Line0->end()->x_index()) += f;
f = scale * (vx0 - dot * vx1) / d1;
......
......@@ -7,13 +7,9 @@ void CircularArc::update(Eigen::MatrixXd &J, Eigen::VectorXd &r) const {
Normalize linearized equations by multiplying through by radius()
*/
// #TODO: Clean this method
double dx, dy, rad;
// Jacobian, Residual
dx = Start->x() - Center->x();
dy = Start->y() - Center->y();
rad = sqrt(dx * dx + dy * dy);
double dx = Start->x() - Center->x();
double dy = Start->y() - Center->y();
double rad = sqrt(dx * dx + dy * dy);
r(EquationIndex) = radius() - rad;
......@@ -242,18 +238,12 @@ bool CircularArc::is_coincident(std::shared_ptr<Curve const> const &c) const {
if (cc.get() == nullptr) {
return false;
} else {
// #TODO: Extract input curve center, put the rest of the method in subroutine
// #TODO: Then, can call subroutine for implementation of rotation version
double xc = cc->Center->x();
double yc = cc->Center->y();
double rc = cc->radius();
double tol = radius() * FLT_EPSILON;
if (abs(Center->x() - xc) < tol && abs(Center->y() - yc) < tol && abs(radius() - rc) < tol) {
return true;
} else {
return false;
}
return (abs(Center->x() - xc) < tol && abs(Center->y() - yc) < tol && abs(radius() - rc) < tol);
}
}
......
......@@ -7,7 +7,6 @@
class CircularArc final : public Curve {
public:
using Curve::on_manifold;
using Curve::on_segment;
public:
......
......@@ -132,7 +132,8 @@ std::vector<std::shared_ptr<Contour>> Constellation::contours() {
if (find_closed_contour(contour_curves, orientation)) {
contours.push_back(std::make_shared<Contour>(contour_curves, orientation));
} else {
return std::vector<std::shared_ptr<Contour>>(); // TODO: Ugly multiple return points
std::cerr << "Failed to find closed contour" << std::endl;
break;
}
}
......
......@@ -136,7 +136,7 @@ MatchOrientation LineSegment::is_identical(double x0, double y0, double x1, doub
double xe = end()->x();
double ye = end()->y();
double tol = FLT_EPSILON * std::fmax(abs(xs - xe), abs(ys - ye)); // #TODO: L1 norm is more efficient tolerance strategy
double tol = FLT_EPSILON * std::fmax(abs(xs - xe), abs(ys - ye));
if (abs(xs - x0) < tol && abs(ys - y0) < tol && abs(xe - x1) < tol && abs(ye - y1) < tol) {
return MatchOrientation::Forward;
......@@ -153,11 +153,7 @@ bool LineSegment::is_coincident(std::shared_ptr<Curve const> const &c) const {
if (l.get() == nullptr) {
return false;
} else {
if (on_manifold(c->start()) && on_manifold(c->end())) {
return true;
} else {
return false;
}
return (on_manifold(c->start()) && on_manifold(c->end()));
}
}
......
......@@ -7,7 +7,6 @@
class LineSegment final : public Curve {
public:
using Curve::on_manifold;
using Curve::on_segment;
public:
......
......@@ -4,7 +4,6 @@
RotateCopy::RotateCopy(std::vector<std::shared_ptr<Curve const>> input, std::shared_ptr<Vertex const> center, double angle, size_t copies, bool remove_internal) {
// Creates rotated copies of the input curves about an vertex
// TODO: Three Groups (?): Leading Verticies, Lagging Verticies, Internal Verticies
// TODO: Check for complete elimination of leading/lagging curves
// Assign properties
......
......@@ -15,152 +15,137 @@
#include <iostream>
double Sketch::solve() {
// TODO: Time to finish the bells and whistles here
// TODO: Use sparse matrix representation for Jacobian, (this is lower priority since typical uses will not produce very large matrices)
// TODO: Try extracting sub matrix and solving (Certain constraints are known to eliminate variables (Radius, Fixation), if the matrix can be permuted so that a submatrix on the block diagonal is lower-triangular, then a subsystem of equations can be solved independently of the remaining equations. Moreover, if the jacobian can be permuted into a block lower triangular form, subsets of equations can be solved blockwise)
// TODO: Extract newton_update over Verticies, Curves, and Constraints into a private method
// TODO: Extract armijo_update into a private method
// TODO: Strictly decreasing armijo_update from armijo_int = 0; (Typical expected case)
// TODO: Use variable_feature_size() and equation_feature_size()
if (NumVariables < NumEquations) {
std::cerr << "Variables = " << NumVariables << ", Equations = " << NumEquations << ". Number of variables is less than the number of equations. The sketch may be over-constrained." << std::endl;
}
Eigen::VectorXd r = Eigen::VectorXd::Zero(NumEquations);
// Matrix and vectors
Eigen::MatrixXd J = Eigen::MatrixXd::Zero(NumEquations, NumVariables);
Eigen::VectorXd r = Eigen::VectorXd::Zero(NumEquations);
Eigen::VectorXd lfs = Eigen::VectorXd::Zero(NumEquations);
// Calculate scale parameter for tolerances
// TODO: Sketch::scale();
double scale = 0.0;
if (Curves.size() > 0) {
for (size_t i = 0;i!= Curves.size();++i){
scale = std::max(scale, Curves[i]->length());
}
} else if (Verticies.size() > 0) {
double vmax = DBL_MIN;
double vmin = DBL_MAX;
for (size_t i = 0; i != Verticies.size(); ++i) {
vmax = std::max(vmax, Verticies[i]->x());
vmin = std::min(vmin, Verticies[i]->x());
vmax = std::max(vmax, Verticies[i]->y());
vmin = std::max(vmin, Verticies[i]->y());
}
// Tolerances and iteration limits
double dbl_tol{std::sqrt(NumEquations + 1.0) * DBL_EPSILON};
double flt_tol{std::sqrt(NumEquations + 1.0) * FLT_EPSILON};
if (vmax != vmin) {
scale = vmax - vmin;
} else {
scale = 1.0;
}
}
size_t double_bits = (CHAR_BIT * sizeof(double));
// Initial residual calculation
// TODO: Sketch::update_equation(J,r)
for (size_t i = 0; i != Verticies.size(); ++i) {
Verticies[i]->update(J, r);
}
for (size_t i = 0; i != Curves.size(); ++i) {
Curves[i]->update(J, r);
}
for (size_t i = 0; i != Constraints.size(); ++i) {
Constraints[i]->update(J, r);
}
update_linearization(J, r);
// Setup iteration
double res_norm{r.norm() / scale};
double res_tol{std::sqrt(NumEquations + 1.0) * DBL_EPSILON};
size_t double_bits = (CHAR_BIT * sizeof(double));
long prev_rank{std::min(NumEquations,NumVariables)};
double res_norm{r.norm()};
double res_tol{characteristic_length() * dbl_tol};
long prev_rank(std::min(NumEquations, NumVariables));
for (size_t j = 0; res_norm > res_tol && j != double_bits; ++j) {
Eigen::FullPivHouseholderQR<Eigen::MatrixXd> QR = J.fullPivHouseholderQr();
Eigen::VectorXd delta = QR.solve(r);
long iter_rank = QR.rank();
if (((J * delta) - r).norm() > res_tol + r.norm() * flt_tol) {
std::cerr << "Linearized sketch equation could not be solved with an average error per equation less than FLT_EPSILON. Some elements may be over-constrained." << std::endl;
}
long iter_rank = QR.rank();
if (iter_rank < prev_rank) {
std::cerr << "Rank = " << iter_rank << ", Equations = " << NumEquations << ", Variables = " << NumVariables << ", Iteration = " << j
<< ". The jacobian rank is less than the number of equations. Some sketch elements may be over-constrained." << std::endl;
// Randomly perturb the delta vector to prevent stagnation at inflection point
std::random_device rdev;
std::mt19937 rgen(rdev());
auto rdist = std::uniform_real_distribution<>(-DBL_EPSILON * scale, DBL_EPSILON * scale);
for (size_t i = 0; i != Variables.size(); ++i) {
delta(Variables[i]->get_index()) += rdist(rgen);
}
// perturb(delta, characteristic_length()); // perturb update in case of stagnation at inflection point
}
prev_rank = iter_rank;
// Update variables
// TODO: Sketch::update_variables(delta);
for (size_t i = 0; i < Variables.size(); ++i) {
Variables[i]->update(delta);
}
update_variables(delta);
update_linearization(J, r);
// TODO: Sketch::update_equation
J.setZero();
double iter_norm{r.norm()};
double iter_tol{res_tol};
if (iter_norm > res_norm) { // Line search
delta *= -0.5;
for (size_t i = 0; i != Verticies.size(); ++i) {
Verticies[i]->update(J, r);
}
size_t line_search_iter{0};
while (iter_norm > res_norm && iter_norm > iter_tol && line_search_iter != double_bits) {
update_variables(delta);
update_linearization(J, r);
for (size_t i = 0; i != Curves.size(); ++i) {
Curves[i]->update(J, r);
}
iter_norm = r.norm();
iter_tol = dbl_tol * characteristic_length();
if (iter_norm > res_norm) {
delta *= 0.5;
line_search_iter += 1;
for (size_t i = 0; i != Constraints.size(); ++i) {
Constraints[i]->update(J, r);
if (line_search_iter == double_bits) {
std::cerr << "Residual norm was not reduced for any step-size greater than 2^-(CHAR_BIT * sizeof(double))" << std::endl;
}
}
}
}
res_norm = iter_norm;
res_tol = iter_tol;
}
// Update Jacobian, residual, and take Armijo step
// TODO: Sketch::armijo(res_norm,res_tol);
double iter_norm{r.norm() / scale};
if (iter_norm > res_norm) {
delta *= -0.5;
return res_norm;
}
size_t armijo_iter{0};
while (iter_norm > res_norm && iter_norm > res_tol && res_norm > res_tol && armijo_iter != double_bits) {
// TODO: Sketch::update_variables(delta);
for (size_t i = 0; i < Variables.size(); ++i) {
Variables[i]->update(delta);
}
void Sketch::update_linearization(Eigen::MatrixXd &J, Eigen::VectorXd &r) const {
J.setZero();
// TODO: Sketch::update_equation
J.setZero();
for (size_t i = 0; i != Verticies.size(); ++i) {
Verticies[i]->update(J, r);
}
for (size_t i = 0; i != Verticies.size(); ++i) {
Verticies[i]->update(J, r);
}
for (size_t i = 0; i != Curves.size(); ++i) {
Curves[i]->update(J, r);
}
for (size_t i = 0; i != Curves.size(); ++i) {
Curves[i]->update(J, r);
}
for (size_t i = 0; i != Constraints.size(); ++i) {
Constraints[i]->update(J, r);
}
}
for (size_t i = 0; i != Constraints.size(); ++i) {
Constraints[i]->update(J, r);
}
void Sketch::update_variables(Eigen::VectorXd &delta) const {
for (size_t i = 0; i != Variables.size(); ++i) {
Variables[i]->update(delta);
}
}
iter_norm = r.norm() / scale;
void Sketch::perturb(Eigen::VectorXd &delta, double scale) const {
// Randomly perturb the delta vector to prevent stagnation at inflection point
std::random_device rdev;
std::mt19937 rgen(rdev());
std::uniform_real_distribution<double> rdist = std::uniform_real_distribution<double>(-DBL_EPSILON * scale, DBL_EPSILON * scale);
for (size_t i = 0; i != Variables.size(); ++i) {
size_t j = Variables[i]->get_index();
delta(j) += rdist(rgen) * scale;
}
}
if (iter_norm > res_norm) {
delta *= 0.5;
armijo_iter += 1;
double Sketch::characteristic_length() const {
double scale{0.0};
if (armijo_iter == double_bits) {
std::cerr << "Residual norm was not reduced for Armijo step-sizes greater than 2^-(CHAR_BIT * sizeof(double))" << std::endl;
}
}
}
if (Curves.size() > 0) {
for (size_t i = 0; i != Curves.size(); ++i) {
scale = std::max(scale, Curves[i]->length());
}
} else if (Verticies.size() > 0) {
double xmax = DBL_MIN;
double xmin = DBL_MAX;
double ymax = DBL_MIN;
double ymin = DBL_MAX;
for (size_t i = 0; i != Verticies.size(); ++i) {
xmax = std::max(xmax, Verticies[i]->x());
xmin = std::min(xmin, Verticies[i]->x());
ymax = std::max(ymax, Verticies[i]->y());
ymin = std::max(ymin, Verticies[i]->y());
}
scale = std::max(xmax - xmin, ymax-ymin);
if(scale == 0) {
scale = std::max(xmax,ymax);
}
res_norm = iter_norm;
}
return res_norm;
return scale;
}
bool Sketch::build() { // TOOD: Instrumentation in tests
bool Sketch::build() {
Constellation c = Constellation(this);
Boundary = c.boundary();
......
......@@ -8,6 +8,7 @@
#include <string>
#include <utility>
#include <vector>
#include <iostream>
#include <boost/filesystem.hpp>
......@@ -40,6 +41,8 @@ public:
bool build(); // TODO: Detailed return enum
double characteristic_length() const;
double solve(); // TODO: Detailed return enum
void add_element(std::shared_ptr<Constraint> c);
......@@ -52,6 +55,12 @@ public:
void add_parameter(std::shared_ptr<Variable> v) { add_parameter(Variables, v); };
void perturb(Eigen::VectorXd &delta, double scale) const;
void update_linearization(Eigen::MatrixXd &J, Eigen::VectorXd &r) const;
void update_variables(Eigen::VectorXd &delta) const;
template<SaveMethod S>
void save_as(std::string path, std::string file_name) const;
......
......@@ -15,7 +15,7 @@ public:
Vertex(std::shared_ptr<Vertex const> const &v) : X(v->X), Y(v->Y) {};
size_t set_equation_index(size_t i) override { //
size_t set_equation_index(size_t i) override {
EquationIndex = i;
return 0;
};
......
......@@ -314,8 +314,6 @@ TEST(Constraint, Distance_Vertex_Vertex) {
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, "Constraint__Distance_Vertex_Vertex");
EXPECT_NEAR(M_LN2, hypot(v0->x() - v1->x(), v0->y() - v1->y()), TOL * M_LN2);
}
......
......@@ -29,5 +29,6 @@ TEST(Curve, supremum) {
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
bool result = s.build();
ASSERT_TRUE(result);
}
\ No newline at end of file
......@@ -81,7 +81,8 @@ TEST(MirrorCopy, nonoverlapping) {
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
bool result = s.build();
ASSERT_TRUE(result);
// Run Tests
test_sketch_size(s, 10, 9, 6, 2);
......@@ -94,7 +95,8 @@ TEST(MirrorCopy, nonoverlapping) {
res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
result = s.build();
ASSERT_TRUE(result);
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, "Pattern__Mirror_nonoverlapping_trapezoid");
......@@ -130,7 +132,8 @@ TEST(MirrorCopy, overlapping) {
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
bool result = s.build();
ASSERT_TRUE(result);
// Run Tests
test_sketch_size(s, 6, 7, 4, 2 - remove_internal);
......@@ -143,7 +146,8 @@ TEST(MirrorCopy, overlapping) {
res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
result = s.build();
ASSERT_TRUE(result);
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("Mirror__overlapping_trapezoid") + std::to_string(remove_internal));
......@@ -182,7 +186,8 @@ TEST(MirrorCopy, multiple_overlapping) {
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
bool result = s.build();
ASSERT_TRUE(result);
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("Mirror__multiple_overlapping_0_") + std::to_string(remove_internal));
......@@ -198,7 +203,8 @@ TEST(MirrorCopy, multiple_overlapping) {
res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
result = s.build();
ASSERT_TRUE(result);
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("Mirror__multiple_overlapping_1_") + std::to_string(remove_internal));
......
......@@ -36,7 +36,8 @@ TEST(RotateCopy, nonoverlapping) {
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
bool result = s.build();
ASSERT_TRUE(result);
// Run Tests
test_sketch_size(s, 13, 12, 14, 3);
......@@ -51,7 +52,8 @@ TEST(RotateCopy, nonoverlapping) {
res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.build();
result = s.build();
ASSERT_TRUE(result);
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, "Rotate__non_overlapping_1");
......@@ -100,13 +102,18 @@ TEST(RotateCopy, overlapping) {
auto a1 = s.new_element<Angle>(l0, l1, a_deg);
// Test base sketch
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("Rotate__overlapping_base_0_") + std::to_string(remove_internal));
bool result = s.build();
ASSERT_TRUE(result);
test_sketch_size(s, 6, 5, 8, 1);
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("Rotate__overlapping_base_0_") + std::to_string(remove_internal));
// Create rotate copy
std::vector<std::shared_ptr<Curve const>> rvec = {l0, l1, c0, c1};
auto r0 = s.new_element<RotateCopy>(rvec, v0, 360.0 / N, N - 2, remove_internal);
......@@ -118,7 +125,6 @@ TEST(RotateCopy, overlapping) {
result = s.build();
ASSERT_TRUE(result);
// Run Tests
test_sketch_size(s, 10, 11, 12, 3 - 2 * remove_internal);
test_rotated_verticies(s, {1, 2, 3, 4}, v0, 360.0 / N, N - 2);
test_rotated_curves(s, {0, 1, 3, 4}, v0, 360.0 / N, N - 2);
......
......@@ -174,6 +174,7 @@ TEST(Star, Suite_2) {
TEST(Star, find_closed_contour_0) {
Sketch sketch;
auto v0 = sketch.new_element<Vertex>(0.0, 0.0);
auto v1 = sketch.new_element<Vertex>(1.0, 0.0);
auto v2 = sketch.new_element<Vertex>(2.0, 0.0);
......@@ -197,30 +198,26 @@ TEST(Star, find_closed_contour_0) {
auto l8 = sketch.new_element<LineSegment>(v4, v0);
// Manual contour creation
{
std::vector<Star> stars;
std::vector<Star> stars;
stars.push_back(Star{v0, &sketch});
stars.push_back(Star{v1, &sketch});
stars.push_back(Star{v4, &sketch});
stars.push_back(Star{v0, &sketch});
stars.push_back(Star{v1, &sketch});
stars.push_back(Star{v4, &sketch});
star_angle_sum_equals_2pi(stars);
}
star_angle_sum_equals_2pi(stars);
// Sketch internal contour parsing
{
double res_norm = sketch.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
EXPECT_TRUE(sketch.build());
EXPECT_EQ(sketch.size_contours(), 1);
EXPECT_TRUE(*sketch.contour(0) == *sketch.boundary());
auto contour = sketch.contour(0);
EXPECT_TRUE(contour->size() == 3);
EXPECT_TRUE(l0 == contour->curve(0) || l0 == contour->curve(1) || l0 == contour->curve(2));
EXPECT_TRUE(l5 == contour->curve(0) || l5 == contour->curve(1) || l5 == contour->curve(2));
EXPECT_TRUE(l8 == contour->curve(0) || l8 == contour->curve(1) || l8 == contour->curve(2));
}
double res_norm = sketch.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
ASSERT_TRUE(sketch.build());
EXPECT_EQ(sketch.size_contours(), 1);
EXPECT_TRUE(*sketch.contour(0) == *sketch.boundary());
auto contour = sketch.contour(0);
EXPECT_TRUE(contour->size() == 3);
EXPECT_TRUE(l0 == contour->curve(0) || l0 == contour->curve(1) || l0 == contour->curve(2));
EXPECT_TRUE(l5 == contour->curve(0) || l5 == contour->curve(1) || l5 == contour->curve(2));
EXPECT_TRUE(l8 == contour->curve(0) || l8 == contour->curve(1) || l8 == contour->curve(2));
}
TEST(Star, find_closed_contour_1) {
......@@ -240,29 +237,25 @@ TEST(Star, find_closed_contour_1) {
auto c1 = sketch.new_element<CircularArc>(v3, v1, v0, 1.0);
// Manual contour construction
{
std::vector<Star> stars;
stars.push_back(Star{v0, &sketch});
stars.push_back(Star{v1, &sketch});
stars.push_back(Star{v2, &sketch});
std::vector<Star> stars;
stars.push_back(Star{v0, &sketch});
stars.push_back(Star{v1, &sketch});
stars.push_back(Star{v2, &sketch});
star_angle_sum_equals_2pi(stars);
}
star_angle_sum_equals_2pi(stars);
// Sketch internal contour parsing
{
double res_norm = sketch.solve();