Sketch.cpp 6.57 KB
 Jason Pries committed Dec 13, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 ``````#include #include "Sketch.h" #include "Vertex.h" #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 26 `````` `````` Jason Pries committed Dec 13, 2016 27 28 29 `````` // 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 30 `````` `````` Jason Pries committed Dec 13, 2016 31 `````` size_t double_bits = (CHAR_BIT * sizeof(double)); `````` Jason Pries committed Oct 11, 2016 32 `````` `````` Jason Pries committed Dec 13, 2016 33 34 `````` // Initial residual calculation update_linearization(J, r); `````` Jason Pries committed Oct 11, 2016 35 `````` `````` Jason Pries committed Dec 13, 2016 36 37 38 39 40 41 42 `````` // 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 43 `````` `````` Jason Pries committed Dec 13, 2016 44 45 `````` 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 46 47 `````` } `````` Jason Pries committed Dec 13, 2016 48 49 50 `````` 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 51 `````` } `````` Jason Pries committed Dec 13, 2016 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 `````` 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 78 `````` } `````` Jason Pries committed Dec 13, 2016 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 `````` 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 101 `````` `````` Jason Pries committed Dec 13, 2016 102 103 104 105 106 ``````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 107 `````` `````` Jason Pries committed Dec 13, 2016 108 109 110 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 ``````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 142 143 `````` } } `````` Jason Pries committed Dec 13, 2016 144 145 `````` return scale; `````` Jason Pries committed Oct 11, 2016 146 147 148 149 150 ``````} bool Sketch::build() { Constellation c = Constellation(this); `````` Jason Pries committed Dec 13, 2016 151 152 `````` Boundary = c.boundary(); bool success = (Boundary->size() > 0); `````` Jason Pries committed Oct 11, 2016 153 `````` `````` Jason Pries committed Dec 13, 2016 154 155 `````` Contours = c.contours(); success = success && (Contours.size() > 0); `````` Jason Pries committed Oct 11, 2016 156 157 158 159 `````` return success; } `````` Jason Pries committed Dec 13, 2016 160 ``````void Sketch::add_element(std::shared_ptr v) { `````` Jason Pries committed Oct 11, 2016 161 162 163 `````` add_element(Verticies, v); } `````` Jason Pries committed Dec 13, 2016 164 ``````void Sketch::add_element(std::shared_ptr c) { `````` Jason Pries committed Oct 11, 2016 165 166 167 `````` add_element(Curves, c); } `````` Jason Pries committed Dec 13, 2016 168 ``````void Sketch::add_element(std::shared_ptr c) { `````` Jason Pries committed Oct 11, 2016 169 170 171 `````` add_element(Constraints, c); } `````` Jason Pries committed Dec 13, 2016 172 ``````void Sketch::add_element(std::shared_ptr p) { `````` Jason Pries committed Oct 11, 2016 173 `````` add_element(Patterns, p); `````` Jason Pries committed Dec 13, 2016 174 `````` p->register_elements(this); `````` Jason Pries committed Oct 11, 2016 175 176 177 ``````} template<> `````` Jason Pries committed Dec 13, 2016 178 ``````void Sketch::save_as(std::string path, std::string file_name) const { `````` Jason Pries committed Oct 11, 2016 179 180 181 `````` /* This is a stub for visualization */ `````` Jason Pries committed Dec 13, 2016 182 183 184 185 186 `````` if (!boost::filesystem::exists(path)) { boost::filesystem::create_directories(path); } `````` Jason Pries committed Oct 11, 2016 187 `````` std::fstream fs; `````` Jason Pries committed Dec 13, 2016 188 189 `````` fs.open(path + file_name + ".oesk", std::fstream::out); `````` Jason Pries committed Oct 11, 2016 190 `````` for (size_t i = 0; i < Curves.size(); ++i) { `````` Jason Pries committed Dec 13, 2016 191 `````` if (!Curves[i]->for_construction()) { `````` Jason Pries committed Oct 11, 2016 192 `````` for (size_t j = 0; j <= 10; ++j) { `````` Jason Pries committed Dec 13, 2016 193 194 `````` double2 v = Curves[i]->point(j / 10.0); fs << v.X << ',' << v.Y << '\n'; `````` Jason Pries committed Oct 11, 2016 195 `````` } `````` Jason Pries committed Dec 13, 2016 196 `````` fs << "NaN" << ',' << "NaN" << '\n'; `````` Jason Pries committed Oct 11, 2016 197 198 `````` } } `````` Jason Pries committed Dec 13, 2016 199 `````` `````` Jason Pries committed Oct 11, 2016 200 201 `````` fs.close(); }``````