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/Matrix.h"
20 #include "polymake/Vector.h"
21 #include "polymake/PowerSet.h"
22 #include "polymake/Array.h"
23 #include "polymake/Rational.h"
24 #include "polymake/Map.h"
25 #include "polymake/Graph.h"
26 #include "polymake/list"
27 #include <string>
28
29 namespace polymake { namespace topaz {
30 namespace {
31
32 // m_set is an m-element set of {0,...,m+n-1}
shuffle_from_m_set(const Int m,const Int n,const Set<Int> & m_set)33 std::vector<bool> shuffle_from_m_set(const Int m, const Int n, const Set<Int>& m_set)
34 {
35 std::vector<bool> shuffle(m+n-2, false);
36 for (const Int i : m_set)
37 shuffle[i] = true;
38 return shuffle;
39 }
40
41 // m_set is an m-element set of {0,...,m+n-1}
facet_from_m_set(const std::list<Int> & s1,const std::list<Int> & s2,const Matrix<Int> & vert_map,const Set<Int> & m_set)42 Set<Int> facet_from_m_set(const std::list<Int>& s1, const std::list<Int>& s2,
43 const Matrix<Int>& vert_map, const Set<Int>& m_set)
44 {
45 auto v1 = s1.begin();
46 auto v2 = s2.begin();
47 Set<Int> facet;
48 facet += vert_map(*v1, *v2);
49 // cout << *v1 << "," << *v2 << " -- " << vert_map(*v1,*v2) << endl;
50 const std::vector<bool> shuffle = shuffle_from_m_set(s1.size(), s2.size(), m_set);
51 for (bool s : shuffle) {
52 if (s) ++v1;
53 else ++v2;
54 facet += vert_map(*v1, *v2);
55 }
56
57 return facet;
58 }
59
60
color_cons_ordering(const Array<Int> & C)61 Array<Int> color_cons_ordering(const Array<Int>& C)
62 {
63 Array<Int> ordering(C.size());
64
65 Set<Int> colors;
66 Map<Int, std::list<Int>> vert_of_color;
67 for (Int i = 0; i < C.size(); ++i) {
68 colors += C[i];
69 vert_of_color[ C[i] ].push_back(i);
70 }
71
72 Int i = 0;
73 for (auto c = entire(colors); !c.at_end(); ++c)
74 for (auto v = entire(vert_of_color[*c]); !v.at_end(); ++v, ++i)
75 ordering[i] = *v;
76
77 return ordering;
78 }
79
combinatorial_simplicial_product_impl(BigObject p_in1,BigObject p_in2,BigObject & p_out,Array<Int> & order1,Array<Int> & order2,OptionSet options)80 void combinatorial_simplicial_product_impl(BigObject p_in1, BigObject p_in2, BigObject& p_out, Array<Int>& order1, Array<Int>& order2, OptionSet options)
81 {
82 const bool no_labels = options["no_labels"];
83 const Array<Set<Int>> C1 = p_in1.give("FACETS");
84 const Array<Set<Int>> C2 = p_in2.give("FACETS");
85 const Int n_vert1 = p_in1.give("N_VERTICES");
86 const Int n_vert2 = p_in2.give("N_VERTICES");
87
88 // read orderings or set default
89 order1 = Array<Int>(n_vert1);
90 order2 = Array<Int>(n_vert2);
91
92 if (!(options["vertex_order1"] >> order1)) {
93 if (options["color_cons"]) {
94 const bool is_foldable = p_in1.give("FOLDABLE");
95 if (!is_foldable)
96 throw std::runtime_error("simplicial_product: Complex 1 is not FOLDABLE");
97 const Array<Int> coloring = p_in1.give("COLORING");
98 order1 = color_cons_ordering(coloring);
99 } else {
100 for (Int i = 0; i < n_vert1; ++i)
101 order1[i] = i;
102 }
103 }
104
105 if (!(options["vertex_order2"] >> order2)) {
106 if (options["color_cons"]) {
107 const bool is_foldable = p_in2.give("FOLDABLE");
108 if (!is_foldable)
109 throw std::runtime_error("simplicial_product: Complex 2 is not FOLDABLE");
110 const Array<Int> coloring = p_in2.give("COLORING");
111 order2 = color_cons_ordering(coloring);
112 } else {
113 for (Int i = 0; i < n_vert2; ++i)
114 order2[i] = i;
115 }
116 }
117
118 // compute vertex map
119 Matrix<Int> vert_map(n_vert1, n_vert2);
120 Int c = 0;
121 for (auto v2 = entire(order2); !v2.at_end(); ++v2)
122 for (auto v1=entire(order1); !v1.at_end(); ++v1, ++c)
123 vert_map(*v1,*v2) = c;
124
125 std::list<Set<Int>> F;
126 for (auto f1 = entire(C1); !f1.at_end(); ++f1)
127 for (auto f2 = entire(C2); !f2.at_end(); ++f2) {
128 const Int m = f1->size()-1;
129 const Int n = f2->size()-1;
130
131 std::list<Int> s1, s2;
132 for (Int i = 0; i < n_vert1; ++i)
133 if (f1->contains(order1[i]))
134 s1.push_back(order1[i]);
135 for (Int i = 0; i < n_vert2; ++i)
136 if (f2->contains(order2[i]))
137 s2.push_back(order2[i]);
138
139 // compute faces
140 for (auto m_set=entire(all_subsets_of_k(range(0,m+n-1), m)); !m_set.at_end(); ++m_set)
141 F.push_back(facet_from_m_set(s1, s2, vert_map, *m_set));
142 }
143
144 p_out.set_description() << "Simplicial product of " << p_in1.name() << " and " << p_in2.name() << "."<<endl;
145 p_out.take("FACETS") << as_array(F);
146
147 if (!no_labels) {
148 const Array<std::string> L1 = p_in1.give("VERTEX_LABELS");
149 const Array<std::string> L2 = p_in2.give("VERTEX_LABELS");
150 Array<std::string> L(L1.size()*L2.size());
151
152 c=0;
153 for (auto v2=entire(order2); !v2.at_end(); ++v2)
154 for (auto v1=entire(order1); !v1.at_end(); ++v1, ++c)
155 L[c] = '(' + L1[*v1] + ',' + L2[*v2] + ')';
156
157 p_out.take("VERTEX_LABELS") << L;
158 }
159 }
160
combinatorial_simplicial_product(BigObject p_in1,BigObject p_in2,OptionSet options)161 BigObject combinatorial_simplicial_product (BigObject p_in1, BigObject p_in2, OptionSet options)
162 {
163 BigObject p_out("SimplicialComplex");
164 Array<Int> order1, order2;
165 combinatorial_simplicial_product_impl(p_in1, p_in2, p_out, order1, order2, options);
166 return p_out;
167 }
168
169 template <typename Scalar>
simplicial_product(BigObject p_in1,BigObject p_in2,OptionSet options)170 BigObject simplicial_product(BigObject p_in1, BigObject p_in2, OptionSet options)
171 {
172 const bool realize = options["geometric_realization"];
173 BigObjectType result_type = realize
174 ? BigObjectType("GeometricSimplicialComplex", mlist<Scalar>())
175 : BigObjectType("SimplicialComplex");
176 BigObject p_out(result_type);
177 Array<Int> order1, order2;
178 combinatorial_simplicial_product_impl(p_in1, p_in2, p_out, order1, order2, options);
179
180 if (realize) {
181 const Matrix<Scalar> GR1 = p_in1.give("COORDINATES");
182 const Matrix<Scalar> GR2 = p_in2.give("COORDINATES");
183 Matrix<Scalar> GR(GR1.rows()*GR2.rows(),GR1.cols()+GR2.cols());
184 Int c = 0;
185 for (auto v2=entire(order2); !v2.at_end(); ++v2)
186 for (auto v1=entire(order1); !v1.at_end(); ++v1, ++c)
187 GR[c] = GR1[*v1] | GR2[*v2];
188
189 p_out.take("COORDINATES") << GR;
190 }
191
192 return p_out;
193 }
194
195 } // end anonymous namespace
196
197 UserFunction4perl("# @category Producing a new simplicial complex from others\n"
198 "# Computes the __simplicial product__ of two complexes.\n"
199 "# Vertex orderings may be given as options.\n"
200 "# @param SimplicialComplex complex1"
201 "# @param SimplicialComplex complex2"
202 "# @option Array<Int> vertex_order1"
203 "# @option Array<Int> vertex_order2"
204 "# @option Bool geometric_realization default 0"
205 "# @option Bool color_cons"
206 "# @option Bool no_labels Do not create [[VERTEX_LABELS]]. default: 0"
207 "# @return SimplicialComplex"
208 "# @example The following returns the product of two edges."
209 "# > $s = simplicial_product(simplex(1), simplex(1));"
210 "# > print $s -> F_VECTOR;"
211 "# | 4 5 2",
212 &combinatorial_simplicial_product, "simplicial_product(SimplicialComplex, SimplicialComplex, {vertex_order1 => undef, vertex_order2 => undef, geometric_realization => 0, color_cons => 0, no_labels => 0})");
213
214 UserFunctionTemplate4perl("# @category Producing a new simplicial complex from others\n"
215 "# Computes the __simplicial product__ of two complexes.\n"
216 "# Vertex orderings may be given as options.\n"
217 "# @param GeometricSimplicialComplex complex1"
218 "# @param GeometricSimplicialComplex complex2"
219 "# @tparam Scalar"
220 "# @option Array<Int> vertex_order1"
221 "# @option Array<Int> vertex_order2"
222 "# @option Bool geometric_realization default 1"
223 "# @option Bool color_cons"
224 "# @option Bool no_labels Do not create [[VERTEX_LABELS]]. default: 0"
225 "# @return GeometricSimplicialComplex<Scalar>"
226 "# @example The following returns the product of the edges (0, 0)--(1, 0) and (0, 0) -- (2, 0)."
227 "# > $C = new GeometricSimplicialComplex(COORDINATES => [[0, 0], [1, 0]], FACETS => [[0, 1]]);"
228 "# > $C1 = new GeometricSimplicialComplex(COORDINATES => [[0, 2], [0, 0]], FACETS => [[0, 1]]);"
229 "# > $s = simplicial_product($C, $C1);"
230 "# > print $s -> COORDINATES;"
231 "# | 0 0 0 2"
232 "# | 1 0 0 2"
233 "# | 0 0 0 0"
234 "# | 1 0 0 0",
235 "simplicial_product<Scalar>(GeometricSimplicialComplex<Scalar>, GeometricSimplicialComplex<Scalar>, {vertex_order1 => undef, vertex_order2 => undef, geometric_realization => 1, color_cons => 0, no_labels => 0})");
236
237
238 } }
239
240 // Local Variables:
241 // mode:C++
242 // c-basic-offset:3
243 // indent-tabs-mode:nil
244 // End:
245