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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
|
/* Copyright (c) 1997-2024
Ewgenij Gawrilow, Michael Joswig, and the polymake team
Technische Universität Berlin, Germany
https://polymake.org
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version: http://www.gnu.org/licenses/gpl.txt.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
--------------------------------------------------------------------------------
*/
#pragma once
#include "polymake/client.h"
#include "polymake/Matrix.h"
#include "polymake/ReverseSearch.h"
#include "polymake/group/action.h"
#include "polymake/Set.h"
namespace polymake {
namespace fan {
namespace reverse_search_chamber_decomposition {
template<typename Scalar>
Vector<Scalar> signature_to_vertex(const Matrix<Scalar>& hyp, const Bitset& signature, const Vector<Scalar>& apv){
Vector<Scalar> result = ones_vector<Scalar>(hyp.rows());
result.slice(~signature) *= -1;
result = T(hyp) * result;
return -apv*result | result;
}
template<typename Scalar, typename CacheType>
class Node {
private:
const Matrix<Scalar>& hyperplanes;
Bitset signature;
CacheType& cache;
Vector<Scalar> vertex;
Map<Vector<Scalar>, Bitset> upNeighbors, downNeighbors;
Bitset neighbor_signature_from_facet(const Vector<Scalar>& facet, bool& facet_is_hyperplane){
// This should be done more efficiently:
// We could normalize the hyperplanes and check for equality. Although
// we would have to modify signatures in this case.
Bitset result(signature);
Int i = 0;
Matrix<Scalar> tmp(cache.get_support_eq());
tmp /= facet;
Int tmprk = rank(tmp);
for(const auto& f : rows(hyperplanes)){
if(rank(tmp/f) == tmprk){
facet_is_hyperplane = true;
result ^= i;
}
i++;
}
return result;
}
void populate_neighbors(){
const Matrix<Scalar> F = cache.get_facets(signature);
const Vector<Scalar>& apv(cache.all_positive_eq());
for(const auto& f : rows(F)){
if(!cache.facet_belongs_to_support(f)){
bool facet_is_hyperplane = false;
Bitset neighborS = neighbor_signature_from_facet(f, facet_is_hyperplane);
if(facet_is_hyperplane){
Vector<Scalar> neighborV = signature_to_vertex(hyperplanes, neighborS, apv);
Scalar apvdiff = neighborV[0]-vertex[0];
if(apvdiff<0){
upNeighbors[neighborV] = neighborS;
} else if (apvdiff==0){
if(lex_compare(neighborV, vertex) == 1){
upNeighbors[neighborV] = neighborS;
} else {
downNeighbors[neighborV] = neighborS;
}
} else {
downNeighbors[neighborV] = neighborS;
}
}
}
}
}
public:
Node& operator=(const Node& in){
signature = in.signature;
vertex = in.vertex;
upNeighbors = in.upNeighbors;
downNeighbors = in.downNeighbors;
return *this;
}
Node(const Matrix<Scalar>& hyp, const Bitset& sig, CacheType& c):
hyperplanes(hyp), signature(sig), cache(c) {
vertex = signature_to_vertex(hyperplanes, signature, cache.all_positive_eq());
populate_neighbors();
}
bool operator==(const Node& other) const {
return signature == other.signature;
}
bool has_jth_child(Int j) const {
return j < downNeighbors.size();
}
bool has_predecessor(const Node& pred) const {
if(!has_upneighbor()){ return false; }
const auto& front = upNeighbors.front();
return front.first == pred.vertex;
}
Node get_jth_child(Int j) const {
Int i = 0;
for(const auto& neighbor : downNeighbors){
if(i == j){
return Node(hyperplanes, neighbor.second, cache);
}
i++;
}
return *this;
}
Node get_predecessor(Int& j) const {
const auto& front = upNeighbors.front();
Node result(hyperplanes, front.second, cache);
j = 0;
for(const auto& neighbor : result.downNeighbors){
if(neighbor.second == signature){
break;
}
j++;
}
return result;
}
Int get_Delta() const {
return downNeighbors.size();
}
bool has_upneighbor() const {
return upNeighbors.size() > 0;
}
Matrix<Scalar> get_rays() const {
return cache.get_rays(signature);
}
const Bitset& get_signature() const {
return signature;
}
};
// Get a generic point in the [[SUPPORT]] of a HyperplaneArrangement, i.e. a
// point that does not lie on any of the [[HYPERPLANES]].
template<typename Scalar>
Vector<Scalar> get_generic_point(const Matrix<Scalar>& hyp, BigObject support){
// TODO: Check whether some of the hyperplanes contain support and ignore these.
// Make sure that generic point is not contained in one of the hyperplanes.
Matrix<Scalar> rays = support.give("RAYS | INPUT_RAYS");
Matrix<Scalar> lineality = support.give("LINEALITY_SPACE | INPUT_LINEALITY");
Matrix<Scalar> gens(rays / lineality);
Vector<Scalar> result(gens.cols());
for(const auto& row : rows(gens)){
result += rand() * row;
}
return result;
}
template<typename Scalar>
Bitset point_to_signature(const Vector<Scalar>& point, const Matrix<Scalar>& hyp, BigObject support){
Bitset signature;
const Int n = hyp.rows();
for( Int i = 0; i < n; ++i){
if(hyp[i] * point > 0){
signature.insert(i);
}
}
return signature;
}
} // namespace reverse_search_chamber_decomposition
} // namespace fan
} // namespace polymake
|
