source:git/gfanlib/gfanlib_mixedvolume.cpp@c987db

spielwiese
Last change on this file since c987db was 5443c1, checked in by Hans Schoenemann <hannes@…>, 4 years ago
porting: gcc-5 (Ubuntu 16.04.5 LTS)
• Property mode set to `100644`
File size: 8.0 KB
Line
1/*
2 * gfanlib_mixedvolume.cpp
3 *
4 *  Created on: Apr 10, 2016
5 *      Author: anders
6 */
7#include <stddef.h>
8#include "gfanlib_mixedvolume.h"
9
10#include "gfanlib_circuittableint.h"
11#include "gfanlib_tropicalhomotopy.h"
12
13#include "gfanlib_z.h"
14
15using namespace std;
16
17typedef gfan::CircuitTableInt32 mvtyp;
18
19namespace gfan{
20
21static Integer stringToInteger(char const *s)
22{
23        Integer ret;
24        while(*s)ret=Integer(10)*ret+((*(s++))-'0');
25        return ret;
26}
27
28
29class MVExceptionRethrow: public exception
30{
31  virtual const char* what() const throw()
32  {
33    return "Exception: Rethrown after detecting aborting flag in tropical homotopy.";
34  }
35}MVExceptionRethrow;
36
37
38class MVAmbientDimensionMismatch: public std::exception
39{
40  virtual const char* what() const throw()
41  {
42        return "Exception: The number of polytopes passed to the mixed volume computation does not match the ambient dimension.";
43  }
44}MVAmbientDimensionMismatch;
45
46
47static std::vector<Matrix<mvtyp> > convertTuple(std::vector<IntMatrix> const &tuple)
48{
49        vector<Matrix<mvtyp> > ret;
50        ret.reserve(tuple.size());
51        for(int i=0;i<tuple.size();i++)
52        {
53                Matrix<mvtyp> a(tuple[i].getHeight(),tuple[i].getWidth());
54                for(int j=0;j<a.getHeight();j++)
55                        for(int k=0;k<a.getWidth();k++)
56                                a[j][k]=tuple[i][j][k];
57                ret.push_back(a);
58        }
59        return ret;
60}
61
62//        mvtyp::Double
63        Integer mixedVolume(vector<IntMatrix> const &tuple, int nthreads, int steps)
64        {
65                //Zero-dimensional volumes are obtained by taking determinants of 0x0 matrices. Such determinants are 1. Hence the volume polynomial is constant 1. The coefficient of the product of the variables is 1.
66                //Alternative explanation. A generic polynomial system with zero variables has exactly 1 solution.
67                if(tuple.size()==0){return Integer(1);}
68
69                // For now, let's make in an error to pass the wrong number of polytopes. Of course we could define the
70                // mixed volume to be 0 in that case, but the code calling will almost certainly be wrong. It is safer to throw.
71                // for(const auto &x:tuple) // Range based for is not supported until gcc 4.6/
72                for(auto x=tuple.begin();x!=tuple.end();x++)
73                        if(x->getHeight()!=tuple.size())
74                                {throw MVAmbientDimensionMismatch;}
75
76                std::vector<Matrix<mvtyp> > tuple2=convertTuple(tuple);
77                //                assert(tuple2.size());
78                typedef SpecializedRTraverser<mvtyp,mvtyp::Double,mvtyp::Divisor> MySpecializedRTraverser;
79                vector<MySpecializedRTraverser> T1;
81                if(N==0)N=1;// Even if we do not parallelize, we still need one traverser.
82                T1.reserve(N);
83                vector<Traverser*> I;
84                for(int i=0;i<N;i++)T1.emplace_back(tuple2);
85                for(int i=0;i<N;i++)I.push_back(&(T1[i]));
86
87//                debug<<"Initialized.\n";
88
91                else
92                        traverse_simple(I[0]);
93
94                bool aborted=false;
95                for(int i=0;i<N;i++)aborted|=I[i]->aborting;
96
97                if(aborted)
98                {
99                        throw MVExceptionRethrow;
100                }
101
102                mvtyp::Double total;
103                int totalSteps=0;
104                for(int i=0;i<N;i++)
105                {
106//                        debug<<"#"<<((MySpecializedRTraverser*)(I[i]))->mixedVolume.toString()<<"Steps:"<<((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps<<"\n";
107                        //debug<<((SpecializedMTraverser*)(I[i]))->T.counter;
109                        totalSteps+=((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps;
110                }
111//                debug<<"Total:"<<(int)total.v<<"\n";
112//                debug<<"Totalsteps:"<<totalSteps<<"\n";
113
114
115                return stringToInteger(total.toString().c_str());
116        }
117
118
119        namespace MixedVolumeExamples{
120                vector<IntMatrix> cyclic(int n)
121                {
122                        vector<IntMatrix> ret;
123                        for(int i=1;i<n;i++)
124                        {
125                                IntMatrix m(n,n);
126                                for(int y=0;y<n;y++)
127                                        for(int x=0;x<n;x++)
128                                                m[y][x]=((x-y+n)%n)<i;
129                                ret.push_back(m);
130                        }
131
132                        IntMatrix last(n,2);
133                        for(int y=0;y<n;y++)last[y][0]=1;
134                        ret.push_back(last);
135
136                        return ret;
137                }
138                vector<IntMatrix> noon(int n)
139                {
140                        vector<IntMatrix> ret;
141                        for(int i=0;i<n;i++)
142                        {
143                                IntMatrix m(n,n+1);
144                                for(int y=0;y<n-1;y++)
145                                        m[y+(i<=y)][y]=2;
146                                for(int x=0;x<n;x++)
147                                        m[i][x]=1;
148                                ret.push_back(m);
149                        }
150                        return ret;
151                }
152                vector<IntMatrix> chandra(int n)
153                {
154                        vector<IntMatrix> ret;
155                        for(int i=0;i<n;i++)
156                        {
157                                IntMatrix m(n,n+1);
158                                for(int y=0;y<n-1;y++)
159                                        m[y][y+1]=1;
160                                for(int x=0;x<n;x++)
161                                        m[i][x]+=1;
162                                ret.push_back(m);
163                        }
164                        return ret;
165                }
166                vector<IntMatrix> katsura(int n)
167                {
168                        n++;
169                        vector<IntMatrix> ret;
170                        for(int i=0;i<n-1;i++)
171                        {
172                                IntMatrix m(n,n+1-((i+1)/2));
173                                for(int y=0;y<n-((i+1)/2);y++)
174                                {
175                                        m[n-1-y][y]=1;
176                                        m[(n-1-y-i)>0 ? (n-1-y-i) : -(n-1-y-i)][y]+=1;
177                                }
178                                m[i][m.getWidth()-1]=1;
179                                ret.push_back(m);
180                        }
181                        ret.push_back(combineLeftRight(IntMatrix::identity(n),IntMatrix(n,1)));
182                        return ret;
183                }
184                vector<IntMatrix> gaukwa(int n)
185                {
186                        vector<IntMatrix> ret;
187                        for(int i=0;i<2*n;i++)
188                                ret.push_back(combineLeftRight(combineOnTop(IntMatrix::identity(n),i*IntMatrix::identity(n)),IntMatrix(2*n,1)));
189                        return ret;
190                }
191                vector<IntMatrix> eco(int n)
192                {
193                        vector<IntMatrix> ret;
194                        for(int i=0;i<n-1;i++)
195                        {
196                                IntMatrix m(n,n-i);
197                                for(int y=0;y<n-(i+1);y++)
198                                {
199                                        m[y+i][y]=1;
200                                        m[n-1][y]=1;
201                                        if(y)m[y-1][y]=1;
202                                }
203                                ret.push_back(m);
204                        }
205                        ret.push_back(combineLeftRight(combineOnTop(IntMatrix::identity(n-1),IntMatrix(1,n-1)),IntMatrix(n,1)));
206                        return ret;
207                }
208        }
209
210}
Note: See TracBrowser for help on using the repository browser.