1 #include "parser.h" 2 #include "printer.h" 3 #include "polynomial.h" 4 #include "division.h" 5 #include "buchberger.h" 6 #include "wallideal.h" 7 #include "lp.h" 8 #include "reversesearch.h" 9 #include "termorder.h" 10 #include "gfanapplication.h" 11 #include "wallideal.h" 12 #include "polyhedralfan.h" 13 #include "linalg.h" 14 #include "polymakefile.h" 15 #include "tropical_weildivisor.h" 16 #include "log.h" 17 #include "symmetry.h" 18 19 class TropicalWeilDivisorApplication : public GFanApplication 20 { 21 StringOption inputOption1; 22 StringOption inputOption2; 23 public: helpText()24 const char *helpText() 25 { 26 return "This program computes the tropical Weil divisor of piecewise linear (or tropical rational) function on a tropical k-cycle. See the Gfan manual for more information.\n"; 27 } TropicalWeilDivisorApplication()28 TropicalWeilDivisorApplication(): 29 inputOption1("-i1","Specify the name of the Polymake input file containing the k-cycle.","polymake1.out"), 30 inputOption2("-i2","Specify the name of the Polymake input file containing the piecewise linear function.","polymake2.out") 31 { 32 registerOptions(); 33 } 34 name()35 const char *name() 36 { 37 return "_tropicalweildivisor"; 38 } 39 readLinealitySpace(const char * filename)40 IntegerVectorList readLinealitySpace(const char *filename) 41 { 42 PolymakeFile inFile; 43 inFile.open(filename); 44 int n=inFile.readCardinalProperty("AMBIENT_DIM"); 45 int homog=inFile.readCardinalProperty("LINEALITY_DIM"); 46 IntegerMatrix rays=inFile.readMatrixProperty("LINEALITY_SPACE",homog,n); 47 return rays.getRows(); 48 } readRays(const char * filename)49 IntegerVectorList readRays(const char *filename) 50 { 51 PolymakeFile inFile; 52 inFile.open(filename); 53 int n=inFile.readCardinalProperty("AMBIENT_DIM"); 54 int nRays=inFile.readCardinalProperty("N_RAYS"); 55 IntegerMatrix rays=inFile.readMatrixProperty("RAYS",nRays,n); 56 return rays.getRows(); 57 } readFan(const char * filename,bool readRayValues,bool readLinearForms,bool readMultiplicities)58 PolyhedralFan readFan(const char *filename, bool readRayValues, bool readLinearForms, bool readMultiplicities) 59 { 60 PolymakeFile inFile; 61 inFile.open(filename); 62 63 int n=inFile.readCardinalProperty("AMBIENT_DIM"); 64 int nRays=inFile.readCardinalProperty("N_RAYS"); 65 int lindim=inFile.readCardinalProperty("LINEALITY_DIM"); 66 IntegerMatrix rays=inFile.readMatrixProperty("RAYS",nRays,n); 67 vector<list<int> > cones=inFile.readMatrixIncidenceProperty("MAXIMAL_CONES"); 68 69 IntegerMatrix multiplicities; 70 if(readMultiplicities) 71 multiplicities=inFile.readMatrixProperty("MULTIPLICITIES",cones.size(),1); 72 73 PolyhedralFan ret(n); 74 75 IntegerVectorList equations=inFile.readMatrixProperty("LINEALITY_SPACE",lindim,n).getRows(); 76 77 IntegerMatrix linearForms(cones.size(),n); 78 if(readLinearForms)linearForms=inFile.readMatrixProperty("LINEAR_FORMS",cones.size(),n); 79 80 IntegerMatrix rayValues(nRays,1); 81 if(readRayValues)rayValues=inFile.readMatrixProperty("RAY_VALUES",nRays,1); 82 83 IntegerMatrix linealityValues(lindim,1); 84 if(readRayValues)linealityValues=inFile.readMatrixProperty("LINEALITY_VALUES",lindim,1); 85 86 87 int I=0; 88 for(vector<list<int> >::const_iterator i=cones.begin();i!=cones.end();i++,I++) 89 { 90 IntegerVectorList coneGenerators; 91 FieldMatrix LFA(Q,lindim+i->size(),n); 92 FieldVector LFB(Q,lindim+i->size()+n); 93 int J=0; 94 for(list<int>::const_iterator j=i->begin();j!=i->end();j++,J++) 95 { 96 coneGenerators.push_back(rays[*j]); 97 LFA[J]=integerVectorToFieldVector(rays[*j],Q); 98 LFB[J]=Q.zHomomorphism(rayValues[*j][0]); 99 } 100 int K=0; 101 for(IntegerVectorList::const_iterator j=equations.begin();j!=equations.end();j++,J++,K++) 102 { 103 LFA[J]=integerVectorToFieldVector(*j,Q); 104 LFB[J]=Q.zHomomorphism(linealityValues[K][0]); 105 } 106 if(readRayValues) 107 { 108 AsciiPrinter P(Stderr);//log0 P<<LFA<<LFB; 109 110 FieldVector LFX=LFA.solver().canonicalize(LFB); 111 if(LFX.subvector(0,LFX.size()-n).isZero()) 112 { 113 linearForms[I]=fieldVectorToIntegerVector(LFX.subvector(LFX.size()-n,LFX.size())); 114 } 115 else 116 { 117 cerr<<"Values on cone are not linear" <<endl; 118 assert(0); 119 } 120 } 121 //AsciiPrinter P(Stderr); 122 PolyhedralCone d(coneGenerators,equations,n); 123 d.canonicalize(); 124 //d.print(P); 125 IntegerVectorList inequalities=d.extremeRays(); 126 IntegerVectorList equations2; 127 IntegerMatrix A=rowsToIntegerMatrix(coneGenerators,n); 128 IntegerMatrix B=rowsToIntegerMatrix(equations,n);///!!!! 129 FieldMatrix M=combineOnTop(integerMatrixToFieldMatrix(A,Q),integerMatrixToFieldMatrix(B,Q)); 130 FieldMatrix kerM=M.reduceAndComputeKernel(); 131 for(int i=0;i<kerM.getHeight();i++) 132 { 133 equations2.push_back(kerM[i].primitive()); 134 } 135 PolyhedralCone c(inequalities,equations2,n); 136 c.canonicalize(); 137 if(readLinearForms || readRayValues)c.setLinearForm(linearForms[I]); 138 139 if(readMultiplicities) 140 c.setMultiplicity(multiplicities[I][0]); 141 142 ret.insert(c); 143 } 144 return ret; 145 } 146 main()147 int main() 148 { 149 // LpSolver::printList(Stderr); 150 // lpSetSolver("cddgmp"); 151 152 153 PolyhedralFan f=readFan(inputOption1.getValue(),false,false,true); 154 PolyhedralFan g=readFan(inputOption2.getValue(),true,false,false); 155 156 IntegerVectorList originalRays=readRays(inputOption1.getValue()); 157 IntegerVectorList linealitySpace=readLinealitySpace(inputOption1.getValue()); 158 /* FileParser P(Stdin); 159 160 PolynomialRing R=P.parsePolynomialRing(); 161 162 Polynomial p=P.parsePolynomial(R); 163 */ 164 165 PolyhedralFan divisor=weilDivisor(f,g); 166 167 AsciiPrinter Q(stdout); 168 169 SymmetryGroup sym(f.getAmbientDimension()); 170 vector<string> comments=f.renamingStrings(divisor.getRaysInPrintingOrder(&sym),originalRays,linealitySpace,&sym); 171 // divisor.printWithIndices(&Q,true,0,false,false,false,false,&comments); 172 divisor.printWithIndices(&Q,FPF_default|FPF_multiplicities,0,&comments); 173 174 /* FileParser P(Stdin); 175 176 PolynomialRing R=P.parsePolynomialRing(); 177 PolynomialSet s=P.parsePolynomialSet(R); 178 179 Polynomial sum(R); 180 for(PolynomialSet::const_iterator i=s.begin();i!=s.end();i++) 181 sum+=*i; 182 183 AsciiPrinter(Stdout).printPolynomial(sum); 184 */ 185 186 187 return 0; 188 } 189 }; 190 191 static TropicalWeilDivisorApplication theApplication; 192