//
// $Source: /cvsroot/gambit/gambit/sources/tools/enumpoly/pelclass.cc,v $
// $Date: 2006/08/17 03:05:00 $
// $Revision: 1.7 $
//
// DESCRIPTION:
// Implementation of interface to Pelican
//
// This file is part of Gambit
// Copyright (c) 2002, The Gambit Project
//
// 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 of the License, or
// (at your option) any later version.
//
// 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.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//
#include "pelclass.h"
/*
#include "glist.h"
#include "gvector.h"
#include "complex.h"
*/
/*************************************************************/
/************** Implementation of class Pelview **************/
/*************************************************************/
node SaveList=0;
void PelView::InitializePelicanMemory() const
{
// mimicking main in Shell.c
LOCS(1);
init_symbol_table();
node_init_store();
PUSH_LOC(SaveList);
SaveList=node_new();
}
Pring PelView::MakePring(const int num) const
{
char* q[] = {"","n1","n2" , "n3", "n4", "n5", "n6", "n7", "n8",
"n9", "n10", "n11", "n12", "n13", "n14", "n15", "n16",
"n17", "n18","n19", "n20", "n21", "n22", "n23", "n24",
"n25", "n26", "n27", "n28","n29", "n30"};
Pring R = makePR(num); // makePR is a memory malloc routine
for(int j=1; j<=R->n; j++)
R->vars[j-1] = q[j];
R->def = "t";
R->n = num;
return R;
}
void PelView::PrintPring(const Pring &ring) const
{
// gout << "The Pring has " << ring->n << " variables: ";
for (int i = 1; i <= ring->n; i++) {
// gout << ring->vars[i-1];
//if (i < ring->n)
// gout << ", ";
//else
// gout << ".\n";
}
// gout << " The homotopy variable is " << ring->def << ".\n";
}
void PelView::Initialize_Idf_T_Gen_node(const Gen_node &node,
const char * label) const
{
node->type=Idf_T;
node->Genval.gval=Copy_String_NQ((char *)label);
node->Genval.idval=Copy_String_NQ((char *)label);
}
Gen_node PelView::CreateRing(const int numvar) const
{
char* q[] = {" ", "","n1","n2" , "n3", "n4", "n5", "n6", "n7", "n8",
"n9", "n10", "n11", "n12", "n13", "n14", "n15", "n16",
"n17", "n18","n19", "n20", "n21", "n22", "n23", "n24",
"n25", "n26", "n27", "n28","n29", "n30"};
Gen_node a1 = gen_node();
Initialize_Idf_T_Gen_node(a1,"n1");
Gen_node atemp;
atemp = gen_node();
if (numvar == 1) a1->next= atemp;
else {
a1->next=NULL;
atemp = a1;
for(int j=2; j<=numvar; j++)
{
Gen_node a = gen_node();
Initialize_Idf_T_Gen_node(a,q[j]);
a->next = NULL;
atemp->next = a;
atemp = a;
}
}
atemp->next=NULL;
Initialize_Idf_T_Gen_node(atemp,"t");
return a1;
}
polynomial1 PelView::GamPolyToPelPoly(const gPoly<double> &p,
const int n,
const Pring ring) const
{
if ((p.MonomialList().Length() == 0)) {
polynomial1 P;
P = makeP(ring);
P->coef.r=0;
P->coef.i=0;
for (int i=0;i<P->R->n;i++)
P->exps[i]=0;
P->next= 0;
return P;
}
else {
polynomial1 P,Ptemp;
P = makeP(ring);
Ptemp= P;
if (p.MonomialList()[1].IsConstant()) {
Ptemp->coef.r= p.MonomialList()[1].Coef();
Ptemp->coef.i= 0;
for (int i=0;i<Ptemp->R->n;i++)
Ptemp->exps[i]=0;
Ptemp->next=0;
}
else {
Ptemp->coef.r= p.MonomialList()[1].Coef();
Ptemp->coef.i= 0;
for (int i=0;i<Ptemp->R->n;i++)
Ptemp->exps[i]= p.MonomialList()[1].ExpV() [i+1];
Ptemp->next=0;
}
for (int j = 2; j <=p.MonomialList().Length(); j++) {
polynomial1 a;
a = makeP(ring);
if (p.MonomialList()[j].IsConstant()) {
a->coef.r= p.MonomialList()[j].Coef();
a->coef.i= 0;
for (int i=0;i<a->R->n;i++)
a->exps[i]=0;
a->next = 0;
Ptemp->next = a;
Ptemp = a;
}
else {
a->coef.r= p.MonomialList()[j].Coef();
a->coef.i= 0;
for (int i=0;i<a->R->n;i++)
a->exps[i]= p.MonomialList()[j].ExpV() [i+1];
// printP(a);
// cout;
a->next = 0;
Ptemp->next = a;
Ptemp = a;
}
}
return P;
}
}
Gen_node
PelView::CreatePelicanVersionOfSystem(const gPolyList<double> &input,
const Pring ring) const
{
Gen_node a;
a = gen_node();
a->type= Mtx_T;
a->next= NULL;
Gmatrix V;
V = Gmatrix_new(1, input.Length());
V->store = input.Length();
V->topc = input.Length();
V->ncols = input.Length();
V->topr=1;
for(int j=1; j<=input.Length(); j++)
{
Gen_node b;
b = gen_node();
b->type = Ply_T;
b->next = NULL;
b->Genval.pval=(GamPolyToPelPoly(input[j], input.Length(), ring));
b->Genval.gval=(char*) (GamPolyToPelPoly(input[j], input.Length(), ring));
b->Genval.idval=(char*) (GamPolyToPelPoly(input[j], input.Length(), ring));
(V->coords[j-1])= b;
}
a->Genval.gval=(char *)V;
return G_System(a);
}
int PelView::GutsOfGetMixedVolume( node A,
node norms,
const Imatrix T) const
{
node ptr = 0, ptc = 0, res = 0;
Imatrix M = 0, Tp = 0;
int v, mv = 0, t = 0;
LOCS(5);
PUSH_LOC(A);
PUSH_LOC(res);
PUSH_LOC(ptc);
PUSH_LOC(ptr);
PUSH_LOC(norms);
ptr = norms;
while (ptr != 0) {
ptc = aset_face(A, (Imatrix) Car((node) Car(ptr)));
Tp = aset_type(ptc, Tp);
if (T != 0 && Imatrix_equal(Tp, T) == TRUE) {
t = 1;
list_insert(Car(ptr),&res, &(list_Imatrix_comp),FALSE);
} else
t = 0;
M = aset_M(ptc, M);
if (
(Imatrix_rref(M, &v) == aset_dim(A) - 1) &&
(IMrows(M) == aset_dim(A) - 1)
) {
if (t == 1)
mv += abs(v);
}
ptr = Cdr(ptr);
}
Imatrix_free(Tp);
Imatrix_free(M);
POP_LOCS();
return mv;
}
int PelView::GetMixedVolume(const Gen_node g) const
{
aset A=0;
Ivector T=0;
int r = 0;
int CP;
int nargs;
LOCS(2);
PUSH_LOC(A);
nargs=Gen_length(g);
//Something happens right here and is necessary for the code,
//but don't ask me what...
if ((nargs ==0) ||
(nargs==1 && Can_Be_Aset(Gen_elt(g,1))!=TRUE )||
(nargs==2 && (r=Can_Be_Vector(Gen_elt(g,2),Int_T))<0)){
}
A=Gen_aset(Gen_elt(g,1));
if (nargs==2 )
if (r==aset_r(A))
T=Gen_to_Imatrix(Gen_elt(g,2));
CP= GutsOfGetMixedVolume(A,aset_lower_facets(A),T);
return CP;
}
Gen_node PelView::Make_scl_Gen_node() const
{
Gen_node g= gen_node();
g->type= Idf_T;
g->next= NULL;
g->Genval.gval= "scl";
g->Genval.idval= "scl";
return g;
}
Gen_node PelView::ToDmatrixGen_node(const Gen_node g) const
{
Gen_node a;
psys PS;
Dmatrix S;
PS=Gen_to_psys(g);
S=psys_scale(PS);
a= Dmatrix_to_Gen(S);
return a;
}
polynomial1 PelView::IdentityElementPoly(const Pring ring) const
{
polynomial1 P;
P = makeP(ring);
P->coef.r=1;
P->coef.i=0;
for (int i=0;i<P->R->n;i++)
P->exps[i]=0;
P->def =0;
P->next= 0;
return P;
}
polynomial1 PelView::HomotopyVariableMonomialPoly(const Pring ring,
const int comp) const
{
polynomial1 P;
P = makeP(ring);
P->coef.r=1;
P->coef.i=0;
P->remaining=0;
P->homog=0;
for (int i=0;i<P->R->n;i++)
P->exps[i]=0;
P->def = comp;
P->next= 0;
return P;
}
Gen_node PelView::SolutionsDerivedFromContinuation(const Pring &ring,
const Gen_node &Genpoly,
const Gen_node &Solve,
const Gen_node &pel_system,
int tweak) const
{
pel_system->next = 0;
Gen_node temp= Link(pel_system ,Make_scl_Gen_node());
Gen_node scl= ToDmatrixGen_node(temp);
Gen_node fs = G_Scale(temp);
// gout<< " Step 12 \n";
// polynomials used for the homothopy
int unity = 1;
Gen_node tee, one, tee1;
tee = Ply_To_Gen(HomotopyVariableMonomialPoly(ring, unity));
one=Ply_To_Gen(IdentityElementPoly(ring));
tee1=Ply_To_Gen(HomotopyVariableMonomialPoly(ring, unity));
tee->next=0;
one->next=0;
// Define the homothopy
// gout<< " Step 13 \n";
Gen_node h = PROC_MUL(Link(tee, fs));
Gen_node h1 = PROC_SUB(Link(one, tee1));
//A little trick
Genpoly->next=0;
Gen_node h2 = G_UnLift(Genpoly);
// gout<< " Step 14 \n";
Gen_node h4 = PROC_MUL(Link(h1, h2));
Gen_node h5 = PROC_ADD(Link(h,h4));
Gen_node cs = G_Cont(Link(h5,G_UnLift(Solve)), tweak);
// gout<< " Step 15 \n";
Gen_node sols = G_Affine(G_UnScale(Link(cs, scl)));
free_Gen_list(scl);
return sols;
}
template <class T> bool HasARedundancy(const Gambit::List<T> &p_list)
{
int i = 1; int j = 2;
while (i < p_list.Length()) {
if (p_list[i] == p_list[j])
return true;
else
j++;
if (j > p_list.Length()) { i++; j = i+1; }
}
return false;
}
//#define PELVIEW_DEBUG
Gen_node PelView::SolveCheckMaybeTryAgain(const Pring &ring,
const Gen_node &Genpoly,
const Gen_node &Qtrig,
const Gen_node &pel_system)
{
Gen_node sols, Solve;
bool done(false);
for (int tweak = 0; tweak <= 3 && !done; tweak ++)
{
Solve = G_Solve(Link(Genpoly, Qtrig), tweak);
#ifdef PELVIEW_DEBUG
print_Gen_list(Solve);
// gout<< " \n Step 11 \n";
#endif
sols = SolutionsDerivedFromContinuation(ring,
Genpoly,
Solve,
pel_system,
tweak);
#ifdef PELVIEW_DEBUG
// gout<<" \n The solution list produced by Pelican is:\n";
print_Gen_list(sols);
#endif
//The same little trick
pel_system->next= 0;
#ifdef PELVIEW_DEBUG
// gout << "\n Step 12 \n";
#endif
solutionsarecorrect = CheckSolutions(G_Verify(Link(pel_system,sols)));
if (solutionsarecorrect &&
!HasARedundancy(GambitRootsFromPelRoots(sols)))
done = true;
if (!done) free_Gen_list(sols);
#ifdef PELVIEW_DEBUG
// gout << "The solutions are ";
// if (!solutionsarecorrect)
// gout << "not ";
// gout << "correct.\n";
#endif
}
return sols;
// gout<< "\n Solve \n"; print_Gen_list(Solve);
// free_Gen_list(Solve);
// Solve doesn't seem to exist at this point, yet this is very ambiguous.
}
Gambit::List<Gambit::Vector<gComplex> >
PelView::GambitRootsFromPelRoots(const Gen_node g) const
{
Gambit::List<Gambit::Vector<gComplex> > alist;
node ptr;
ptr = Gen_to_Dvector_list(Gen_lval(copy_Gen_node(g)));
Dmatrix P;
int i;
while(ptr!=0) {
P=(Dmatrix)(Car(Car(ptr)));
int numbervar;
numbervar= (int)(DVlength(P)-3)/2;
Gambit::Vector<gComplex> vector(1, numbervar);
int j = 1;
for(i=3; i<DVlength(P);i++) {
if (i%2==0 && i<DVlength(P)) {
double re,im;
re = DVref(P,i-1);
im = DVref(P,i);
gComplex complexsol(re, im);
vector[j] = complexsol;
j+=1;
}
}
alist.Append(vector);
ptr=Cdr(ptr);
}
return alist;
}
void PelView::DisplayComplexRootList(const Gambit::List<Gambit::Vector<gComplex> >
&complexroots) const
{
#ifdef UNUSED
for (int k = 1; k <= complexroots.Length(); k++)
for (int m = 1; m <= complexroots[k].Length(); m++) {
if (m ==1)
gout << "{";
for (int n = 1; n <= m; n++)
gout << " ";
gout << complexroots[k][m];
if (m < complexroots[k].Length())
gout << ",\n";
else
gout << " }\n";
}
#endif // UNUSED
}
int PelView::Dmnsn() const
{
return input.Dmnsn();
}
Gambit::List<Gambit::Vector<double> >
PelView::RealRoots(const Gambit::List<Gambit::Vector<gComplex> > &clist) const
{
Gambit::List<Gambit::Vector<double> > answer;
for (int i = 1; i <= clist.Length(); i++) {
bool is_real = true;
for (int j = 1; j <= Dmnsn(); j++)
if (Gambit::abs(clist[i][j].ImaginaryPart()) > 0.0001)
is_real = false;
if (is_real) {
Gambit::Vector<double> next(Dmnsn());
for (int j = 1; j <= Dmnsn(); j++)
next[j] = (double)clist[i][j].RealPart();
answer.Append(next);
}
}
return answer;
}
bool PelView::CheckSolutions(const Gen_node g) const
{
Gen_node goo;
goo = g->Genval.lval;
while (goo!=0) {
if (Gambit::abs(goo->Genval.dval) > 0.01)
return 0;
goo = goo->next;
}
return 1;
}
PelView::PelView(const gPolyList<double> &mylist):input(mylist)
{
InitializePelicanMemory();
#ifdef PELVIEW_DEBUG
// gout << "We begin with the polynomial list\n" << mylist << "\n";
#endif
Pring ring = MakePring(input.Length());
#ifdef PELVIEW_DEBUG
// gout<< " Step 1 \n";
PrintPring(ring);
#endif
Gen_node vic = Set_Ring((CreateRing(input.Length()+1)));
#ifdef PELVIEW_DEBUG
// gout<< "\n Step 2 \n";
print_Gen_node(vic);
#endif
// Gen_node pel_system = G_System(translate(input,ring));
Gen_node pel_system = CreatePelicanVersionOfSystem(input,ring);
#ifdef PELVIEW_DEBUG
// gout<< "\n\n Step 3 \n";
// gout << "The translated system is:\n";
print_Gen_node(pel_system);
#endif
Gen_node Atype= G_AType(pel_system); // Atype is the vector of numbers of
#ifdef PELVIEW_DEBUG
// gout<< " Step 4 \n";
print_Gen_list(Atype); // polys with each support type
#endif
Gen_node Aset = G_Aset(pel_system);
#ifdef PELVIEW_DEBUG
// gout<< "\n Step 5 \n";
// gout << "The list of unlifted vertex tuples is:\n";
print_Gen_node(Aset);
#endif
Gen_node Randlift = G_RandLift(Aset);
#ifdef PELVIEW_DEBUG
// gout<< "\n\n Step 6 \n";
// gout<< "After the random lift, the vertex tuples are:\n";
print_Gen_list(Randlift);
// gout << "\n\n Step 7 \n";
// gout << "To see the cells in the subdivision, define ACTUALLY_PRINT in pelclqhl.cc:\n";
#endif
Gen_node Qtrig = G_Qtrig(Link(Randlift, Atype));
/* - There is a commented out error in node_push_local that was triggered
by this in various conditions. In particular, commented out the final print
avoids the error, somehow!!
// gout<< "\n Step 8 \n";
// gout<< "Before Randlift is ";
// gout << Randlift << "\n";
silent_print_Gen_list(Randlift);
// gout<< "After Randlift is " << Randlift << "\n";
*/
mixedvolume = GetMixedVolume(Randlift);
#ifdef PELVIEW_DEBUG
// gout<< "\n Step 8 \n";
// gout << "The mixed volume is " << mixedvolume << ".\n";
#endif
Gen_node Genpoly = G_Gen_Poly(Randlift);
#ifdef PELVIEW_DEBUG
// gout<< "\n Step 9 \n";
// gout << "The homotopy system is:\n";
print_Gen_list(Genpoly);
// gout<< "\n Step 10 \n";
// gout << "The Gen_node Solve (solutions of start system) is: \n";
#endif
Gen_node sols = SolveCheckMaybeTryAgain(ring, Genpoly, Qtrig, pel_system);
complexroots = GambitRootsFromPelRoots(sols);
#ifdef PELVIEW_DEBUG
// gout << "\n Step 13 \n";
// gout<<" After conversion to Gambit, the complex roots are...\n";
DisplayComplexRootList(complexroots);
#endif
realroots = RealRoots(complexroots);
#ifdef PELVIEW_DEBUG
// gout << "\n Step 14 \n";
// gout<<" The real solutions are...\n";
// for (int k = 1; k <= realroots.Length(); k++)
// gout << realroots[k] << "\n";
// gout<< "\n Step 15 \n";
// gout<< "\n Memory testing \n";
print_Gen_list(Genpoly);
#endif
free_Gen_list(Genpoly);
#ifdef PELVIEW_DEBUG
// gout<< "\n Qtrig \n"; print_Gen_list(Qtrig);
#endif
free_Gen_list(Qtrig);
free_Gen_list(vic);
free_Gen_list(sols);
#ifdef PELVIEW_DEBUG
// gout << "And the memory should be clean...\n";
#endif
}
PelView::PelView(const PelView & given)
: input(given.input),
complexroots(given.complexroots),
realroots(given.realroots),
mixedvolume(given.mixedvolume),
solutionsarecorrect(given.solutionsarecorrect)
{
}
PelView::~PelView()
{
}
PelView& PelView::operator =(const PelView &rhs)
{
// gout << "For (eventual) const'ness, operator = not allowed for PelView\n";
exit (0);
return *this;
}
bool PelView::operator ==(const PelView &rhs) const
{
return (input == rhs.input &&
complexroots == rhs.complexroots &&
realroots == rhs.realroots &&
mixedvolume == rhs.mixedvolume &&
solutionsarecorrect == rhs.solutionsarecorrect);
}
bool PelView::operator !=(const PelView &rhs) const
{
return !(*this == rhs);
}
Gambit::List<Gambit::Vector<gComplex> > PelView::ComplexRoots() const
{
return complexroots;
}
Gambit::List<Gambit::Vector<double> > PelView::RealRoots() const
{
return realroots;
}
int PelView::MixedVolume() const
{
return mixedvolume;
}
int PelView::NumComplexRoots() const
{
return complexroots.Length();
}
bool PelView::FoundAllRoots() const
{
return (NumComplexRoots() == MixedVolume());
}
int PelView::NumRealRoots() const
{
return realroots.Length();
}
/**********************************************************************/
// was the original main in the pelican driver
int old_main()
{
//stuff for Gambit-text to Gambit-data
for (int loop = 1; loop <= 3; loop++) {
gSpace Space(4);
ORD_PTR ptr = &lex;
term_order Lex(&Space, ptr);
ptr = &reversedeglex;
term_order ReverseDegLex(&Space, ptr);
ptr = &reverselex;
term_order ReverseLex(&Space, ptr);
//Default system
std::string gx = " 2 + n2 ";
std::string gy = " 1 + 78 * n1 + 2 * n2 + n4 * n1^2";
std::string gz = " 3 + n3 + n4";
std::string gu = " 4 * n1 - n2 * n3 + 6 * n1 * n4^3";
gPoly<double> px(&Space,gx,&Lex);
gPoly<double> py(&Space,gy,&Lex);
gPoly<double> pz(&Space,gz,&Lex);
gPoly<double> pu(&Space,gu,&Lex);
gPolyList<double> mylist(&Space, &ReverseDegLex);
mylist += px;
mylist += py;
mylist += pz;
mylist += pu;
PelView atlast(mylist);
gSpace NewSpace(3);
ptr = &lex;
term_order NewLex(&NewSpace, ptr);
ptr = &reversedeglex;
term_order NewReverseDegLex(&NewSpace, ptr);
ptr = &reverselex;
term_order NewReverseLex(&NewSpace, ptr);
//Default system
std::string newgx = " 2 + n2 ";
std::string newgy = " 1 + 78 * n1 + 2 * n2 + n3 * n1^2";
std::string newgz = " 3 + n3";
gPoly<double> newpx(&NewSpace,newgx,&NewLex);
gPoly<double> newpy(&NewSpace,newgy,&NewLex);
gPoly<double> newpz(&NewSpace,newgz,&NewLex);
gPolyList<double> mynewlist(&NewSpace, &NewReverseDegLex);
mynewlist += newpx;
mynewlist += newpy;
mynewlist += newpz;
PelView newatlast(mynewlist);
}
return 0;
}
syntax highlighted by Code2HTML, v. 0.9.1