Sketch.cpp 8.9 KB
 Jason Pries committed Dec 13, 2016 1 ``````#include `````` JasonPries committed Feb 27, 2017 2 ``````#include `````` Jason Pries committed Dec 13, 2016 3 4 5 `````` #include "Sketch.h" #include "Vertex.h" `````` JasonPries committed Feb 27, 2017 6 7 ``````#include "CircularArc.h" #include "LineSegment.h" `````` Jason Pries committed Dec 13, 2016 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 ``````#include "Curve.h" #include "Constraint.h" #include "Pattern.h" #include "Constellation.h" #include "Contour.h" #include "Star.h" #include "Branch.h" #include "Variable.h" #include "doublen.h" #include double Sketch::solve() { // 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: Use variable_feature_size() and equation_feature_size() // 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); `````` Jason Pries committed Oct 11, 2016 29 `````` `````` Jason Pries committed Dec 13, 2016 30 31 32 `````` // Tolerances and iteration limits double dbl_tol{std::sqrt(NumEquations + 1.0) * DBL_EPSILON}; double flt_tol{std::sqrt(NumEquations + 1.0) * FLT_EPSILON}; `````` Jason Pries committed Oct 11, 2016 33 `````` `````` Jason Pries committed Dec 13, 2016 34 `````` size_t double_bits = (CHAR_BIT * sizeof(double)); `````` Jason Pries committed Oct 11, 2016 35 `````` `````` Jason Pries committed Dec 13, 2016 36 37 `````` // Initial residual calculation update_linearization(J, r); `````` Jason Pries committed Oct 11, 2016 38 `````` `````` Jason Pries committed Dec 13, 2016 39 40 41 42 43 44 45 `````` // Setup iteration 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 QR = J.fullPivHouseholderQr(); Eigen::VectorXd delta = QR.solve(r); `````` Jason Pries committed Oct 11, 2016 46 `````` `````` Jason Pries committed Dec 13, 2016 47 48 `````` 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; `````` Jason Pries committed Oct 11, 2016 49 50 `````` } `````` Jason Pries committed Dec 13, 2016 51 52 53 `````` long iter_rank = QR.rank(); if (iter_rank < prev_rank) { // perturb(delta, characteristic_length()); // perturb update in case of stagnation at inflection point `````` Jason Pries committed Oct 11, 2016 54 `````` } `````` Jason Pries committed Dec 13, 2016 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 `````` prev_rank = iter_rank; update_variables(delta); update_linearization(J, r); double iter_norm{r.norm()}; double iter_tol{res_tol}; if (iter_norm > res_norm) { // Line search delta *= -0.5; 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); iter_norm = r.norm(); iter_tol = dbl_tol * characteristic_length(); if (iter_norm > res_norm) { delta *= 0.5; line_search_iter += 1; 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; } } } `````` Jason Pries committed Oct 11, 2016 81 `````` } `````` Jason Pries committed Dec 13, 2016 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 `````` res_norm = iter_norm; res_tol = iter_tol; } return res_norm; } void Sketch::update_linearization(Eigen::MatrixXd &J, Eigen::VectorXd &r) const { J.setZero(); 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); } } `````` Jason Pries committed Oct 11, 2016 104 `````` `````` Jason Pries committed Dec 13, 2016 105 106 107 108 109 ``````void Sketch::update_variables(Eigen::VectorXd &delta) const { for (size_t i = 0; i != Variables.size(); ++i) { Variables[i]->update(delta); } } `````` Jason Pries committed Oct 11, 2016 110 `````` `````` Jason Pries committed Dec 13, 2016 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 ``````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 rdist = std::uniform_real_distribution(-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; } } double Sketch::characteristic_length() const { 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 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); `````` Jason Pries committed Oct 11, 2016 145 146 `````` } } `````` Jason Pries committed Dec 13, 2016 147 148 `````` return scale; `````` Jason Pries committed Oct 11, 2016 149 150 151 152 153 ``````} bool Sketch::build() { Constellation c = Constellation(this); `````` Jason Pries committed Dec 13, 2016 154 155 `````` Boundary = c.boundary(); bool success = (Boundary->size() > 0); `````` Jason Pries committed Oct 11, 2016 156 `````` `````` Jason Pries committed Dec 13, 2016 157 158 `````` Contours = c.contours(); success = success && (Contours.size() > 0); `````` Jason Pries committed Oct 11, 2016 159 160 161 162 `````` return success; } `````` Jason Pries committed Dec 13, 2016 163 ``````void Sketch::add_element(std::shared_ptr v) { `````` Jason Pries committed Oct 11, 2016 164 165 166 `````` add_element(Verticies, v); } `````` Jason Pries committed Dec 13, 2016 167 ``````void Sketch::add_element(std::shared_ptr c) { `````` Jason Pries committed Oct 11, 2016 168 169 170 `````` add_element(Curves, c); } `````` Jason Pries committed Dec 13, 2016 171 ``````void Sketch::add_element(std::shared_ptr c) { `````` Jason Pries committed Oct 11, 2016 172 173 174 `````` add_element(Constraints, c); } `````` Jason Pries committed Dec 13, 2016 175 ``````void Sketch::add_element(std::shared_ptr p) { `````` Jason Pries committed Oct 11, 2016 176 `````` add_element(Patterns, p); `````` Jason Pries committed Dec 13, 2016 177 `````` p->register_elements(this); `````` Jason Pries committed Oct 11, 2016 178 179 180 ``````} template<> `````` Jason Pries committed Dec 13, 2016 181 ``````void Sketch::save_as(std::string path, std::string file_name) const { `````` Jason Pries committed Oct 11, 2016 182 183 184 `````` /* This is a stub for visualization */ `````` Jason Pries committed Dec 13, 2016 185 186 187 188 189 `````` if (!boost::filesystem::exists(path)) { boost::filesystem::create_directories(path); } `````` Jason Pries committed Oct 11, 2016 190 `````` std::fstream fs; `````` Jason Pries committed Dec 13, 2016 191 192 `````` fs.open(path + file_name + ".oesk", std::fstream::out); `````` Jason Pries committed Oct 11, 2016 193 `````` for (size_t i = 0; i < Curves.size(); ++i) { `````` Jason Pries committed Dec 13, 2016 194 `````` if (!Curves[i]->for_construction()) { `````` Jason Pries committed Oct 11, 2016 195 `````` for (size_t j = 0; j <= 10; ++j) { `````` Jason Pries committed Dec 13, 2016 196 197 `````` double2 v = Curves[i]->point(j / 10.0); fs << v.X << ',' << v.Y << '\n'; `````` Jason Pries committed Oct 11, 2016 198 `````` } `````` Jason Pries committed Dec 13, 2016 199 `````` fs << "NaN" << ',' << "NaN" << '\n'; `````` Jason Pries committed Oct 11, 2016 200 201 `````` } } `````` Jason Pries committed Dec 13, 2016 202 `````` `````` Jason Pries committed Oct 11, 2016 203 `````` fs.close(); `````` JasonPries committed Feb 27, 2017 204 205 ``````} `````` JasonPries committed Feb 28, 2017 206 207 ``````std::vector Sketch::select_periodic_boundary_pairs(std::shared_ptr const &v0, double_t angle) const { std::vector cbp; `````` JasonPries committed Feb 27, 2017 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 `````` std::vector> bc = Boundary->curves(); while (bc.begin() != bc.end()) { auto c0 = bc.begin(); auto c1 = --bc.end(); MatchOrientation m{MatchOrientation::None}; while (c0 != c1) { auto cc0 = (*c0).get(); auto cc1 = (*c1).get(); m = (*c1)->is_identical(*c0, v0, angle); if (m == MatchOrientation::None) { m = (*c0)->is_identical(*c1, v0, angle); if (m != MatchOrientation::None) { std::swap(c0, c1); } } if (m == MatchOrientation::Forward) { `````` JasonPries committed Feb 28, 2017 229 `````` cbp.emplace_back(*c0, *c1, true, v0, angle); `````` JasonPries committed Feb 27, 2017 230 231 `````` break; } else if (m == MatchOrientation::Reverse) { `````` JasonPries committed Feb 28, 2017 232 `````` cbp.emplace_back(*c0, *c1, false, v0, angle); `````` JasonPries committed Feb 27, 2017 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 `````` break; } else { ++c0; } } if (m != MatchOrientation::None) { // if match found, erase c1 bc.erase(c0); bc.erase(c1); } else { bc.erase(c1); } } return cbp; } `````` JasonPries committed Feb 28, 2017 250 ``````std::vector> Sketch::select_radial_boundary(std::shared_ptr const &v0, double_t radius) const { `````` JasonPries committed Feb 27, 2017 251 252 253 254 255 256 257 258 259 260 261 262 `````` std::vector> rb; std::vector> bc = Boundary->curves(); auto test_circle = std::make_shared(std::make_shared(), std::make_shared(), v0, radius); for (std::shared_ptr c : bc) { if (test_circle->is_coincident(c)) { rb.push_back(c); } } return rb; `````` JasonPries committed Feb 28, 2017 263 264 265 266 267 268 269 270 271 272 273 ``````} std::shared_ptr Sketch::select_periodic_vertex(std::shared_ptr const &v, std::shared_ptr const &origin, double_t angle) const { for(auto const &v_test : Verticies) { if (v_test->is_identical(v, origin, angle)) { return v_test; } } std::shared_ptr nullv; return nullv; `````` Jason Pries committed Oct 11, 2016 274 ``}``