~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

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
#include "traverser_stableintersection.h"
#include <iostream>
#include "halfopencone.h"
#include "tropical.h"
#include "division.h"
#include "wallideal.h"
#include "groebnerengine.h"
#include "tropical2.h"
#include "multiplicity.h"
#include "mixedvolume.h"
#include "log.h"

StableIntersectionTraverser::StableIntersectionTraverser(PolynomialSet const &coneGroebnerBasis_, PolynomialSet const &idealGroebnerBasis_):
	coneGroebnerBasis(coneGroebnerBasis_),
	idealGroebnerBasis(idealGroebnerBasis_),
	theCone(coneGroebnerBasis_.getRing().getNumberOfVariables())
	{
		n=coneGroebnerBasis_.getRing().getNumberOfVariables();
		updatePolyhedralCone(mixedVolume(coneGroebnerBasis));
		PolyhedralCone homogeneitySpac=homogeneitySpace(coneGroebnerBasis);
		d=homogeneitySpac.dimensionOfLinealitySpace();
	}

void StableIntersectionTraverser::updatePolyhedralCone(int multiplicity)
{
	theCone=PolyhedralCone(fastNormals(wallInequalities(idealGroebnerBasis)),wallInequalities(coneGroebnerBasis),n);
	theCone.canonicalize();
	theCone.setMultiplicity(multiplicity);
}

static void checkSameLeadingTerms(PolynomialSet const &a, PolynomialSet const &b)
{
  assert(a.size()==b.size());
  PolynomialSet::const_iterator A=a.begin();
  for(PolynomialSet::const_iterator B=b.begin();B!=b.end();B++,A++)
    assert(A->getMarked().m.exponent==B->getMarked().m.exponent);
}

void StableIntersectionTraverser::changeCone(IntegerVector const &ridgeVector, IntegerVector const &rayVector)
{
	IntegerVectorList m;
	m.push_back(ridgeVector);
	m.push_back(rayVector);
	MatrixTermOrder T(m);

	idealGroebnerBasis.markAndScale(T);
	coneGroebnerBasis=initialFormsAssumeMarked(initialFormsAssumeMarked(idealGroebnerBasis,ridgeVector),rayVector);

	updatePolyhedralCone(mixedVolume(coneGroebnerBasis));
	checkSameLeadingTerms(idealGroebnerBasis,coneGroebnerBasis);
}

IntegerVectorList StableIntersectionTraverser::link(IntegerVector const &ridgeVector)
{
	assert(idealGroebnerBasis.containsInClosedGroebnerCone(ridgeVector));
	PolynomialSet tempIdeal=initialFormsAssumeMarked(idealGroebnerBasis,ridgeVector);

	tempIdeal.simplestPolynomialsFirst();

	PolyhedralFan theLink=tropicalHyperSurfaceIntersectionClosed(n,tempIdeal);

	IntegerVectorList rays1=theLink.getRaysInPrintingOrder(0);

	log2 {
	  cerr<<"Ray candidates:"<<endl;
	  AsciiPrinter(Stderr)<<rays1;
	}

	IntegerVectorList rays2;
	for(IntegerVectorList::const_iterator i=rays1.begin();i!=rays1.end();i++)
		if(mixedVolumePositive(initialForms(tempIdeal,*i)))
			{
			PolyhedralCone theRidge=theCone.faceContaining(ridgeVector);
                        IntegerVectorList linealityGenerators=theRidge.generatorsOfLinealitySpace();
//                        IntegerVectorList linealityGenerators=theRidge.linealitySpace().dualCone().getEquations();
			IntegerVectorList generators;generators.push_back(*i);
			PolyhedralCone ray=PolyhedralCone::givenByRays(generators,linealityGenerators,n);
			ray.canonicalize();
			rays2.push_back(ray.getUniquePoint());
			//			rays2.push_back(*i);
			}
	log2{
	cerr<<"True rays:"<<endl;
	AsciiPrinter(Stderr)<<rays2;
	}

	return rays2;
}

PolyhedralCone & StableIntersectionTraverser::refToPolyhedralCone()
{
	return theCone;
}