Sketch.cpp 8.49 KB
Newer Older
1
#include <random>
2
#include <utility>
3 4 5

#include "Sketch.h"
#include "Vertex.h"
6 7
#include "CircularArc.h"
#include "LineSegment.h"
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 <iostream>

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's avatar
Jason Pries committed
29

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's avatar
Jason Pries committed
33

34
    size_t double_bits = (CHAR_BIT * sizeof(double));
Jason Pries's avatar
Jason Pries committed
35

36 37
    // Initial residual calculation
    update_linearization(J, r);
Jason Pries's avatar
Jason Pries committed
38

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<Eigen::MatrixXd> QR = J.fullPivHouseholderQr();
        Eigen::VectorXd delta = QR.solve(r);
Jason Pries's avatar
Jason Pries committed
46

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's avatar
Jason Pries committed
49 50
        }

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's avatar
Jason Pries committed
54
        }
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's avatar
Jason Pries committed
81
        }
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's avatar
Jason Pries committed
104

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's avatar
Jason Pries committed
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 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<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;
    }
}

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's avatar
Jason Pries committed
145 146
        }
    }
147 148

    return scale;
Jason Pries's avatar
Jason Pries committed
149 150 151 152 153
}

bool Sketch::build() {
    Constellation c = Constellation(this);

154 155
    Boundary = c.boundary();
    bool success = (Boundary->size() > 0);
Jason Pries's avatar
Jason Pries committed
156

157 158
    Contours = c.contours();
    success = success && (Contours.size() > 0);
Jason Pries's avatar
Jason Pries committed
159 160 161 162

    return success;
}

163
void Sketch::add_element(std::shared_ptr<Vertex> v) {
Jason Pries's avatar
Jason Pries committed
164 165 166
    add_element(Verticies, v);
}

167
void Sketch::add_element(std::shared_ptr<Curve> c) {
Jason Pries's avatar
Jason Pries committed
168 169 170
    add_element(Curves, c);
}

171
void Sketch::add_element(std::shared_ptr<Constraint> c) {
Jason Pries's avatar
Jason Pries committed
172 173 174
    add_element(Constraints, c);
}

175
void Sketch::add_element(std::shared_ptr<Pattern> p) {
Jason Pries's avatar
Jason Pries committed
176
    add_element(Patterns, p);
177
    p->register_elements(this);
Jason Pries's avatar
Jason Pries committed
178 179 180
}

template<>
181
void Sketch::save_as<SaveMethod::Rasterize>(std::string path, std::string file_name) const {
Jason Pries's avatar
Jason Pries committed
182 183 184
    /*
        This is a stub for visualization
    */
185 186 187 188 189

    if (!boost::filesystem::exists(path)) {
        boost::filesystem::create_directories(path);
    }

Jason Pries's avatar
Jason Pries committed
190
    std::fstream fs;
191 192
    fs.open(path + file_name + ".oesk", std::fstream::out);

Jason Pries's avatar
Jason Pries committed
193
    for (size_t i = 0; i < Curves.size(); ++i) {
194
        if (!Curves[i]->for_construction()) {
Jason Pries's avatar
Jason Pries committed
195
            for (size_t j = 0; j <= 10; ++j) {
196 197
                double2 v = Curves[i]->point(j / 10.0);
                fs << v.X << ',' << v.Y << '\n';
Jason Pries's avatar
Jason Pries committed
198
            }
199
            fs << "NaN" << ',' << "NaN" << '\n';
Jason Pries's avatar
Jason Pries committed
200 201
        }
    }
202

Jason Pries's avatar
Jason Pries committed
203
    fs.close();
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
}

std::vector<ContinuousBoundaryPair> Sketch::select_periodic_boundary_pairs(std::shared_ptr<Vertex> v0, double_t angle) const {
    std::vector<ContinuousBoundaryPair> cbp;

    std::vector<std::shared_ptr<Curve const>> 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) {
                cbp.emplace_back(*c0, *c1, true);
                break;
            } else if (m == MatchOrientation::Reverse) {
                cbp.emplace_back(*c0, *c1, false);
                break;
            } else {
                ++c0;
            }
        }

        if (m != MatchOrientation::None) { // if match found, erase c1
            bc.erase(c0);
            bc.erase(c1);
        } else {
            bc.erase(c1);
        }
    }

    return cbp;
}

std::vector<std::shared_ptr<Curve const>> Sketch::select_radial_boundary(std::shared_ptr<Vertex> v0, double_t radius) const {
    std::vector<std::shared_ptr<Curve const>> rb;

    std::vector<std::shared_ptr<Curve const>> bc = Boundary->curves();

    auto test_circle = std::make_shared<CircularArc>(std::make_shared<Vertex>(), std::make_shared<Vertex>(), v0, radius);
    for (std::shared_ptr<Curve const> c : bc) {
        if (test_circle->is_coincident(c)) {
            rb.push_back(c);
        }
    }

    return rb;
Jason Pries's avatar
Jason Pries committed
263
}