Ticket #289: question.sing

// Clade distribution parameterizations
//
// balanced tree ((ab)c,(de))  
// pseudo-cat ((ab)(de),c)
// caterpillar ((((ab)c)d)e)
//
// modified from clade-dist-homog.sing on August 30, 2010
//
// In particular, the species trees were changed!

LIB "primdec.lib";
LIB "qhmoduli.lib";
LIB "solve.lib";
LIB "poly.lib";

system("--browser", "mac");
option("redSB");

// 1 = wordy; 0 = quiet
int wordy = 0;

// common ring for elimination ideals with parameters
ring r = 0,(x,y,z,H,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25), dp;

// Parameterization for the three non-caterpillar maps

// balanced ((ab)c,(de))
poly bal1 = 1-1/18*x*y^3*z-2/3*x-1/9*x*y^3;
poly bal2 = -1/18*x*y^3*z+1/3*x-1/9*x*y^3;
poly bal3 = 1/9*y*x*z+1/18*x*y^3*z;
poly bal4 = 1/9*y*x*z+1/18*x*y^3*z;
poly bal5 = -1/18*x*y^3*z+1/3*x-1/9*x*y^3;
poly bal6 = 1/9*y*x*z+1/18*x*y^3*z;
poly bal7 = 1/9*y*x*z+1/18*x*y^3*z;
poly bal8 = -1/9*y*x*z+1/18*x*y^3*z+2/9*y*z;
poly bal9 = -1/9*y*x*z+1/18*x*y^3*z+2/9*y*z;
poly bal10 = 1-2/3*z-1/18*y*x*z-1/9*y*z;
poly bal11 = 1-1/18*y*x*z+1/12*x*y^3*z-2/3*y-1/9*y*z-1/3*y*x+1/6*x*y^3;
poly bal12 = -1/9*y*x*z-1/36*x*y^3*z+2/9*y*z;
poly bal13 = -1/9*y*x*z-1/36*x*y^3*z+2/9*y*z;
poly bal14 = 1/9*y*x*z-1/36*x*y^3*z;
poly bal15 = 1/9*y*x*z-1/36*x*y^3*z;
poly bal16 = -1/12*y*x*z-1/18*x*y^3*z+1/6*y*x+1/18*x*y^3;
poly bal17 = 1/9*y*x*z-1/36*x*y^3*z;
poly bal18 = 1/9*y*x*z-1/36*x*y^3*z;
poly bal19 = -1/12*y*x*z-1/18*x*y^3*z+1/6*y*x+1/18*x*y^3;
poly bal20 = 1/12*y*x*z-1/18*x*y^3*z+1/3*y-1/6*y*z-1/6*y*x+1/18*x*y^3;
poly bal21 = 1/3*z-1/12*y*x*z+1/60*x*y^3*z-1/6*y*z;
poly bal22 = 1/3*z-1/12*y*x*z+1/60*x*y^3*z-1/6*y*z;
poly bal23 = 1/12*y*x*z+1/60*x*y^3*z+1/3*y-1/6*y*z-1/6*y*x;
poly bal24 = -1/12*y*x*z+1/60*x*y^3*z+1/6*y*x;
poly bal25 = -1/12*y*x*z+1/60*x*y^3*z+1/6*y*x;

ideal Bi = c1-bal1,c2-bal2,c3-bal3,c4-bal4,c5-bal5,c6-bal6,c7-bal7,c8-bal8,c9-bal9,c10-bal10,
           c11-bal11,c12-bal12,c13-bal13,c14-bal14,c15-bal15,c16-bal16,c17-bal17,c18-bal18,c19-bal19,c20-bal20,
           c21-bal21,c22-bal22,c23-bal23,c24-bal24,c25-bal25;

" ";
"Working to homogenize parameterizations ....";


ideal B = bal1,bal2,bal3,bal4,bal5,bal6,bal7,bal8,bal9,bal10,
          bal11,bal12,bal13,bal14,bal15,bal16,bal17,bal18,bal19,bal20,
          bal21,bal22,bal23,bal24,bal25;

// get homogeneous versions
list balDeg = deg(bal1), deg(bal2), deg(bal3), deg(bal8), deg(bal10), deg(bal11), deg(bal12), deg(bal14), 
             deg(bal16), deg(bal20), deg(bal21), deg(bal23), deg(bal24);
"   balDeg";
"    ",balDeg[1..size(balDeg)];  


" ";
"Working on 5-taxon balanced  ((ab)c,(de))";
ideal Hbal = B;
Hbal = homog(Hbal,H);
int sze = size(Hbal);

int m = Max(balDeg);
if (wordy==1){
  "The maximum degree before homogenization is", m;
}

//"Before making all the same degree...";
//"deg,  poly";
//for (int i=1; i<=sze; i=i+1){deg(Hbal[i]), "   ", Hbal[i];}
if (wordy==1)
  {
  "Those needing powers of H";
  for (int i=1; i<=sze; i=i+1)
Hcat    {
     if (m-deg(Hbal[i]))
       {
        i, "   ", Hbal[i];
       }
     }
  }

// homogenize polynomials
for (int i=1; i<=sze; i=i+1)
  {
   if (m-deg(Hbal[i]))
     {
      Hbal[i] = Hbal[i]*H^(m-deg(Hbal[i]));
     }
  }

if (wordy==1){
  "After making all the same degree.", deg(Hps[1..sze]); 
}

poly gg1 = c1-Hbal[1];
poly gg2 = c2-Hbal[2];
poly gg3 = c3-Hbal[3];
poly gg8 = c8-Hbal[8];
poly gg10 = c10-Hbal[10];
poly gg11 = c11-Hbal[11];
poly gg12 = c12-Hbal[12];
poly gg14 = c14-Hbal[14];
poly gg16 = c16-Hbal[16];
poly gg20 = c20-Hbal[20];
poly gg21 = c21-Hbal[21];
poly gg23 = c23-Hbal[23];
poly gg24 = c24-Hbal[24];

ideal JHbal = gg1, gg2, gg3, gg8, gg10, gg11, gg12, gg14, gg16, gg20, gg21, gg23, gg24;

// option(prot);
degBound=20;
"  For the balanced tree:  The degree bound is ", degBound;

// JHbal;
ideal stdJHbal = elim1(JHbal,x*y*z*H);
"  The size of stdJHbal is", size(stdJHbal);

" ";
"Summary:  ";
"  The size of stdJHbal is", size(stdJHbal);

"The degrees of the std's are ";
"  ", deg(stdJHbal[1..size(stdJHbal)]);


" ";
// option(prot);
degBound=10;
"**** Resetting degBound ****";
"  For the balanced tree:  The degree bound is ", degBound;

// JHbal;
ideal stdJHbal = elim1(JHbal,x*y*z*H);
"  The size of stdJHbal is", size(stdJHbal);

" ";
"Summary:  ";
"  The size of stdJHbal is", size(stdJHbal);

"The degrees of the std's are ";
"  ", deg(stdJHbal[1..size(stdJHbal)]);


" ";
// option(prot);
degBound=5;
"**** Resetting degBound ****";
"  For the balanced tree:  The degree bound is ", degBound;

// JHbal;
ideal stdJHbal = elim1(JHbal,x*y*z*H);
"  The size of stdJHbal is", size(stdJHbal);

" ";
"Summary:  ";
"  The size of stdJHbal is", size(stdJHbal);

"The degrees of the std's are ";
"  ", deg(stdJHbal[1..size(stdJHbal)]);


" ";
"THE END";

string cmd;
cmd = read("question.sing");