Sketch.cpp 6.57 KB
Newer Older
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 <random>

#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 <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
26

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

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

33
34
    // Initial residual calculation
    update_linearization(J, r);
Jason Pries's avatar
Jason Pries committed
35

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

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

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's avatar
Jason Pries committed
51
        }
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's avatar
Jason Pries committed
78
        }
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's avatar
Jason Pries committed
101

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

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<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
142
143
        }
    }
144
145

    return scale;
Jason Pries's avatar
Jason Pries committed
146
147
148
149
150
}

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

151
152
    Boundary = c.boundary();
    bool success = (Boundary->size() > 0);
Jason Pries's avatar
Jason Pries committed
153

154
155
    Contours = c.contours();
    success = success && (Contours.size() > 0);
Jason Pries's avatar
Jason Pries committed
156
157
158
159

    return success;
}

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

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

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

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

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

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

Jason Pries's avatar
Jason Pries committed
187
    std::fstream fs;
188
189
    fs.open(path + file_name + ".oesk", std::fstream::out);

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

Jason Pries's avatar
Jason Pries committed
200
201
    fs.close();
}