Skip to content

Commit b1bf722

Browse files
rjc8237rjc8237
committed
Abstracted flip jacobian constructor
1 parent 60f1ff8 commit b1bf722

File tree

2 files changed

+129
-95
lines changed

2 files changed

+129
-95
lines changed

‎src/optimization/cone_metric.cpp‎

Lines changed: 35 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -119,19 +119,13 @@ int DifferentiableConeMetric::n_reduced_coordinates() const
119119
}
120120

121121
PennerConeMetric::PennerConeMetric(const Mesh<Scalar>& m, const VectorX& metric_coords)
122-
: DifferentiableConeMetric(m)
122+
: DifferentiableConeMetric(m),
123+
m_need_jacobian(true),
124+
m_transition_jacobian_lol(m.n_edges())
123125
{
124126
build_refl_proj(m, he2e, e2he, m_proj, m_embed);
125127
build_refl_matrix(m_proj, m_embed, m_projection);
126128
expand_metric_coordinates(metric_coords);
127-
128-
// Initialize jacobian to the identity
129-
int num_edges = e2he.size();
130-
m_transition_jacobian_lol =
131-
std::vector<std::vector<std::pair<int, Scalar>>>(num_edges, std::vector<std::pair<int, Scalar>>());
132-
for (int e = 0; e < num_edges; ++e) {
133-
m_transition_jacobian_lol[e].push_back(std::make_pair(e, 1.0));
134-
}
135129
}
136130

