source: git/gfanlib/gfanlib_zcone.cpp @ 26b713

jengelh-datetimespielwiese
Last change on this file since 26b713 was 26b713, checked in by Yue Ren <ren@…>, 10 years ago
chg: updated gfanlib package to newest version in master
  • Property mode set to 100644
File size: 33.4 KB
Line 
1/*
2 * lib_cone.cpp
3 *
4 *  Created on: Sep 29, 2010
5 *      Author: anders
6 */
7
8#include "gfanlib_zcone.h"
9
10#include <vector>
11#include <set>
12
13#include "setoper.h"
14#include "cdd.h"
15
16namespace gfan{
17
18  static void cddinitGmp()
19  {
20    static bool initialized;
21    if(!initialized)
22      {
23        dd_set_global_constants();  /* First, this must be called. */
24        initialized=true;
25      }
26  }
27
28
29class LpSolver
30{
31  static dd_MatrixPtr ZMatrix2MatrixGmp(ZMatrix const &g, dd_ErrorType *Error)
32  {
33    int n=g.getWidth();
34    dd_MatrixPtr M=NULL;
35    dd_rowrange m_input,i;
36    dd_colrange d_input,j;
37    dd_RepresentationType rep=dd_Inequality;
38    dd_boolean found=dd_FALSE, newformat=dd_FALSE, successful=dd_FALSE;
39    char command[dd_linelenmax], comsave[dd_linelenmax];
40    dd_NumberType NT;
41
42    (*Error)=dd_NoError;
43
44    rep=dd_Inequality; newformat=dd_TRUE;
45
46    m_input=g.getHeight();
47    d_input=n+1;
48
49    NT=dd_Rational;
50
51    M=dd_CreateMatrix(m_input, d_input);
52    M->representation=rep;
53    M->numbtype=NT;
54
55    for (i = 0; i < m_input; i++) {
56      dd_set_si(M->matrix[i][0],0);
57      for (j = 1; j < d_input; j++) {
58        g[i][j-1].setGmp(mpq_numref(M->matrix[i][j]));
59        mpz_init_set_ui(mpq_denref(M->matrix[i][j]), 1);
60        mpq_canonicalize(M->matrix[i][j]);
61      }
62    }
63
64    successful=dd_TRUE;
65
66    return M;
67  }
68  static dd_MatrixPtr ZMatrix2MatrixGmp(ZMatrix const &inequalities, ZMatrix const &equations, dd_ErrorType *err)
69  {
70    ZMatrix g=inequalities;
71    g.append(equations);
72    int numberOfInequalities=inequalities.getHeight();
73    int numberOfRows=g.getHeight();
74    dd_MatrixPtr A=NULL;
75    cddinitGmp();
76    A=ZMatrix2MatrixGmp(g, err);
77    for(int i=numberOfInequalities;i<numberOfRows;i++)
78      set_addelem(A->linset,i+1);
79    return A;
80  }
81  static ZMatrix getConstraints(dd_MatrixPtr A, bool returnEquations)
82  {
83    int rowsize=A->rowsize;
84    int n=A->colsize-1;
85
86    ZMatrix ret(0,n);
87    for(int i=0;i<rowsize;i++)
88      {
89        bool isEquation=set_member(i+1,A->linset);
90        if(isEquation==returnEquations)
91          {
92            QVector v(n);
93            for(int j=0;j<n;j++)v[j]=Rational(A->matrix[i][j+1]);
94            ret.appendRow(QToZVectorPrimitive(v));
95          }
96      }
97    return ret;
98  }
99  static bool isFacet(ZMatrix const &g, int index)
100  {
101    bool ret;
102    dd_MatrixPtr M=NULL,M2=NULL,M3=NULL;
103    dd_colrange d;
104    dd_ErrorType err=dd_NoError;
105    dd_rowset redrows,linrows,ignoredrows, basisrows;
106    dd_colset ignoredcols, basiscols;
107    dd_DataFileType inputfile;
108    FILE *reading=NULL;
109
110    cddinitGmp();
111
112    M=ZMatrix2MatrixGmp(g, &err);
113    if (err!=dd_NoError) goto _L99;
114
115    d=M->colsize;
116
117    static dd_Arow temp;
118    dd_InitializeArow(g.getWidth()+1,&temp);
119
120    ret= !dd_Redundant(M,index+1,temp,&err);
121
122    dd_FreeMatrix(M);
123    dd_FreeArow(g.getWidth()+1,temp);
124
125    if (err!=dd_NoError) goto _L99;
126    return ret;
127   _L99:
128    assert(0);
129    return false;
130  }
131
132  /*
133    Heuristic for checking if inequality of full dimensional cone is a
134    facet. If the routine returns true then the inequality is a
135    facet. If it returns false it is unknown.
136   */
137  static bool fastIsFacetCriterion(ZMatrix const &normals, int i)
138  {
139    int n=normals.getWidth();
140    for(int j=0;j<n;j++)
141      if(normals[i][j].sign()!=0)
142        {
143          int sign=normals[i][j].sign();
144          bool isTheOnly=true;
145          for(int k=0;k<normals.getHeight();k++)
146            if(k!=i)
147              {
148                if(normals[i][j].sign()==sign)
149                  {
150                    isTheOnly=false;
151                    break;
152                  }
153              }
154          if(isTheOnly)return true;
155        }
156    return false;
157  }
158
159  static bool fastIsFacet(ZMatrix const &normals, int i)
160  {
161    if(fastIsFacetCriterion(normals,i))return true;
162    return isFacet(normals,i);
163  }
164
165  class MyHashMap
166  {
167    typedef std::vector<std::set<ZVector> > Container;
168    Container container;
169    int tableSize;
170  public:
171    class iterator
172    {
173      class MyHashMap &hashMap;
174      int index; // having index==-1 means that we are before/after the elements.
175      std::set<ZVector>::iterator i;
176    public:
177      bool operator++()
178      {
179        if(index==-1)goto search;
180        i++;
181        while(i==hashMap.container[index].end())
182          {
183            search:
184            index++;
185            if(index>=hashMap.tableSize){
186              index=-1;
187              return false;
188            }
189            i=hashMap.container[index].begin();
190          }
191        return true;
192      }
193      ZVector const & operator*()const
194      {
195        return *i;
196      }
197      ZVector operator*()
198      {
199        return *i;
200      }
201      iterator(MyHashMap &hashMap_):
202        hashMap(hashMap_)
203        {
204          index=-1;
205        }
206    };
207    unsigned int function(const ZVector &v)
208    {
209      unsigned int ret=0;
210      int n=v.size();
211      for(int i=0;i<n;i++)
212        ret=(ret<<3)+(ret>>29)+v.UNCHECKEDACCESS(i).hashValue();
213      return ret%tableSize;
214    }
215    MyHashMap(int tableSize_):
216      container(tableSize_),
217      tableSize(tableSize_)
218      {
219        assert(tableSize_>0);
220      }
221    void insert(const ZVector &v)
222    {
223      container[function(v)].insert(v);
224    }
225    void erase(ZVector const &v)
226    {
227      container[function(v)].erase(v);
228    }
229    iterator begin()
230    {
231      iterator ret(*this);
232      ++    ret;
233      return ret;
234    }
235    int size()
236    {
237      iterator i=begin();
238      int ret=0;
239      do{ret++;}while(++i);
240      return ret;
241    }
242  };
243
244
245  static ZMatrix normalizedWithSumsAndDuplicatesRemoved(ZMatrix const &a)
246  {
247    // TODO: write a version of this function which will work faster if the entries fit in 32bit
248    if(a.getHeight()==0)return a;
249    int n=a.getWidth();
250    ZVector temp1(n);
251//    ZVector temp2(n);
252    ZMatrix ret(0,n);
253    MyHashMap b(a.getHeight());
254
255    for(int i=0;i<a.getHeight();i++)
256      {
257        assert(!(a[i].isZero()));
258        b.insert(a[i].normalized());
259      }
260
261      {
262        MyHashMap::iterator i=b.begin();
263
264        do
265          {
266            MyHashMap::iterator j=i;
267            while(++j)
268              {
269                ZVector const &I=*i;
270                ZVector const &J=*j;
271                for(int k=0;k<n;k++)temp1[k]=I.UNCHECKEDACCESS(k)+J.UNCHECKEDACCESS(k);
272//                normalizedLowLevel(temp1,temp2);
273//                b.erase(temp2);//this can never remove *i or *j
274                b.erase(temp1.normalized());//this can never remove *i or *j
275              }
276          }
277          while(++i);
278      }
279    ZMatrix original(0,n);
280    {
281      MyHashMap::iterator i=b.begin();
282      do
283        {
284          original.appendRow(*i);
285        }
286      while(++i);
287    }
288
289    for(int i=0;i!=original.getHeight();i++)
290      for(int j=0;j!=a.getHeight();j++)
291        if(!dependent(original[i],a[j]))
292            {
293              ZVector const &I=original[i];
294              ZVector const &J=a[j];
295              for(int k=0;k<n;k++)temp1[k]=I.UNCHECKEDACCESS(k)+J.UNCHECKEDACCESS(k);
296//              normalizedLowLevel(temp1,temp2);
297//              b.erase(temp2);//this can never remove *i or *j
298              b.erase(temp1.normalized());//this can never remove *i or *j
299            }
300        {
301          MyHashMap::iterator i=b.begin();
302          do
303          {
304            ZVector temp=*i;
305            ret.appendRow(temp);
306          }
307          while(++i);
308      }
309    return ret;
310  }
311public:
312  static ZMatrix fastNormals(ZMatrix const &inequalities)
313  {
314    ZMatrix normals=normalizedWithSumsAndDuplicatesRemoved(inequalities);
315    for(int i=0;i!=normals.getHeight();i++)
316      if(!fastIsFacet(normals,i))
317        {
318          normals[i]=normals[normals.getHeight()-1];
319          normals.eraseLastRow();
320          i--;
321        }
322    return normals;
323  }
324  void removeRedundantRows(ZMatrix &inequalities, ZMatrix &equations, bool removeInequalityRedundancies)
325  {
326    cddinitGmp();
327    int numberOfEqualities=equations.getHeight();
328    int numberOfInequalities=inequalities.getHeight();
329    int numberOfRows=numberOfEqualities+numberOfInequalities;
330
331    if(numberOfRows==0)return;//the full space, so description is already irredundant
332
333    dd_rowset r=NULL;
334    ZMatrix g=inequalities;
335    g.append(equations);
336
337    dd_LPSolverType solver=dd_DualSimplex;
338    dd_MatrixPtr A=NULL;
339    dd_ErrorType err=dd_NoError;
340
341    A=ZMatrix2MatrixGmp(g,&err);
342    if (err!=dd_NoError) goto _L99;
343
344    for(int i=numberOfInequalities;i<numberOfRows;i++)
345      set_addelem(A->linset,i+1);
346
347    A->objective=dd_LPmax;
348
349    dd_rowset impl_linset;
350    dd_rowset redset;
351    dd_rowindex newpos;
352
353    if(removeInequalityRedundancies)
354      dd_MatrixCanonicalize(&A, &impl_linset, &redset, &newpos, &err);
355    else
356      dd_MatrixCanonicalizeLinearity(&A, &impl_linset, &newpos, &err);
357
358    if (err!=dd_NoError) goto _L99;
359
360    {
361      int n=A->colsize-1;
362      equations=ZMatrix(0,n);     //TODO: the number of rows needed is actually known
363      inequalities=ZMatrix(0,n);  //is known by set_card(). That might save some copying.
364
365      {
366        int rowsize=A->rowsize;
367        QVector point(n);
368        for(int i=0;i<rowsize;i++)
369          {
370            for(int j=0;j<n;j++)point[j]=Rational(A->matrix[i][j+1]);
371            ((set_member(i+1,A->linset))?equations:inequalities).appendRow(QToZVectorPrimitive(point));
372          }
373      }
374      assert(set_card(A->linset)==equations.getHeight());
375      assert(A->rowsize==equations.getHeight()+inequalities.getHeight());
376
377      set_free(impl_linset);
378      if(removeInequalityRedundancies)
379        set_free(redset);
380      free(newpos);
381
382      dd_FreeMatrix(A);
383      return;
384    }
385   _L99:
386    assert(!"Cddlib reported error when called by Gfanlib.");
387  }
388  ZVector relativeInteriorPoint(const ZMatrix &inequalities, const ZMatrix &equations)
389  {
390    QVector retUnscaled(inequalities.getWidth());
391    cddinitGmp();
392    int numberOfEqualities=equations.getHeight();
393    int numberOfInequalities=inequalities.getHeight();
394    int numberOfRows=numberOfEqualities+numberOfInequalities;
395
396    dd_rowset r=NULL;
397    ZMatrix g=inequalities;
398    g.append(equations);
399
400    dd_LPSolverType solver=dd_DualSimplex;
401    dd_MatrixPtr A=NULL;
402    dd_ErrorType err=dd_NoError;
403
404    A=ZMatrix2MatrixGmp(g,&err);
405    if (err!=dd_NoError) goto _L99;
406    {
407      dd_LPSolutionPtr lps1;
408      dd_LPPtr lp,lp1;
409
410      for(int i=0;i<numberOfInequalities;i++)
411        dd_set_si(A->matrix[i][0],-1);
412      for(int i=numberOfInequalities;i<numberOfRows;i++)
413        set_addelem(A->linset,i+1);
414
415      A->objective=dd_LPmax;
416      lp=dd_Matrix2LP(A, &err);
417      if (err!=dd_NoError) goto _L99;
418
419      lp1=dd_MakeLPforInteriorFinding(lp);
420      dd_LPSolve(lp1,solver,&err);
421      if (err!=dd_NoError) goto _L99;
422
423      lps1=dd_CopyLPSolution(lp1);
424
425      assert(!dd_Negative(lps1->optvalue));
426
427      for (int j=1; j <(lps1->d)-1; j++)
428        retUnscaled[j-1]=Rational(lps1->sol[j]);
429
430      dd_FreeLPData(lp);
431      dd_FreeLPSolution(lps1);
432      dd_FreeLPData(lp1);
433      dd_FreeMatrix(A);
434      return QToZVectorPrimitive(retUnscaled);
435    }
436_L99:
437    assert(0);
438    return QToZVectorPrimitive(retUnscaled);
439  }
440  void dual(ZMatrix const &inequalities, ZMatrix const &equations, ZMatrix &dualInequalities, ZMatrix &dualEquations)
441  {
442    int result;
443
444    dd_MatrixPtr A=NULL;
445    dd_ErrorType err=dd_NoError;
446
447    cddinitGmp();
448
449    A=ZMatrix2MatrixGmp(inequalities, equations, &err);
450
451    dd_PolyhedraPtr poly;
452    poly=dd_DDMatrix2Poly2(A, dd_LexMin, &err);
453
454    if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) assert(0);
455
456    dd_MatrixPtr      A2=dd_CopyGenerators(poly);
457
458    dualInequalities=getConstraints(A2,false);
459    dualEquations=getConstraints(A2,true);
460
461    dd_FreeMatrix(A2);
462    dd_FreeMatrix(A);
463    dd_FreePolyhedra(poly);
464
465    return;
466   _L99:
467    assert(0);
468  }
469  // this procedure is take from cddio.c.
470  static void dd_ComputeAinc(dd_PolyhedraPtr poly)
471  {
472  /* This generates the input incidence array poly->Ainc, and
473     two sets: poly->Ared, poly->Adom.
474  */
475    dd_bigrange k;
476    dd_rowrange i,m1;
477    dd_colrange j;
478    dd_boolean redundant;
479    dd_MatrixPtr M=NULL;
480    mytype sum,temp;
481
482    dd_init(sum); dd_init(temp);
483    if (poly->AincGenerated==dd_TRUE) goto _L99;
484
485    M=dd_CopyOutput(poly);
486    poly->n=M->rowsize;
487    m1=poly->m1;
488
489    /* this number is same as poly->m, except when
490        poly is given by nonhomogeneous inequalty:
491        !(poly->homogeneous) && poly->representation==Inequality,
492        it is poly->m+1.   See dd_ConeDataLoad.
493     */
494    poly->Ainc=(set_type*)calloc(m1, sizeof(set_type));
495    for(i=1; i<=m1; i++) set_initialize(&(poly->Ainc[i-1]),poly->n);
496    set_initialize(&(poly->Ared), m1);
497    set_initialize(&(poly->Adom), m1);
498
499    for (k=1; k<=poly->n; k++){
500      for (i=1; i<=poly->m; i++){
501        dd_set(sum,dd_purezero);
502        for (j=1; j<=poly->d; j++){
503          dd_mul(temp,poly->A[i-1][j-1],M->matrix[k-1][j-1]);
504          dd_add(sum,sum,temp);
505        }
506        if (dd_EqualToZero(sum)) {
507          set_addelem(poly->Ainc[i-1], k);
508        }
509      }
510      if (!(poly->homogeneous) && poly->representation==dd_Inequality){
511        if (dd_EqualToZero(M->matrix[k-1][0])) {
512          set_addelem(poly->Ainc[m1-1], k);  /* added infinity inequality (1,0,0,...,0) */
513        }
514      }
515    }
516
517    for (i=1; i<=m1; i++){
518      if (set_card(poly->Ainc[i-1])==M->rowsize){
519        set_addelem(poly->Adom, i);
520      }
521    }
522    for (i=m1; i>=1; i--){
523      if (set_card(poly->Ainc[i-1])==0){
524        redundant=dd_TRUE;
525        set_addelem(poly->Ared, i);
526      }else {
527        redundant=dd_FALSE;
528        for (k=1; k<=m1; k++) {
529          if (k!=i && !set_member(k, poly->Ared)  && !set_member(k, poly->Adom) &&
530              set_subset(poly->Ainc[i-1], poly->Ainc[k-1])){
531            if (!redundant){
532              redundant=dd_TRUE;
533            }
534            set_addelem(poly->Ared, i);
535          }
536        }
537      }
538    }
539    dd_FreeMatrix(M);
540    poly->AincGenerated=dd_TRUE;
541  _L99:;
542    dd_clear(sum);  dd_clear(temp);
543  }
544
545
546  std::vector<std::vector<int> > extremeRaysInequalityIndices(const ZMatrix &inequalities)
547  {
548    int dim2=inequalities.getHeight();
549    if(dim2==0)return std::vector<std::vector<int> >();
550    int dimension=inequalities.getWidth();
551
552    dd_MatrixPtr A=NULL;
553    dd_ErrorType err=dd_NoError;
554
555    cddinitGmp();
556    A=ZMatrix2MatrixGmp(inequalities, &err);
557
558    dd_PolyhedraPtr poly;
559    poly=dd_DDMatrix2Poly2(A, dd_LexMin, &err);
560
561    if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) assert(0);
562    if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
563
564    std::vector<std::vector<int> > ret;
565
566    /*
567      How do we interpret the cddlib output?  For a long ting gfan has
568      been using poly->n as the number of rays of the cone and thus
569      returned sets of indices that actually gave the lineality
570      space. The mistake was then caught later in PolyhedralCone. On Feb
571      17 2009 gfan was changed to check the length of each set to make
572      sure that it does not specify the lineality space and only return
573      those sets giving rise to rays.  This does not seem to be the best
574      strategy and might even be wrong.
575     */
576
577
578    for (int k=1; k<=poly->n; k++)
579      {
580        int length=0;
581        for (int i=1; i<=poly->m1; i++)
582          if(set_member(k,poly->Ainc[i-1]))length++;
583        if(length!=dim2)
584          {
585            std::vector<int> v(length);
586            int j=0;
587            for (int i=1; i<=poly->m1; i++)
588              if(set_member(k,poly->Ainc[i-1]))v[j++]=i-1;
589            ret.push_back(v);
590          }
591      }
592
593    dd_FreeMatrix(A);
594    dd_FreePolyhedra(poly);
595
596    return ret;
597   _L99:
598    assert(0);
599    return std::vector<std::vector<int> >();
600  }
601
602};
603
604LpSolver lpSolver;
605
606bool ZCone::isInStateMinimum(int s)const
607{
608  return state>=s;
609}
610
611
612bool operator<(ZCone const &a, ZCone const &b)
613{
614  assert(a.state>=3);
615  assert(b.state>=3);
616
617  if(a.n<b.n)return true;
618  if(a.n>b.n)return false;
619
620  if(a.equations<b.equations)return true;
621  if(b.equations<a.equations)return false;
622
623  if(a.inequalities<b.inequalities)return true;
624  if(b.inequalities<a.inequalities)return false;
625
626  return false;
627}
628
629
630bool operator!=(ZCone const &a, ZCone const &b)
631{
632  return (a<b)||(b<a);
633}
634
635
636void ZCone::ensureStateAsMinimum(int s)const
637{
638  if((state<1) && (s==1))
639    {
640     {
641        QMatrix m=ZToQMatrix(equations);
642        m.reduce();
643        m.removeZeroRows();
644
645        ZMatrix newInequalities(0,inequalities.getWidth());
646        for(int i=0;i<inequalities.getHeight();i++)
647          {
648            QVector w=ZToQVector(inequalities[i]);
649            w=m.canonicalize(w);
650            if(!w.isZero())
651              newInequalities.appendRow(QToZVectorPrimitive(w));
652          }
653
654        inequalities=newInequalities;
655        inequalities.sortAndRemoveDuplicateRows();
656        equations=QToZMatrixPrimitive(m);
657      }
658
659      if(!(preassumptions&PCP_impliedEquationsKnown))
660      if(inequalities.getHeight()>1)//there can be no implied equation unless we have at least two inequalities
661        lpSolver.removeRedundantRows(inequalities,equations,false);
662
663      assert(inequalities.getWidth()==equations.getWidth());
664      }
665  if((state<2) && (s>=2) && !(preassumptions&PCP_facetsKnown))
666    {
667/*       if(inequalities.size()>25)
668         {
669          IntegerVectorList h1;
670          IntegerVectorList h2;
671          bool a=false;
672          for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
673            {
674              if(a)
675                h1.push_back(*i);
676              else
677                h2.push_back(*i);
678              a=!a;
679            }
680          PolyhedralCone c1(h1,equations);
681          PolyhedralCone c2(h2,equations);
682          c1.ensureStateAsMinimum(2);
683          c2.ensureStateAsMinimum(2);
684          inequalities=c1.inequalities;
685          for(IntegerVectorList::const_iterator i=c2.inequalities.begin();i!=c2.inequalities.end();i++)
686            inequalities.push_back(*i);
687        }
688*/
689      if(equations.getHeight())
690        {
691          QMatrix m=ZToQMatrix(equations);
692          m.reduce();
693          m.REformToRREform();
694          ZMatrix inequalities2(0,equations.getWidth());
695          for(int i=0;i<inequalities.getHeight();i++)
696            {
697              inequalities2.appendRow(QToZVectorPrimitive(m.canonicalize(ZToQVector(inequalities[i]))));
698            }
699          inequalities=LpSolver::fastNormals(inequalities2);
700          goto noFallBack;
701        fallBack://alternativ (disabled)
702          lpSolver.removeRedundantRows(inequalities,equations,true);
703        noFallBack:;
704        }
705      else
706        inequalities=LpSolver::fastNormals(inequalities);
707    }
708  if((state<3) && (s>=3))
709    {
710      QMatrix equations2=ZToQMatrix(equations);
711      equations2.reduce(false,false,true);
712      equations2.REformToRREform();
713      for(int i=0;i<inequalities.getHeight();i++)
714        {
715          inequalities[i]=QToZVectorPrimitive(equations2.canonicalize(ZToQVector(inequalities[i])));
716        }
717      inequalities.sortRows();
718      equations=QToZMatrixPrimitive(equations2);
719    }
720  if(state<s)
721    state=s;
722}
723
724std::ostream &operator<<(std::ostream &f, ZCone const &c)
725{
726  f<<"Ambient dimension:"<<c.n<<std::endl;
727  f<<"Inequalities:"<<std::endl;
728  f<<c.inequalities<<std::endl;
729  f<<"Equations:"<<std::endl;
730  f<<c.equations<<std::endl;
731}
732
733
734ZCone::ZCone(int ambientDimension):
735  inequalities(ZMatrix(0,ambientDimension)),
736  equations(ZMatrix(0,ambientDimension)),
737  n(ambientDimension),
738  state(1),
739  preassumptions(PCP_impliedEquationsKnown|PCP_facetsKnown),
740  multiplicity(1),
741  haveExtremeRaysBeenCached(false),
742  linearForms(ZMatrix(0,ambientDimension))
743{
744}
745
746
747ZCone::ZCone(ZMatrix const &inequalities_, ZMatrix const &equations_, int preassumptions_):
748  inequalities(inequalities_),
749  equations(equations_),
750  state(0),
751  preassumptions(preassumptions_),
752  multiplicity(1),
753  haveExtremeRaysBeenCached(false),
754  n(inequalities_.getWidth()),
755  linearForms(ZMatrix(0,inequalities_.getWidth()))
756  {
757  assert(preassumptions_<4);//OTHERWISE WE ARE DOING SOMETHING STUPID LIKE SPECIFYING AMBIENT DIMENSION
758  assert(equations_.getWidth()==n);
759  ensureStateAsMinimum(1);
760}
761
762void ZCone::canonicalize()
763{
764  ensureStateAsMinimum(3);
765}
766
767void ZCone::findFacets()
768{
769  ensureStateAsMinimum(2);
770}
771
772ZMatrix ZCone::getFacets()const
773{
774  ensureStateAsMinimum(2);
775  return inequalities;
776}
777
778void ZCone::findImpliedEquations()
779{
780  ensureStateAsMinimum(1);
781}
782
783ZMatrix ZCone::getImpliedEquations()const
784{
785  ensureStateAsMinimum(1);
786  return equations;
787}
788
789ZVector ZCone::getRelativeInteriorPoint()const
790{
791  ensureStateAsMinimum(1);
792//  assert(state>=1);
793
794  return lpSolver.relativeInteriorPoint(inequalities,equations);
795}
796
797ZVector ZCone::getUniquePoint()const
798{
799  ZMatrix rays=extremeRays();
800  ZVector ret(n);
801  for(int i=0;i<rays.getHeight();i++)
802    ret+=rays[i];
803
804  return ret;
805}
806
807ZVector ZCone::getUniquePointFromExtremeRays(ZMatrix const &extremeRays)const
808{
809  ZVector ret(n);
810  for(int i=0;i<extremeRays.getHeight();i++)
811    if(contains(extremeRays[i]))ret+=extremeRays[i];
812  return ret;
813}
814
815
816int ZCone::ambientDimension()const
817{
818  return n;
819}
820
821
822int ZCone::codimension()const
823{
824  return ambientDimension()-dimension();
825}
826
827
828int ZCone::dimension()const
829{
830//  assert(state>=1);
831  ensureStateAsMinimum(1);
832  return ambientDimension()-equations.getHeight();
833}
834
835
836int ZCone::dimensionOfLinealitySpace()const
837{
838  ZMatrix temp=inequalities;
839  temp.append(equations);
840  ZCone temp2(ZMatrix(0,n),temp);
841  return temp2.dimension();
842}
843
844
845bool ZCone::isOrigin()const
846{
847  return dimension()==0;
848}
849
850
851bool ZCone::isFullSpace()const
852{
853  for(int i=0;i<inequalities.getHeight();i++)
854    if(!inequalities[i].isZero())return false;
855  for(int i=0;i<equations.getHeight();i++)
856    if(!equations[i].isZero())return false;
857  return true;
858}
859
860
861ZCone intersection(const ZCone &a, const ZCone &b)
862{
863  assert(a.ambientDimension()==b.ambientDimension());
864  ZMatrix inequalities=a.inequalities;
865  inequalities.append(b.inequalities);
866  ZMatrix equations=a.equations;
867  equations.append(b.equations);
868
869  equations.sortAndRemoveDuplicateRows();
870  inequalities.sortAndRemoveDuplicateRows();
871
872  {
873    ZMatrix Aequations=a.equations;
874    ZMatrix Ainequalities=a.inequalities;
875    Aequations.sortAndRemoveDuplicateRows();
876    Ainequalities.sortAndRemoveDuplicateRows();
877    if((Ainequalities.getHeight()==inequalities.getHeight()) && (Aequations.getHeight()==equations.getHeight()))return a;
878    ZMatrix Bequations=b.equations;
879    ZMatrix Binequalities=b.inequalities;
880    Bequations.sortAndRemoveDuplicateRows();
881    Binequalities.sortAndRemoveDuplicateRows();
882    if((Binequalities.getHeight()==inequalities.getHeight()) && (Bequations.getHeight()==equations.getHeight()))return b;
883  }
884
885  return ZCone(inequalities,equations);
886}
887
888/*
889PolyhedralCone product(const PolyhedralCone &a, const PolyhedralCone &b)
890{
891  IntegerVectorList equations2;
892  IntegerVectorList inequalities2;
893
894  int n1=a.n;
895  int n2=b.n;
896
897  for(IntegerVectorList::const_iterator i=a.equations.begin();i!=a.equations.end();i++)
898    equations2.push_back(concatenation(*i,IntegerVector(n2)));
899  for(IntegerVectorList::const_iterator i=b.equations.begin();i!=b.equations.end();i++)
900    equations2.push_back(concatenation(IntegerVector(n1),*i));
901  for(IntegerVectorList::const_iterator i=a.inequalities.begin();i!=a.inequalities.end();i++)
902    inequalities2.push_back(concatenation(*i,IntegerVector(n2)));
903  for(IntegerVectorList::const_iterator i=b.inequalities.begin();i!=b.inequalities.end();i++)
904    inequalities2.push_back(concatenation(IntegerVector(n1),*i));
905
906  PolyhedralCone ret(inequalities2,equations2,n1+n2);
907  ret.setMultiplicity(a.getMultiplicity()*b.getMultiplicity());
908  ret.setLinearForm(concatenation(a.getLinearForm(),b.getLinearForm()));
909
910  ret.ensureStateAsMinimum(a.state);
911  ret.ensureStateAsMinimum(b.state);
912
913  return ret;
914}*/
915
916
917ZCone ZCone::positiveOrthant(int dimension)
918{
919  return ZCone(ZMatrix::identity(dimension),ZMatrix(0,dimension));
920}
921
922
923ZCone ZCone::givenByRays(ZMatrix const &generators, ZMatrix const &linealitySpace)
924{
925  ZCone dual(generators,linealitySpace);
926  ZMatrix inequalities=dual.extremeRays();
927  ZMatrix equations=dual.generatorsOfLinealitySpace();
928
929  return ZCone(inequalities,equations,3);
930}
931
932
933bool ZCone::containsPositiveVector()const
934{
935  ZCone temp=intersection(*this,ZCone::positiveOrthant(n));
936  return temp.getRelativeInteriorPoint().isPositive();
937}
938
939
940bool ZCone::contains(ZVector const &v)const
941{
942  for(int i=0;i<equations.getHeight();i++)
943    {
944      if(!dot(equations[i],v).isZero())return false;
945    }
946  for(int i=0;i<inequalities.getHeight();i++)
947    {
948      if(dot(inequalities[i],v).sign()<0)return false;
949    }
950  return true;
951}
952
953
954bool ZCone::containsRowsOf(ZMatrix const &m)const
955{
956  for(int i=0;i<m.getHeight();i++)
957    if(!contains(m[i]))return false;
958  return true;
959}
960
961
962bool ZCone::contains(ZCone const &c)const
963{
964  ZCone c2=intersection(*this,c);
965  ZCone c3=c;
966  c2.canonicalize();
967  c3.canonicalize();
968  return !(c2!=c3);
969}
970
971
972bool ZCone::containsRelatively(ZVector const &v)const
973{
974  ensureStateAsMinimum(1);
975//  assert(state>=1);
976  for(int i=0;i<equations.getHeight();i++)
977    {
978      if(!dot(equations[i],v).isZero())return false;
979    }
980  for(int i=0;i<inequalities.getHeight();i++)
981    {
982      if(dot(inequalities[i],v).sign()<=0)return false;
983    }
984  return true;
985}
986
987
988bool ZCone::isSimplicial()const
989{
990//  assert(state>=2);
991  ensureStateAsMinimum(2);
992  return codimension()+inequalities.getHeight()+dimensionOfLinealitySpace()==n;
993}
994
995
996ZCone ZCone::linealitySpace()const
997{
998  ZCone ret(ZMatrix(0,n),combineOnTop(equations,inequalities));
999//  ret.ensureStateAsMinimum(state);
1000  return ret;
1001}
1002
1003
1004ZCone ZCone::dualCone()const
1005{
1006  ensureStateAsMinimum(1);
1007//  assert(state>=1);
1008
1009  ZMatrix dualInequalities,dualEquations;
1010  lpSolver.dual(inequalities,equations,dualInequalities,dualEquations);
1011  ZCone ret(dualInequalities,dualEquations);
1012  ret.ensureStateAsMinimum(state);
1013
1014  return ret;
1015}
1016
1017
1018ZCone ZCone::negated()const
1019{
1020  ZCone ret(-inequalities,equations,(areFacetsKnown()?PCP_facetsKnown:0)|(areImpliedEquationsKnown()?PCP_impliedEquationsKnown:0));
1021//  ret.ensureStateAsMinimum(state);
1022  return ret;
1023}
1024
1025
1026ZMatrix ZCone::extremeRays(ZMatrix const *generatorsOfLinealitySpace)const
1027{
1028//  assert((dimension()==ambientDimension()) || (state>=3));
1029  if(dimension()!=ambientDimension())
1030    ensureStateAsMinimum(3);
1031
1032  if(haveExtremeRaysBeenCached)return cachedExtremeRays;
1033  ZMatrix ret(0,n);
1034  std::vector<std::vector<int> > indices=lpSolver.extremeRaysInequalityIndices(inequalities);
1035
1036  for(int i=0;i<indices.size();i++)
1037    {
1038      /* At this point we know lineality space, implied equations and
1039         also inequalities for the ray. To construct a vector on the
1040         ray which is stable under (or indendent of) angle and
1041         linarity preserving transformation we find the dimension 1
1042         subspace orthorgonal to the implied equations and the
1043         lineality space and pick a suitable primitive generator */
1044
1045          /* To be more precise,
1046           * let E be the set of equations, and v the inequality defining a  ray R.
1047           * We wish to find a vector satisfying these, but it must also be orthogonal
1048           * to the lineality space of the cone, that is, in the span of {E,v}.
1049           * One way to get such a vector is to project v to E an get a vector p.
1050           * Then v-p is in the span of {E,v} by construction.
1051           * The vector v-p is also in the orthogonal complement to E by construction,
1052           * that is, the span of R.
1053           * We wish to argue that it is not zero.
1054           * That would imply that v=p, meaning that v is in the span of the equations.
1055           * However, that would contradict that R is a ray.
1056           * In case v-p does not satisfy the inequality v (is this possible?), we change the sign.
1057           *
1058           * As a consequence we need the following procedure
1059           * primitiveProjection():
1060           *    Input: E,v
1061           *    Output: A primitive representation of the vector v-p, where p is the projection of v onto E
1062           *
1063           * Notice that the output is a Q linear combination of the input and that p is
1064           * a linear combination of E. The check that p has been computed correctly,
1065           * it suffices to check that v-p satisfies the equations E.
1066           * The routine will actually first compute a multiple of v-p.
1067           * It will do this using floating point arithmetics. It will then transform
1068           * the coefficients to get the multiple of v-p into integers. Then it
1069           * verifies in exact arithmetics, that with these coefficients we get a point
1070           * satisfying E. It then returns the primitive vector on the ray v-p.
1071           * In case of a failure it falls back to an implementation using rational arithmetics.
1072           */
1073
1074
1075          std::vector<int> asVector(inequalities.getHeight());
1076          for(int j=0;j<indices[i].size();j++){asVector[indices[i][j]]=1;}
1077          ZMatrix equations=this->equations;
1078          ZVector theInequality;
1079
1080          for(int j=0;j<asVector.size();j++)
1081            if(asVector[j])
1082              equations.appendRow(inequalities[j]);
1083            else
1084              theInequality=inequalities[j];
1085
1086          assert(!theInequality.isZero());
1087
1088          ZVector thePrimitiveVector;
1089          if(generatorsOfLinealitySpace)
1090          {
1091            QMatrix temp=ZToQMatrix(combineOnTop(equations,*generatorsOfLinealitySpace));
1092            thePrimitiveVector=QToZVectorPrimitive(temp.reduceAndComputeVectorInKernel());
1093          }
1094          else
1095          {
1096            QMatrix linealitySpaceOrth=ZToQMatrix(combineOnTop(this->equations,inequalities));
1097
1098
1099            QMatrix temp=combineOnTop(linealitySpaceOrth.reduceAndComputeKernel(),ZToQMatrix(equations));
1100            thePrimitiveVector=QToZVectorPrimitive(temp.reduceAndComputeVectorInKernel());
1101          }
1102          if(!contains(thePrimitiveVector))thePrimitiveVector=-thePrimitiveVector;
1103          ret.appendRow(thePrimitiveVector);
1104    }
1105
1106  cachedExtremeRays=ret;
1107  haveExtremeRaysBeenCached=true;
1108
1109  return ret;
1110}
1111
1112
1113Integer ZCone::getMultiplicity()const
1114{
1115  return multiplicity;
1116}
1117
1118
1119void ZCone::setMultiplicity(Integer const &m)
1120{
1121  multiplicity=m;
1122}
1123
1124
1125ZMatrix ZCone::getLinearForms()const
1126{
1127  return linearForms;
1128}
1129
1130
1131void ZCone::setLinearForms(ZMatrix const &linearForms_)
1132{
1133  linearForms=linearForms_;
1134}
1135
1136
1137ZMatrix ZCone::quotientLatticeBasis()const
1138{
1139//  assert(isInStateMinimum(1));// Implied equations must have been computed in order to know the span of the cone
1140  ensureStateAsMinimum(1);
1141
1142
1143  int a=equations.getHeight();
1144  int b=inequalities.getHeight();
1145
1146  // Implementation below could be moved to nonLP part of code.
1147
1148  // small vector space defined by a+b equations.... big by a equations.
1149
1150  ZMatrix M=combineLeftRight(combineLeftRight(
1151                                                  equations.transposed(),
1152                                                  inequalities.transposed()),
1153                                 ZMatrix::identity(n));
1154  M.reduce(false,true);
1155  /*
1156    [A|B|I] is reduced to [A'|B'|C'] meaning [A'|B']=C'[A|B] and A'=C'A.
1157
1158    [A'|B'] is in row echelon form, implying that the rows of C' corresponding to zero rows
1159    of [A'|B'] generate the lattice cokernel of [A|B] - that is the linealityspace intersected with Z^n.
1160
1161    [A'] is in row echelon form, implying that the rows of C' corresponding to zero rows of [A'] generate
1162    the lattice cokernel of [A] - that is the span of the cone intersected with Z^n.
1163
1164    It is clear that the second row set is a superset of the first. Their difference is a basis for the quotient.
1165   */
1166  ZMatrix ret(0,n);
1167
1168  for(int i=0;i<M.getHeight();i++)
1169    if(M[i].subvector(0,a).isZero()&&!M[i].subvector(a,a+b).isZero())
1170      {
1171        ret.appendRow(M[i].subvector(a+b,a+b+n));
1172      }
1173
1174  return ret;
1175}
1176
1177
1178ZVector ZCone::semiGroupGeneratorOfRay()const
1179{
1180  ZMatrix temp=quotientLatticeBasis();
1181  assert(temp.getHeight()==1);
1182  for(int i=0;i<inequalities.getHeight();i++)
1183    if(dot(temp[0],inequalities[i]).sign()<0)
1184      {
1185        temp[0]=-temp[0];
1186        break;
1187      }
1188  return temp[0];
1189}
1190
1191
1192ZCone ZCone::link(ZVector const &w)const
1193{
1194  /* Observe that the inequalities giving rise to facets
1195   * also give facets in the link, if they are kept as
1196   * inequalities. This means that the state cannot decrease
1197   * when taking links - that is why we specify the PCP flags.
1198   */
1199  ZMatrix inequalities2(0,n);
1200  for(int j=0;j<inequalities.getHeight();j++)
1201    if(dot(w,inequalities[j]).sign()==0)inequalities2.appendRow(inequalities[j]);
1202  ZCone C(inequalities2,equations,(areImpliedEquationsKnown()?PCP_impliedEquationsKnown:0)|(areFacetsKnown()?PCP_facetsKnown:0));
1203  C.ensureStateAsMinimum(state);
1204
1205  C.setLinearForms(getLinearForms());
1206  C.setMultiplicity(getMultiplicity());
1207
1208  return C;
1209}
1210
1211bool ZCone::hasFace(ZCone const &f)const
1212{
1213  if(!contains(f.getRelativeInteriorPoint()))return false;
1214  ZCone temp1=faceContaining(f.getRelativeInteriorPoint());
1215  temp1.canonicalize();
1216  ZCone temp2=f;
1217  temp2.canonicalize();
1218
1219  return !(temp2!=temp1);
1220}
1221
1222ZCone ZCone::faceContaining(ZVector const &v)const
1223{
1224  assert(n==v.size());
1225  assert(contains(v));
1226  ZMatrix newEquations=equations;
1227  ZMatrix newInequalities(0,n);
1228  for(int i=0;i<inequalities.getHeight();i++)
1229    if(dot(inequalities[i],v).sign()!=0)
1230      newInequalities.appendRow(inequalities[i]);
1231    else
1232      newEquations.appendRow(inequalities[i]);
1233
1234  ZCone ret(newInequalities,newEquations,(state>=1)?PCP_impliedEquationsKnown:0);
1235  ret.ensureStateAsMinimum(state);
1236  return ret;
1237}
1238
1239
1240ZMatrix ZCone::getInequalities()const
1241{
1242  return inequalities;
1243}
1244
1245
1246ZMatrix ZCone::getEquations()const
1247{
1248  return equations;
1249}
1250
1251
1252ZMatrix ZCone::generatorsOfSpan()const
1253{
1254  ensureStateAsMinimum(1);
1255  QMatrix l=ZToQMatrix(equations);
1256  return QToZMatrixPrimitive(l.reduceAndComputeKernel());
1257}
1258
1259
1260ZMatrix ZCone::generatorsOfLinealitySpace()const
1261{
1262  QMatrix l=ZToQMatrix(combineOnTop(inequalities,equations));
1263  return QToZMatrixPrimitive(l.reduceAndComputeKernel());
1264}
1265
1266};
Note: See TracBrowser for help on using the repository browser.