1 /* Copyright (c) 1997-2021
2 Ewgenij Gawrilow, Michael Joswig, and the polymake team
3 Technische Universität Berlin, Germany
4 https://polymake.org
5
6 This program is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version: http://www.gnu.org/licenses/gpl.txt.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 --------------------------------------------------------------------------------
16 */
17
18 #include "polymake/client.h"
19 #include "polymake/Set.h"
20 #include "polymake/Graph.h"
21 #include "polymake/Matrix.h"
22 #include "polymake/linalg.h"
23 #include "polymake/Rational.h"
24 #include "polymake/Array.h"
25 #include "polymake/list"
26 #include "polymake/vector"
27 #include <cmath>
28
29 namespace polymake { namespace graph {
30 namespace {
31
32 template <typename Coords> inline
max_norm(const Matrix<Coords> & V,Int i,Int j)33 Coords max_norm(const Matrix<Coords>& V, Int i, Int j)
34 {
35 return accumulate(attach_operation(V[i]-V[j], operations::abs_value()), operations::max());
36 }
37
38 template <typename Coords> inline
square_norm(const Matrix<Coords> & V,Int i,Int j)39 Coords square_norm(const Matrix<Coords>& V, Int i, Int j)
40 {
41 return static_cast<Coords>(std::sqrt(double(sqr(V[i] - V[j]))));
42 }
43 }
44
45
46 template <typename Coords>
tentacle_graph(const Array<Int> & tentacles,const Matrix<Coords> & metric)47 ListReturn tentacle_graph(const Array<Int>& tentacles, const Matrix<Coords>& metric)
48 {
49 const Int n_nodes = tentacles.size();
50 Graph<> G(n_nodes);
51 EdgeMap<Undirected,Coords> weights(G);
52 for (Int i = 1; i < n_nodes; ++i)
53 for (Int j = 0; j < i; ++j)
54 weights(i, j) = metric(tentacles[i], tentacles[j]);
55
56 ListReturn results;
57 results << G << weights;
58 return results;
59 }
60
61
62 template <typename Coords>
bounded_embedder(const Graph<> & BG,const Matrix<Coords> & V,const Set<Int> & far_face,const Array<Int> & fixed_nodes,const Matrix<Coords> & fixed_coord,bool use_max_norm)63 Matrix<Coords> bounded_embedder(const Graph<>& BG, const Matrix<Coords>& V, const Set<Int>& far_face,
64 const Array<Int>& fixed_nodes, const Matrix<Coords>& fixed_coord,
65 bool use_max_norm)
66 {
67 const Int n_nodes = BG.nodes() + far_face.size();
68 Matrix<Coords> GR(n_nodes,3);
69 const Int n_fixed_nodes = fixed_nodes.size();
70 if (n_fixed_nodes < 3)
71 throw std::runtime_error("bounded_embedder: Less than three fixed nodes.");
72
73 GR.minor(fixed_nodes,All) = fixed_coord;
74
75 #if POLYMAKE_DEBUG
76 const bool debug_print = get_debug_level() > 1;
77 #endif
78
79 // compute distances according to metric chosen
80
81 Coords scale=zero_value<Coords>();
82 if (use_max_norm) {
83 for (auto fix1 = entire(fixed_nodes); !fix1.at_end(); ++fix1) {
84 auto fix2 = fix1;
85 while (!(++fix2).at_end())
86 scale += max_norm(V,*fix1,*fix2)/square_norm(GR,*fix1,*fix2);
87 }
88 }
89 else {
90 for (auto fix1 = entire(fixed_nodes); !fix1.at_end(); ++fix1) {
91 auto fix2 = fix1;
92 while (!(++fix2).at_end())
93 scale += square_norm(V,*fix1,*fix2)/square_norm(GR,*fix1,*fix2);
94 }
95 }
96 scale /= double(n_fixed_nodes)*(double(n_fixed_nodes)-1)/2;
97
98 EdgeMap<Undirected, Coords> BGmap(BG);
99 for (auto e=entire(edges(BG)); !e.at_end(); ++e)
100 BGmap[*e] = use_max_norm ? max_norm(V, e.from_node(), e.to_node())/scale : square_norm(V, e.from_node(), e.to_node())/scale;
101
102 #if POLYMAKE_DEBUG
103 if (debug_print) {
104 cout << "Vertices:\n" << V << "\n"
105 "Graph:\n" << BG << "\n"
106 "EdgeMap:\n " << BGmap << " \n "
107 "fixed_nodes:\n" << GR.minor(fixed_nodes,All) << "\n"
108 "scale: " << scale << endl;
109 }
110 #endif
111
112 // compute coordinates of remaining nodes
113 const Set<Int> inner_nodes = sequence(0, n_nodes) - Set<Int>(entire(fixed_nodes)) - far_face;
114 const Int n_inner_nodes = inner_nodes.size();
115 Matrix<Coords> stress_matrix(n_inner_nodes, n_nodes);
116 Matrix<Coords> rhs(n_inner_nodes,3);
117
118 Int m = 0;
119 for (auto n = entire(inner_nodes); !n.at_end(); ++n,++m) {
120 for (auto e = entire(BG.out_edges(*n)); !e.at_end(); ++e) {
121 stress_matrix(m,*n) += 1/BGmap[*e]; // spring constant
122 const Int nn = e.to_node();
123 if ( inner_nodes.contains(nn) )
124 stress_matrix(m,nn) = -1/BGmap[*e]; // spring constant
125 else
126 rhs[m] += GR[nn]/BGmap[*e];
127 }
128 }
129
130 GR.minor(inner_nodes,All) = inv(stress_matrix.minor(All,inner_nodes)) * rhs;
131 #if POLYMAKE_DEBUG
132 if (debug_print) {
133 cout << "inner_nodes: " << inner_nodes << "\n"
134 "stress_matrix:\n" << stress_matrix << "\n"
135 "rhs:\n" << rhs << "\n"
136 "stress_matrix_inv:\n" << inv(stress_matrix.minor(All,inner_nodes)) << "\n"
137 "result:\n" << GR << endl;
138 }
139 #endif
140
141 return GR.minor(~far_face,All);
142 }
143
144 FunctionTemplate4perl("bounded_embedder($ Matrix $$ Matrix; $=1)");
145 FunctionTemplate4perl("tentacle_graph($ Matrix)");
146
147 } }
148
149 // Local Variables:
150 // mode:C++
151 // c-basic-offset:3
152 // indent-tabs-mode:nil
153 // End:
154