Changeset 8833bf in git for Singular


Ignore:
Timestamp:
May 8, 2015, 1:41:28 PM (9 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', '2a584933abf2a2d3082034c7586d38bb6de1a30a')
Children:
dfb374a725f78f1788adb5fd03badc36ec04c6a0
Parents:
52b3cf72ab7fa42d230232612144a2dd8a99aff0
Message:
maxlike.lib: comments, format
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/maxlike.lib

    r52b3cf7 r8833bf  
    196196static proc is_neg_def(matrix H)
    197197{//determines whether the given matrix is negative definite
    198    //returns 1 if it is, 0 if it isnt
     198  //returns 1 if it is, 0 if it isnt
    199199  matrix M=H-diag(var(1),ncols(H));
    200200  poly f=det(M);
    … …  
    300300  }
    301301
    302    //execute("ring outR =(complex,"+string(prec)+"),(x(1.."+string(n)+")),dp;");
     302  //execute("ring outR =(complex,"+string(prec)+"),(x(1.."+string(n)+")),dp;");
    303303  ring outR=(complex,prec),x(1..n),dp;
    304304  list MPOINTS = imap(R,L2);
    … …  
    379379"
    380380{//as opposed to (get)maxpoints, which first eliminates and then solves, this procedure
    381    //solves and then projects
    382    //furthermore, it also creates a list of the values the generators of I have at the
    383    //points in MPOINTS (that is, a list of the probability distributions)
     381  //solves and then projects
     382  //furthermore, it also creates a list of the values the generators of I have at the
     383  //points in MPOINTS (that is, a list of the probability distributions)
    384384  matrix H=logHessian(I,u);
    385385  def r=basering;
    … …  
    511511//   VALS, containing the probability distributions at those points.
    512512";
    513    if(size(#)==0) { print(MPOINTS); print(display); }
    514    return(outR);
     513  if(size(#)==0) { print(MPOINTS); print(display); }
     514  return(outR);
    515515}
    516516example
    … …  
    535535//////////////////////////////////////////////////////////////////////////////////////////
    536536//////////////////////////////////////////////////////////////////////////////////////////
     537//Here, we present an example of a data vector for which the likelihood function has more
     538//than one biologically meaningful local maximum.
     539//You can generate DNA sequence data, which has this data vector, using Seq-Gen with the
     540//following input:
     541//Tree: (Taxon1:0.6074796219,Taxon2:4.7911951859,Taxon3:0.5522879636);
     542//Seq-Gen options: -mHKY -l7647 -n1 -z28503
     543//You can find Seq-Gen at http://tree.bio.ed.ac.uk/software/seqgen/
     544//
     545//-   write(":w ThreeTaxonClaw.tree",
     546//-         "(Taxon1:0.6074796219,Taxon2:4.7911951859,Taxon3:0.5522879636);");
     547//-   int i=system("sh",
     548//-      "seq-gen -mHKY -l7647 -n1 -z28503 -q  < ThreeTaxonClaw.tree > ThreeTaxonC
     549//law.dat");
     550//-   intvec u=getintvec("ThreeTaxonClaw.dat");
     551//
    537552
    538553/*
    … …  
    550565    6*1/3*(1-mu1)*1/3*(1-mu2)*1/3*(1-mu3);
    551566   ideal I = f1,f2,f3,f4,f5;
    552    write(":w ThreeTaxonClaw.tree",
    553          "(Taxon1:0.6074796219,Taxon2:4.7911951859,Taxon3:0.5522879636);");
    554    int i=system("sh",
    555       "seq-gen -mHKY -l7647 -n1 -z28503 -q  < ThreeTaxonClaw.tree > ThreeTaxonClaw.dat");
    556    intvec u=getintvec("ThreeTaxonClaw.dat");
     567   intvec u = 770,2234,1156,2331,1156;
    557568   maxPoints(I,u,50);
    558569}
     570
    559571
    560572proc bad_seq_gen_example2()
    … …  
    616628*/
    617629
     630
    618631//////////////////////////////////////////////////////////////////////////////////////////
    619632//////////////////////////////////////////////////////////////////////////////////////////
    … …  
    623636//These are some of the procedures I used to generate and test examples for my Master
    624637//thesis. To use the ones incorporating Seq-Gen, you may have to adjust the shell command
    625 //with which Singular calls Seq-Gen. Furthermore, as the names of the procedures suggest,
    626 //we always use the Jukes-Cantor model.
     638//with which Singular calls Seq-Gen. You can find Seq-Gen at
     639//http://tree.bio.ed.ac.uk/software/seqgen/
     640//As the names of the procedures suggest, we always use the Jukes-Cantor model.
     641//They also tell you which of them use Seq-Gen.
    627642
    628643/*
    … …  
    687702    if(1)
    688703    {
    689        if(is_neg_def(hessi[k]) == 1)
    690        {
    691          hessi2=hessi2+list(hessi[k]);
    692          L2=L2+list(L[k]);
     704      if(is_neg_def(hessi[k]) == 1)
     705      {
     706        hessi2=hessi2+list(hessi[k]);
     707        L2=L2+list(L[k]);
    693708      }
    694709    }
    … …  
    725740    return(c);
    726741  }
     742
    727743  print(L2);
    728744}
     745
    729746
    730747proc getintvec(string linkstr)
    … …  
    769786  for(i=1; i<=size(taxons[1]); i++)
    770787  {
    771     if((taxons[1][i] == taxons[2][i]) and (taxons[2][i] == taxons[3][i]))
    772     {
    773       u[1]=u[1]+1;
    774       i++;
    775       continue;//continue does not execute the increment statement of the loop
    776     }
    777 
    778     if(taxons[1][i] == taxons[2][i])
    779     {
    780       u[3]=u[3]+1;
    781       i++;
    782       continue;
    783     }
    784 
    785     if(taxons[1][i] == taxons[3][i])
    786     {
    787       u[4]=u[4]+1;
    788       i++;
    789       continue;
    790     }
    791 
    792     if(taxons[2][i] == taxons[3][i])
    793     {
    794       u[5]=u[5]+1;
    795       i++;
    796       continue;
    797     }
    798 
    799     u[2]=u[2]+1;
    800   }
     788     if((taxons[1][i] == taxons[2][i]) and (taxons[2][i] == taxons[3][i]))
     789     {
     790        u[1]=u[1]+1;
     791        i++;
     792        continue;//continue does not execute the increment statement of the loop
     793     }
     794
     795     if(taxons[1][i] == taxons[2][i])
     796     {
     797        u[3]=u[3]+1;
     798        i++;
     799        continue;
     800     }
     801
     802     if(taxons[1][i] == taxons[3][i])
     803     {
     804        u[4]=u[4]+1;
     805        i++;
     806        continue;
     807     }
     808
     809     if(taxons[2][i] == taxons[3][i])
     810     {
     811        u[5]=u[5]+1;
     812        i++;
     813        continue;
     814     }
     815
     816     u[2]=u[2]+1;
     817  }
     818
    801819  return(u);
    802820}
    … …  
    851869    Iu=likeIdeal(I,u);
    852870
    853 
    854871    Iu=std(Iu);
    855872    //Iu=groebner(Iu);
    … …  
    857874    if(d != 0)
    858875    {
    859        nzd++;
    860        display="-*-*-*- not 0-dim. for u= "+string(u)+", i= "+string(i)+" -*-*-*-";
    861        print(display);
    862        write(lf,display);
    863        i++;
    864        continue;
     876      nzd++;
     877      display="-*-*-*- not 0-dim. for u= "+string(u)+", i= "+string(i)+" -*-*-*-";
     878      print(display);
     879      write(lf,display);
     880      i++;
     881      continue;
    865882    }
    866883
    … …  
    869886    if(eatoutput >= 2)
    870887    {
    871        num++;
     888      num++;
    872889      write(lf,writestr+"; number: "+string(eatoutput));
    873890      display="-*-*-*- Failed for u= "+string(u)+", i= "+string(i)+" -*-*-*-";
    … …  
    904921  ideal I=f1,f2,f3,f4,f5;
    905922
    906 
    907923  link lu=":a UsedIntvecsJC69seqgen.txt";
    908    link lf=":a FailedJC69seqgen.txt";
     924  link lf=":a FailedJC69seqgen.txt";
    909925
    910926  string readstr="ThreeTaxonClaw.dat";
    … …  
    922938  {
    923939    shcmd="seq-gen ";
    924      shcmd=shcmd+"-mHKY -l"+string(len)+" -n1 -z"+string(sd);
    925      shcmd=shcmd+" -q  < ThreeTaxonClaw.tree > ThreeTaxonClaw.dat";
    926      sd++;
     940    shcmd=shcmd+"-mHKY -l"+string(len)+" -n1 -z"+string(sd);
     941    shcmd=shcmd+" -q  < ThreeTaxonClaw.tree > ThreeTaxonClaw.dat";
     942    sd++;
    927943    eatoutput=system("sh",shcmd);
    928944    u=getintvec(readstr);
    … …  
    932948    Iu=likeIdeal(I,u);
    933949
    934 
    935950    Iu=std(Iu);
    936951    //Iu=groebner(Iu);
    … …  
    938953    if(d != 0)
    939954    {
    940        nzd++;
    941        display="-*-*-*- not 0-dim. for u= "+string(u)+", i= "+string(i)+" -*-*-*-";
    942        print(display);
    943        write(lf,display);
    944        i++;
    945        continue;
    946     }
    947 
     955      nzd++;
     956      display="-*-*-*- not 0-dim. for u= "+string(u)+", i= "+string(i)+" -*-*-*-";
     957      print(display);
     958      write(lf,display);
     959      i++;
     960      continue;
     961    }
    948962
    949963    eatoutput=getguaranteedmaxPoints(Iu,logHessian(I,u),1);
    … …  
    951965    if(eatoutput >= 2)
    952966    {
    953        num++;
     967      num++;
    954968      write(lf,writestr+"; number: "+string(eatoutput));
    955969      display="-*-*-*- no unique maximum for u= "+
    956            string(u)+", i= "+string(i)+" -*-*-*-";
     970      string(u)+", i= "+string(i)+" -*-*-*-";
    957971      print(display);
    958972      display="";
    … …  
    960974  }
    961975
    962 
    963976  display="-------------- i = "+string(i)+" --------------";
    964    display=display+newline+"not zero-dimensional in "+string(nzd)+" cases"+newline;
    965    display=display+"no unique maximum in "+string(num)+" cases"+newline;
    966    print(display);
    967    write(lf,display);
    968 
    969    close(lu);
     977  display=display+newline+"not zero-dimensional in "+string(nzd)+" cases"+newline;
     978  display=display+"no unique maximum in "+string(num)+" cases"+newline;
     979  print(display);
     980  write(lf,display);
     981
     982  close(lu);
    970983  close(lf);
    971984
    972    return(nzd,num);
     985  return(nzd,num);
    973986}
    974987
    975988proc randclawtree(int a)
    976989{
    977    ring r=(complex,10),x,dp;
    978    number n1,n2,n3;
    979    int n;
    980    int i;
    981    for(i=1; i<=a; i++)
    982    {
    983       n=random(1,1000000);
    984       n1=number(n)/1000000;
    985       n2=number(random(1,1000000-n))/1000000;
    986       n3=1-n1-n2;
    987       print("(Taxon1:"+string(n1)+",Taxon2:"+string(n2)+",Taxon3:"+string(n3)+");");
    988    }
     990  ring r=(complex,10),x,dp;
     991  number n1,n2,n3;
     992  int n;
     993  int i;
     994  for(i=1; i<=a; i++)
     995  {
     996    n=random(1,1000000);
     997    n1=number(n)/1000000;
     998    n2=number(random(1,1000000-n))/1000000;
     999    n3=1-n1-n2;
     1000    print("(Taxon1:"+string(n1)+",Taxon2:"+string(n2)+",Taxon3:"+string(n3)+");");
     1001  }
    9891002}
    9901003
    9911004proc checkrandomJC69writebeginning(int r, int a, int sta, int up)
    9921005{
    993    string writestr=newline+newline+newline+newline+newline+newline;
    994    writestr=writestr+"*****************************************************"+newline;
    995    writestr=writestr+"starting new loop with the following parameters"+newline;
    996    writestr=writestr+"number of runs: "+string(r)+newline;
    997    writestr=writestr+"number of intvecs per run: "+string(a)+newline;
    998    writestr=writestr+"starting random seed: "+string(sta)+newline;
    999    writestr=writestr+"upper bound for the entries of the intvecs: "+string(up)+newline;
    1000    writestr=writestr+"*****************************************************";
    1001 
    1002    link lu=":a UsedIntvecsJC69rand.txt";
    1003    link lf=":a FailedJC69rand.txt";
    1004    write(lu,writestr);
     1006  string writestr=newline+newline+newline+newline+newline+newline;
     1007  writestr=writestr+"*****************************************************"+newline;
     1008  writestr=writestr+"starting new loop with the following parameters"+newline;
     1009  writestr=writestr+"number of runs: "+string(r)+newline;
     1010  writestr=writestr+"number of intvecs per run: "+string(a)+newline;
     1011  writestr=writestr+"starting random seed: "+string(sta)+newline;
     1012  writestr=writestr+"upper bound for the entries of the intvecs: "+string(up)+newline;
     1013  writestr=writestr+"*****************************************************";
     1014
     1015  link lu=":a UsedIntvecsJC69rand.txt";
     1016  link lf=":a FailedJC69rand.txt";
     1017  write(lu,writestr);
    10051018  write(lf,writestr);
    10061019  close(lu);
    … …  
    10101023proc checkrandomJC69writeend(int r, int a, int sta, int up, int s, int t)
    10111024{
    1012    writestr=newline+newline+newline;
    1013    writestr=writestr+"*****************************************************"+newline;
    1014    writestr=writestr+"ending loop with the following parameters"+newline;
    1015    writestr=writestr+"number of runs: "+string(r)+newline;
    1016    writestr=writestr+"number of intvecs per run: "+string(a)+newline;
    1017    writestr=writestr+"starting random seed: "+string(sta)+newline;
    1018    writestr=writestr+"upper bound for the entries of the intvecs: "+string(up)+newline;
    1019    writestr=writestr+newline+"in the whole loop, there were a total of"+newline;
    1020    writestr=writestr+"   "+string(s)+" examples with non-zero-dim. likeideal"+newline;
    1021    writestr=writestr+"   "+string(t)+
     1025  writestr=newline+newline+newline;
     1026  writestr=writestr+"*****************************************************"+newline;
     1027  writestr=writestr+"ending loop with the following parameters"+newline;
     1028  writestr=writestr+"number of runs: "+string(r)+newline;
     1029  writestr=writestr+"number of intvecs per run: "+string(a)+newline;
     1030  writestr=writestr+"starting random seed: "+string(sta)+newline;
     1031  writestr=writestr+"upper bound for the entries of the intvecs: "+string(up)+newline;
     1032  writestr=writestr+newline+"in the whole loop, there were a total of"+newline;
     1033  writestr=writestr+"   "+string(s)+" examples with non-zero-dim. likeideal"+newline;
     1034  writestr=writestr+"   "+string(t)+
    10221035      " examples with more than one biol. meaningful local maximum";
    1023    writestr=writestr+newline+"*****************************************************";
    1024 
    1025 
    1026    write(lu,writestr);
     1036  writestr=writestr+newline+"*****************************************************";
     1037
     1038
     1039  write(lu,writestr);
    10271040  write(lf,writestr);
    10281041  close(lu);
    … …  
    10341047   //r the number of runs, a the number of intvecs per run, sta the starting random
    10351048   //seed, up the upper bound for the entries of the intvecs
    1036   checkrandomJC69writebeginning(r,a,sta,up);
     1049   checkrandomJC69writebeginning(r,a,sta,up);
    10371050
    10381051   int nzd, num, i, s, t;
    … …  
    10641077   link lf=":a FailedJC69seqgen.txt";
    10651078   write(lu,writestr);
     1079   write(lf,writestr);
     1080   close(lu);
     1081   close(lf);
     1082}
     1083
     1084
     1085proc checkseqgenJC69writeend(int r, int a, int sd, int sta,
     1086       int len, int p, int s, int t, intvec ls)
     1087{
     1088  writestr=newline+newline+newline;
     1089  writestr=writestr+"*****************************************************"+newline;
     1090  writestr=writestr+"ending loop with the following parameters"+newline;
     1091  writestr=writestr+"number of runs: "+string(r)+newline;
     1092  writestr=writestr+"number of intvecs per run: "+string(a)+newline;
     1093  writestr=writestr+"starting random seed for seqgen: "+string(sd)+newline;
     1094  writestr=writestr+"starting random seed for random: "+string(sta)+newline;
     1095  writestr=writestr+"length of the generated sequences: "+string(len)+newline;
     1096  writestr=writestr+newline+"in the whole loop, there were a total of"+newline;
     1097  writestr=writestr+"   "+string(s)+" examples with non-zero-dim. likeideal"+newline;
     1098  writestr=writestr+"   "+string(t)+
     1099        " examples with more than one biol. meaningful local maximum";
     1100  writestr=writestr+newline+"*****************************************************";
     1101  writestr=writestr+"used lengths:"+newline+string(ls);
     1102  writestr=writestr+newline+"*****************************************************";
     1103
     1104  write(lu,writestr);
    10661105  write(lf,writestr);
    10671106  close(lu);
    … …  
    10691108}
    10701109
    1071 
    1072 proc checkseqgenJC69writeend(int r, int a, int sd, int sta,
    1073        int len, int p, int s, int t, intvec ls)
    1074 {
    1075    writestr=newline+newline+newline;
    1076    writestr=writestr+"*****************************************************"+newline;
    1077    writestr=writestr+"ending loop with the following parameters"+newline;
    1078    writestr=writestr+"number of runs: "+string(r)+newline;
    1079    writestr=writestr+"number of intvecs per run: "+string(a)+newline;
    1080    writestr=writestr+"starting random seed for seqgen: "+string(sd)+newline;
    1081    writestr=writestr+"starting random seed for random: "+string(sta)+newline;
    1082    writestr=writestr+"length of the generated sequences: "+string(len)+newline;
    1083    writestr=writestr+newline+"in the whole loop, there were a total of"+newline;
    1084    writestr=writestr+"   "+string(s)+" examples with non-zero-dim. likeideal"+newline;
    1085    writestr=writestr+"   "+string(t)+
    1086         " examples with more than one biol. meaningful local maximum";
    1087    writestr=writestr+newline+"*****************************************************";
    1088    writestr=writestr+"used lengths:"+newline+string(ls);
    1089    writestr=writestr+newline+"*****************************************************";
    1090 
    1091 
    1092    write(lu,writestr);
    1093   write(lf,writestr);
    1094   close(lu);
    1095   close(lf);
    1096 }
    1097 
    1098 
    10991110proc checkseqgenJC69loop(int r, int a, int sd, int sta, int len, int p)
    11001111{
    1101    //r the number of runs, a the number of intvecs per run, sd the starting random
    1102    //seed, len the starting length, p the amount len will increase (on average)
    1103    //after each run (via + random(1,2*p-1))
    1104    //sta the random seed for random
    1105 
    1106    checkseqgenJC69writebeginning(r,a,sd,sta,len,p);
    1107 
    1108    system("random",sta);
    1109    intvec ls;
    1110    int nzd, num, i, s, t;
    1111    for(i=1; i<=r; i++)
    1112    {
    1113       (nzd,num)=checkseqgenJC69run(a,sd,len);
    1114       sd++;
    1115       ls[i]=len;
    1116       len=len+random(1,2*p-1);
    1117       s=s+nzd;
    1118       t=t+num;
    1119    }
    1120 
    1121    checkseqgenJC69writeend(r,a,sd,sta,len,p,s,t,ls);
     1112  //r the number of runs, a the number of intvecs per run, sd the starting random
     1113  //seed, len the starting length, p the amount len will increase (on average)
     1114  //after each run (via + random(1,2*p-1))
     1115  //sta the random seed for random
     1116
     1117  checkseqgenJC69writebeginning(r,a,sd,sta,len,p);
     1118
     1119  system("random",sta);
     1120  intvec ls;
     1121  int nzd, num, i, s, t;
     1122  for(i=1; i<=r; i++)
     1123  {
     1124    (nzd,num)=checkseqgenJC69run(a,sd,len);
     1125    sd++;
     1126    ls[i]=len;
     1127    len=len+random(1,2*p-1);
     1128    s=s+nzd;
     1129    t=t+num;
     1130  }
     1131
     1132  checkseqgenJC69writeend(r,a,sd,sta,len,p,s,t,ls);
    11221133}
    11231134*/
Note: See TracChangeset for help on using the changeset viewer.