137131
VectorX PennerConeMetric::reduce_metric_coordinates(const VectorX& metric_coords) const
@@ -177,21 +171,7 @@ std::unique_ptr<DifferentiableConeMetric> PennerConeMetric::scale_conformally(
177171

178172
MatrixX PennerConeMetric::get_transition_jacobian() const
179173
{
180-
std::vector<T> tripletList;
181-
int num_edges = e2he.size();
182-
tripletList.reserve(5 * num_edges);
183-
for (int i = 0; i < num_edges; ++i) {
184-
for (auto it : m_transition_jacobian_lol[i]) {
185-
tripletList.push_back(T(i, it.first, it.second));
186-
}
187-
}
188-
189-
// Create the matrix from the triplets
190-
MatrixX transition_jacobian;
191-
transition_jacobian.resize(num_edges, num_edges);
192-
transition_jacobian.reserve(tripletList.size());
193-
transition_jacobian.setFromTriplets(tripletList.begin(), tripletList.end());
194-
174+
MatrixX transition_jacobian = get_flip_jacobian();
195175
return m_identification * (transition_jacobian * m_projection);
196176
}
197177

@@ -224,55 +204,39 @@ bool PennerConeMetric::flip_ccw(int _h, bool Ptolemy)
224204
// Perform the flip in the base class
225205
bool success = DifferentiableConeMetric::flip_ccw(_h, Ptolemy);
226206

227-
// Get local mesh information near hd
228-
int hd = _h;
229-
int hb = n[hd];
230-
int ha = n[hb];
231-
int hdo = opp[hd];
232-
int hao = n[hdo];
233-
int hbo = n[hao];
234-
235-
// Get edges corresponding to halfedges
236-
int ed = he2e[hd];
237-
int eb = he2e[hb];
238-
int ea = he2e[ha];
239-
int eao = he2e[hao];
240-
int ebo = he2e[hbo];
241-
242-
// Compute the shear for the edge ed
243-
// TODO Check if edge or halfedge coordinate
244-
Scalar la = l[ha];
245-
Scalar lb = l[hb];
246-
Scalar lao = l[hao];
247-
Scalar lbo = l[hbo];
248-
Scalar x = (la * lbo) / (lb * lao);
249-
250-
// The matrix Pd corresponding to flipping edge ed is the identity except for
251-
// the row corresponding to edge ed, which has entries defined by Pd_scalars
252-
// in the column corresponding to the edge with the same index in Pd_edges
253-
std::vector<int> Pd_edges = {ed, ea, ebo, eao, eb};
254-
std::vector<Scalar> Pd_scalars =
255-
{-1.0, x / (1.0 + x), x / (1.0 + x), 1.0 / (1.0 + x), 1.0 / (1.0 + x)};
256-
int num_entries = 0;
257-
for (int i = 0; i < 5; ++i) {
258-
int ei = Pd_edges[i];
259-
num_entries += m_transition_jacobian_lol[ei].size();
207+
if (m_need_jacobian) {
208+
// Get local mesh information near hd
209+
int hd = _h;
210+
int hb = n[hd];
211+
int ha = n[hb];
212+
int hdo = opp[hd];
213+
int hao = n[hdo];
214+
int hbo = n[hao];
215+
216+
// Get edges corresponding to halfedges
217+
int ed = he2e[hd];
218+
int eb = he2e[hb];
219+
int ea = he2e[ha];
220+
int eao = he2e[hao];
221+
int ebo = he2e[hbo];
222+
223+
// Compute the shear for the edge ed
224+
// TODO Check if edge or halfedge coordinate
225+
Scalar la = l[ha];
226+
Scalar lb = l[hb];
227+
Scalar lao = l[hao];
228+
Scalar lbo = l[hbo];
229+
Scalar x = (la * lbo) / (lb * lao);
230+
231+
// The matrix Pd corresponding to flipping edge ed is the identity except for
232+
// the row corresponding to edge ed, which has entries defined by Pd_scalars
233+
// in the column corresponding to the edge with the same index in Pd_edges
234+
std::vector<int> Pd_edges = {ed, ea, ebo, eao, eb};
235+
std::vector<Scalar> Pd_scalars =
236+
{-1.0, x / (1.0 + x), x / (1.0 + x), 1.0 / (1.0 + x), 1.0 / (1.0 + x)};
237+
m_transition_jacobian_lol.multiply_by_matrix(Pd_edges, Pd_scalars, ed);
260238
}
261239

262-
// Compute the new row of J_del corresponding to edge ed, which is the only
263-
// edge that changes
264-
std::vector<std::pair<int, Scalar>> J_del_d_new;
265-
J_del_d_new.reserve(num_entries);
266-
for (int i = 0; i < 5; ++i) {
267-
int ei = Pd_edges[i];
268-
Scalar Di = Pd_scalars[i];
269-
for (auto it : m_transition_jacobian_lol[ei]) {
270-
J_del_d_new.push_back(std::make_pair(it.first, Di * it.second));
271-
}
272-
}
273-
274-
m_transition_jacobian_lol[ed] = J_del_d_new;
275-
276240
return success;
277241
}
278242

‎src/optimization/cone_metric.hh‎

Lines changed: 94 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,96 @@
77

88
namespace CurvatureMetric {
99

10+
class FlipMatrixGenerator
11+
{
12+
public:
13+
FlipMatrixGenerator(int size)
14+
: m_size(size) {
15+
reset();
16+
17+
}
18+
19+
void reset() {
20+
m_list_of_lists =
21+
std::vector<std::vector<std::pair<int, Scalar>>>(m_size, std::vector<std::pair<int, Scalar>>());
22+
for (int i = 0; i < m_size; ++i) {
23+
m_list_of_lists[i].push_back(std::make_pair(i, 1.0));
24+
}
25+
}
26+
27+
void multiply_by_matrix(
28+
const std::vector<int>& matrix_indices,
29+
const std::vector<Scalar>& matrix_scalars,
30+
int ed
31+
) {
32+
int num_entries = 0;
33+
for (int i = 0; i < 5; ++i) {
34+
int ei = matrix_indices[i];
35+
num_entries += m_list_of_lists[ei].size();
36+
}
37+
38+
// Compute the new row of J_del corresponding to edge ed, which is the only
39+
// edge that changes
40+
std::vector<std::pair<int, Scalar>> J_del_d_new;
41+
J_del_d_new.reserve(num_entries);
42+
for (int i = 0; i < 5; ++i) {
43+
int ei = matrix_indices[i];
44+
Scalar Di = matrix_scalars[i];
45+
for (auto it : m_list_of_lists[ei]) {
46+
J_del_d_new.push_back(std::make_pair(it.first, Di * it.second));
47+
}
48+
}
49+
std::sort(J_del_d_new.begin(), J_del_d_new.end());
50+
51+
// Compress vector
52+
m_list_of_lists[ed] = std::vector<std::pair<int, Scalar>>();
53+
m_list_of_lists[ed].reserve(num_entries);
54+
int index = -1;
55+
Scalar value = 0.0;
56+
for (const auto& entry : J_del_d_new)
57+
{
58+
if (index != entry.first)
59+
{
60+
if (index >= 0)
61+
{
62+
m_list_of_lists[ed].push_back(std::make_pair(index, value));
63+
}
64+
index = entry.first;
65+
value = entry.second;
66+
}
67+
else
68+
{
69+
value += entry.second;
70+
}
71+
}
72+
m_list_of_lists[ed].push_back(std::make_pair(index, value));
73+
}
74+
75+
MatrixX build_matrix() const {
76+
// Build triplets from list of lists
77+
typedef Eigen::Triplet<Scalar> T;
78+
std::vector<T> tripletList;
79+
tripletList.reserve(5 * m_size);
80+
for (int i = 0; i < m_size; ++i) {
81+
for (auto it : m_list_of_lists[i]) {
82+
tripletList.push_back(T(i, it.first, it.second));
83+
}
84+
}
85+
86+
// Create the matrix from the triplets
87+
MatrixX matrix;
88+
matrix.resize(m_size, m_size);
89+
matrix.reserve(tripletList.size());
90+
matrix.setFromTriplets(tripletList.begin(), tripletList.end());
91+
92+
return matrix;
93+
}
94+
95+
private:
96+
int m_size;
97+
std::vector<std::vector<std::pair<int, Scalar>>> m_list_of_lists;
98+
};
99+
10100
class DifferentiableConeMetric : public Mesh<Scalar>
11101
{
12102
public:
@@ -95,42 +185,22 @@ public:
95185

96186
MatrixX get_flip_jacobian() const
97187
{
98-
typedef Eigen::Triplet<Scalar> T;
99-
std::vector<T> tripletList;
100-
int num_edges = e2he.size();
101-
tripletList.reserve(5 * num_edges);
102-
for (int i = 0; i < num_edges; ++i) {
103-
for (auto it : m_transition_jacobian_lol[i]) {
104-
tripletList.push_back(T(i, it.first, it.second));
105-
}
106-
}
107-
108-
// Create the matrix from the triplets
109-
MatrixX transition_jacobian;
110-
transition_jacobian.resize(num_edges, num_edges);
111-
transition_jacobian.reserve(tripletList.size());
112-
transition_jacobian.setFromTriplets(tripletList.begin(), tripletList.end());
113-
114-
return transition_jacobian;
188+
return m_transition_jacobian_lol.build_matrix();
115189
}
116190

117191
// Reset Jacobian
118192
void reset()
119193
{
120194
// Initialize jacobian to the identity
121-
int num_edges = e2he.size();
122-
m_transition_jacobian_lol =
123-
std::vector<std::vector<std::pair<int, Scalar>>>(num_edges, std::vector<std::pair<int, Scalar>>());
124-
for (int e = 0; e < num_edges; ++e) {
125-
m_transition_jacobian_lol[e].push_back(std::make_pair(e, 1.0));
126-
}
195+
m_transition_jacobian_lol.reset();
127196
}
128197

129198
protected:
130199
std::vector<int> m_embed;
131200
std::vector<int> m_proj;
132201
MatrixX m_projection;
133-
std::vector<std::vector<std::pair<int, Scalar>>> m_transition_jacobian_lol;
202+
bool m_need_jacobian;
203+
FlipMatrixGenerator m_transition_jacobian_lol;
134204

135205
VectorX reduce_metric_coordinates(const VectorX& metric_coords) const;
136206
void expand_metric_coordinates(const VectorX& metric_coords);

0 commit comments

Comments
 (0)