create_gml.c

#line 16 "create_gml.lw"
#include<LEDA/graph.h>
#include<LEDA/graph_alg.h>
#include<LEDA/sortseq.h>
#include<LEDA/stack.h>
#include<LEDA/set.h>
#include<LEDA/p_queue.h>
#include<LEDA/d2_dictionary.h>
#include<LEDA/map.h>
#include<LEDA/array2.h>
#include<LEDA/node_map2.h>
#include<LEDA/string.h>

#include<ctype.h>

using namespace LEDA;
using namespace std;

int add = 15;
int lp = 3;
int gp = 8;
int sp = 0;
int egp = 0;

class ABA_CPUTIMER {
  float mytime;
  float laststart;

public:
  ABA_CPUTIMER(void*) {
    mytime=0;
  }

  void start() {
    used_time(laststart);
  };

  void stop() {
    laststart=used_time(laststart);
    mytime+=laststart;
  };

  friend ostream& operator<<(ostream& o,const ABA_CPUTIMER& c) { 
    o<<c.mytime;
    return o; 
  };
};

double maxupperlower=10000;
bool fewvariables=true;
bool allnonpruned=false;
double mylowerbound=0;
bool writeoptvalue=true;
bool readoptvalue=false;
int edgebounce=4;
int gapbounce=3;
double inittolerance=0.98;
double pricingtolerance=0.98;
ABA_CPUTIMER heur(nil);
ABA_CPUTIMER extmixed(nil);
ABA_CPUTIMER extpair(nil);
ABA_CPUTIMER trans(nil);
ABA_CPUTIMER pgbuild(nil);
ABA_CPUTIMER prepro(nil);
ABA_CPUTIMER roottimer(nil);
double rootbound;
int allgvars;
int allavars;
int redavars;

LEDA::string filename;

class gap_edge {
 private:
  int s1, s2, l1, l2;
  double v;

 public:
  gap_edge(int s1t, int s2t, int l1t, int l2t) :
    s1(s1t), s2(s2t), l1(l1t), l2(l2t) {
  };

  gap_edge() {};

  int start_seq() const { return s1; };
  int end_seq() const { return s2; };
  int start_pos() const { return l1; };
  int end_pos() const { return l2; };

  double& value() { return v; };

  friend ostream& operator<<(ostream& o, const gap_edge&) {
    return o; }; 
  friend istream& operator>>(istream& i,const gap_edge&) {
    return i; };

  bool operator==(const gap_edge& ge) const {
    if(ge.start_seq()!=start_seq()) return false;
    if(ge.end_seq()!=end_seq()) return false;
    if(ge.start_pos()!=start_pos()) return false;
    if(ge.end_pos()!=end_pos()) return false;
    return true;
  };

  /*
  friend
  int compare(const gap_edge& g1, const gap_edge& g2) {
    if (compare(g1.s1, g2.s1)!=0) return compare(g1.s1, g2.s1);
    if (compare(g1.s2, g2.s2)!=0) return compare(g1.s2, g2.s2);
    if (compare(g1.l1, g2.l1)!=0) return compare(g1.l1, g2.l1);
    return compare(g1.l2, g2.l2);
  };
  */
};

namespace LEDA {
int compare(const gap_edge& g1, const gap_edge& g2) {
    if (compare(g1.start_seq(), g2.start_seq())!=0) return compare(g1.start_seq(), g2.start_seq());
    if (compare(g1.end_seq(), g2.end_seq())!=0) return compare(g1.end_seq(), g2.end_seq());
    if (compare(g1.start_pos(), g2.start_pos())!=0) return compare(g1.start_pos(), g2.start_pos());
    return compare(g1.end_pos(), g2.end_pos());
};
};

int Hash(const gap_edge& ge) {
  return 50000000*ge.start_seq()+1000000*ge.end_seq()+
    1000*ge.start_pos()+ge.end_pos();
};

int ID_Number(const gap_edge& ge) { 
  return Hash(ge);
}

LEDA::string Name[100];

class node_information {
 private:
  int p, s;
  char c;

 public:
  node_information(int pt, int st, char ct) :
    p(pt), s(st), c(ct) {
  };

  node_information() {};

  int pos() { return p; };
  int seq() { return s; };
  char letter() { return c; };

  friend ostream& operator<<(ostream& o, const node_information& i) {
    if(i.p==0) o<<Name[i.s]<<" ";
    o<<i.c;
    return o; }; 
  friend istream& operator>>(istream& i,const node_information&) {
    return i; };

};

edge_array<bool> LPEdge;
set<gap_edge> LPGap;
sortseq<gap_edge, double> GAP_COST;
sortseq<gap_edge, double> UPG;
GRAPH<node_information, two_tuple<int,double> > G;
edge_map<double> UPE(G);
node* Nodes[100];
node_map2<edge> AM(G);
int length[100];
LEDA::string strings[100];

edge AM_(node u, node v) {
  if(G[u].seq()>G[v].seq())
    return AM(u,v);
  else
    return AM(v,u);
}

double gap_cost(int k, int l, int i) {
  if(k>l) leda_swap(k,l);
  if(egp==0) {
    if((k==0) || (l==length[i]-1)) 
      return -sp*sqrt((float) l-k+1)-lp*(l-k+1) +add*(l-k+1);
  }
  return -gp-sp*sqrt((float) l-k+1)-lp*(l-k+1) +add*(l-k+1);
}

double** CM;
map<char,int> CMap;

double align_cost(int p1, int s1, int p2, int s2, const edge_array<two_tuple<int, double> >& EA) {
  if(s1<s2) return align_cost(p2,s2,p1,s1, EA);
  return EA[AM(Nodes[s1][p1],Nodes[s2][p2])].second();
};

double align_cost(char c1, char c2) {
  return CM[CMap[c1]][CMap[c2]]+2*add;
};

double gap_cost(gap_edge ge) {
  seq_item s=GAP_COST.lookup(ge);
  if(s==nil) {
    //cout<<"mil errror\n";
    //cout<<ge.end_pos()<<" "<<ge.start_pos()<<" "<<ge.end_seq()<<" "
    //<<length[ge.end_seq()]<<" "
    //<<gap_cost(ge.end_pos(),ge.start_pos(), ge.end_seq())<<" ge\n";
    return gap_cost(ge.end_pos(),ge.start_pos(), ge.end_seq());
  }
  return GAP_COST[s];
}

bool zero(const double& d) {
  return ((d>-0.001) && (d<0.001));
}

bool lzero(const double& d) {
  return d<-0.001;
}

bool lezero(const double& d) {
  return d<0.001;
}

#line 865 "create_gml.lw"
void make_pairgraph(GRAPH<node_information, double>& G, int i, int j,
                    GRAPH<edge,double>& H, int lgi, int lgj, edge_array<double>& ew,
                    node* lowest, node* highest, list<edge>& L, bool border_nodes,
                    node& u, int wi=-1, int wj=-1) 
{
  // Achtung, i Index wird Vertaucht i=lgi-i-1!

  pgbuild.start();

  H.clear();

  node_array<int> IJ(H, L.size()+2,0);
  node_array<int> IJ1(H, L.size()+2,0);
  node_array<double> V(H, L.size()+2, 0.0);

  edge e;
  node w; u=nil;
  bool s=false, t=false;
  forall(e, L) {
    w=H.new_node(e);
    if(G[source(e)].seq()==i)
      IJ[w]=lgi*G[target(e)].pos()+(lgi-G[source(e)].pos()-1);
    else
      IJ[w]=lgi*G[source(e)].pos()+(lgi-G[target(e)].pos()-1);
    if(IJ[w]==0) s=true;
    if(IJ[w]==lgj*lgi-1) t=true;
    V[w]=ew[e];
  }

  if(s==false) {
    w=H.new_node(nil);
    IJ[w]=0;
    V[w]=0;
  }

  if(t==false) {
    w=H.new_node(nil);
    IJ[w]=lgj*lgi-1;
    V[w]=0;
  }

  H.sort_nodes(IJ);

  for(int k=0; k<lgi; k++) {
    lowest[k]=nil;
    highest[k]=nil;
  }

  int i1,j1;
  forall_nodes(w, H) {   
    i1=IJ[w]%lgi; j1=-1;
    while(i1>=0) {
      if((lowest[lgi-i1-1]!=nil) && (IJ[lowest[lgi-i1-1]]/lgi>j1)) {
        H.new_edge(w, lowest[lgi-i1-1], -V[w]);
        j1=IJ[lowest[lgi-i1-1]]/lgi;
      }
      i1--;
    }
    lowest[lgi-IJ[w]%lgi-1]=w;
    if((highest[lgi-IJ[w]%lgi-1]==nil) && (IJ[w]>=0)) 
      highest[lgi-IJ[w]%lgi-1]=w;
    if((IJ[w]%lgi <= lgi-wi-1) && (IJ[w]/lgi <= wj)) u=w;
  }
  if(wi!=-1) 
    if(IJ[u]!=lgi-wi-1+wj*lgi) {
      cout<<"wrong u\n";
      cout<<IJ[u]%lgi<<" "<<IJ[u]/lgi<<" "<<wi<<" "<<wj<<endl;
    }
  
  pgbuild.stop();
}



#line 638 "create_gml.lw"
double pairwise_align(LEDA::string s1, LEDA::string s2, LEDA::string& s3, LEDA::string& s4, int n1, int n2,
                      const edge_array<two_tuple<int, double> >& EA, array2<double>& V) {
  int i,j,k, l1=s1.length(), l2=s2.length(); 
  array2<double> G(-1,l1,-1,l2), E(-1,l1,-1,l2), F(-1,l1,-1,l2);
  array2<int> TE(-1,l1,-1,l2), TF(-1,l1,-1,l2), TVE(-1,l1,-1,l2), TVF(-1,l1,-1,l2); // Trace 0=V(i-1,j-1);k=V(i,k);-k=V(k,j)

  V(-1,-1)=0;
  for(i=0; i<l1; i++) {
   F(i,-1)=gap_cost(gap_edge(n2,n1,0,i));
   V(i,-1)=F(i,-1);
   TF(i,-1)=i+1;
   TE(i,-1)=0;
   TVE(i,-1)=0;
   TVF(i,-1)=i+1;
   G(i,-1)=-MAXDOUBLE/2;
   E(i,-1)=-MAXDOUBLE/2;
  }
  for(j=0; j<l2; j++) {
   E(-1,j)=gap_cost(gap_edge(n1,n2,0,j));
   V(-1,j)=E(-1,j);
   TE(-1,j)=j+1;
   TF(-1,j)=0;
   TVE(-1,j)=j+1;
   TVF(-1,j)=0;
   G(-1,j)=-MAXDOUBLE/2;
   F(-1,j)=-MAXDOUBLE/2;
  }

  for(i=0; i<l1; i++) { 
    for(j=0; j<l2; j++) {
      G(i,j)=V(i-1,j-1)+align_cost(i,n1,j,n2, EA);
      E(i,j)=G(i,j-1)+gap_cost(gap_edge(n1,n2,j,j)); TE(i,j)=1;
      for(k=-1; k<j-1; k++) {
        if(E(i,j)<G(i,k)+gap_cost(gap_edge(n1,n2,k+1,j))) {
          E(i,j)=G(i,k)+gap_cost(gap_edge(n1,n2,k+1,j)); TE(i,j)=j-k; 
        }
      }
      F(i,j)=G(i-1,j)+gap_cost(gap_edge(n2,n1,i,i)); TF(i,j)=1;
      for(k=-1; k<i-1; k++) {
        if(F(i,j)<G(k,j)+gap_cost(gap_edge(n2,n1,k+1,i))) {
          F(i,j)=G(k,j)+gap_cost(gap_edge(n2,n1,k+1,i)); TF(i,j)=i-k;
        }
      }
      V(i,j)=G(i,j); TVE(i,j)=0; TVF(i,j)=0;
      for(k=-1; k<j; k++) {
        if(V(i,j)<F(i,k)+gap_cost(gap_edge(n1,n2,k+1,j))) {
          V(i,j)=F(i,k)+gap_cost(gap_edge(n1,n2,k+1,j));
          TVF(i,j)=TF(i,k);
          TVE(i,j)=j-k;
        }
      }
      if(E(i,j)>V(i,j)) { V(i,j)=E(i,j); TVE(i,j)=TE(i,j); TVF(i,j)=0; }
      if(F(i,j)>V(i,j)) { V(i,j)=F(i,j); TVF(i,j)=TF(i,j); TVE(i,j)=0; }

    }
  }

  cout<<V(l1-1,l2-1)<<" pairwise\n";

  s3=""; s4="";
  i=l1-1; j=l2-1;
  int de,df;
  
  double d=0;
  while((i>=0) || (j>=0)) {
    if(TVE(i,j)>0) {
      for(k=0; k<TVE(i,j); k++) {
        s3="-"+s3; s4=LEDA::string(s2[j-k])+s4; 
      }
      if(gap_cost(gap_edge(n1,n2,j-k+1,j))>0) cout<<j<<" "<<k<<" pos gap\n";
      d+=gap_cost(gap_edge(n1,n2,j-k+1,j));
      de=k;
    } else de=0;
    if(TVF(i,j)>0) {
      for(k=0; k<TVF(i,j); k++) {
        s3=LEDA::string(s1[i-k])+s3; s4="-"+s4;
      }
      if(gap_cost(gap_edge(n2,n1,i-k+1, i))>0) cout<<"pos gap2\n";
      d+=gap_cost(gap_edge(n2,n1,i-k+1, i));
      df=k;
    } else df=0;
    i-=df; j-=de;
    if((i>=0) && (j>=0)) {
      s3=LEDA::string(s1[i])+s3; s4=LEDA::string(s2[j])+s4;
      d+=align_cost(i,n1,j,n2,EA);
      i--; j--;
    }
  }
 
  cout<<d<<" test value\n";
  cout<<s3<<"\n"<<s4<<"\n";

  return V(l1-1,l2-1);
};

//TODO this is almost the same as the previous, but for an other objective function
// and does not compute the solution strings!
double pairwise_align(LEDA::string s1, LEDA::string s2, array2<double>& V, int i1, int i2) {
  int i,j,k, l1=s1.length(), l2=s2.length(); 
  array2<double> G(-1,l1,-1,l2), E(-1,l1,-1,l2), F(-1,l1,-1,l2);
  array2<int> TE(-1,l1,-1,l2), TF(-1,l1,-1,l2), TVE(-1,l1,-1,l2), TVF(-1,l1,-1,l2); // Trace 0=V(i-1,j-1);k=V(i,k);-k=V(k,j)

  V(-1,-1)=0;
  for(i=0; i<l1; i++) {
   F(i,-1)=gap_cost(0,i,i1);
   V(i,-1)=F(i,-1);
   TF(i,-1)=i+1;
   TE(i,-1)=0;
   TVE(i,-1)=0;
   TVF(i,-1)=i+1;
   G(i,-1)=-MAXDOUBLE/2;
   E(i,-1)=-MAXDOUBLE/2;
  }
  for(j=0; j<l2; j++) {
   E(-1,j)=gap_cost(0,j,i2);
   V(-1,j)=E(-1,j);
   TE(-1,j)=j+1;
   TF(-1,j)=0;
   TVE(-1,j)=j+1;
   TVF(-1,j)=0;
   G(-1,j)=-MAXDOUBLE/2;
   F(-1,j)=-MAXDOUBLE/2;
  }

  for(i=0; i<l1; i++) { 
    for(j=0; j<l2; j++) {
      G(i,j)=V(i-1,j-1)+align_cost(s1[i],s2[j]);
      E(i,j)=G(i,j-1)+gap_cost(j,j,i2); TE(i,j)=1;
      for(k=-1; k<j-1; k++) {
        if(E(i,j)<G(i,k)+gap_cost(j,k+1,i2)) {
          E(i,j)=G(i,k)+gap_cost(j,k+1,i2); TE(i,j)=j-k; 
        }
      }
      F(i,j)=G(i-1,j)+gap_cost(i,i,i1); TF(i,j)=1;
      for(k=-1; k<i-1; k++) {
        if(F(i,j)<G(k,j)+gap_cost(i,k+1,i1)) {
          F(i,j)=G(k,j)+gap_cost(i,k+1,i1); TF(i,j)=i-k;
        }
      }
      V(i,j)=G(i,j); TVE(i,j)=0; TVF(i,j)=0;
      for(k=-1; k<j; k++) {
        if(V(i,j)<F(i,k)+gap_cost(j,k+1,i2)) {
          V(i,j)=F(i,k)+gap_cost(j,k+1,i2);
          TVF(i,j)=TF(i,k);
          TVE(i,j)=j-k;
        }
      }
      if(E(i,j)>V(i,j)) { V(i,j)=E(i,j); TVE(i,j)=TE(i,j); TVF(i,j)=0; }
      if(F(i,j)>V(i,j)) { V(i,j)=F(i,j); TVF(i,j)=TF(i,j); TVE(i,j)=0; }

    }
  }

  cout<<V(l1-1,l2-1)<<" pairwise\n";

  return V(l1-1,l2-1);
};

// Pairwise sequence alignemnent with the cost function as first, but in
// reverse order! ohne rueckrechnung
double pairwise_align_rev(LEDA::string s1, LEDA::string s2, LEDA::string& s3, LEDA::string& s4, 
                          int n1, int n2,
                          const edge_array<two_tuple<int, double> >& EA, array2<double>& V) {
  int i,j,k, l1=s1.length(), l2=s2.length(); 
  array2<double> G(-1,l1,-1,l2), E(-1,l1,-1,l2), F(-1,l1,-1,l2);
  array2<int> TE(-1,l1,-1,l2), TF(-1,l1,-1,l2), TVE(-1,l1,-1,l2), 
    TVF(-1,l1,-1,l2); // Trace 0=V(i-1,j-1);k=V(i,k);-k=V(k,j)

  V(l1,l2)=0;
  for(i=0; i<l1; i++) {
   F(i,l2)=gap_cost(gap_edge(n2,n1,i,l1-1));
   V(i,l2)=F(i,l2);
   TF(i,l2)=l1-i;
   TE(i,l2)=0;
   TVE(i,l2)=0;
   TVF(i,l2)=l1-i;
   G(i,l2)=-MAXDOUBLE/2;
   E(i,l2)=-MAXDOUBLE/2;
  }
  for(j=0; j<l2; j++) {
   E(l1,j)=gap_cost(gap_edge(n1,n2,j,l2-1));
   V(l1,j)=E(l1,j);
   TE(l1,j)=l2-j;
   TF(l1,j)=0;
   TVE(l1,j)=l2-j;
   TVF(l1,j)=0;
   G(l1,j)=-MAXDOUBLE/2;
   F(l1,j)=-MAXDOUBLE/2;
  }

  for(i=l1-1; i>=0; i--) { 
    for(j=l2-1; j>=0; j--) {
      G(i,j)=V(i+1,j+1)+align_cost(i,n1,j,n2, EA);
      E(i,j)=G(i,j+1)+gap_cost(gap_edge(n1,n2,j,j)); TE(i,j)=1;
      for(k=j+1; k<=l2; k++) {
        if(E(i,j)<G(i,k)+gap_cost(gap_edge(n1,n2,j,k-1))) {
          E(i,j)=G(i,k)+gap_cost(gap_edge(n1,n2,j,k-1)); TE(i,j)=k-j; 
        }
      }
      F(i,j)=G(i+1,j)+gap_cost(gap_edge(n2,n1,i,i)); TF(i,j)=1;
      for(k=i+1; k<=l1; k++) {
        if(F(i,j)<G(k,j)+gap_cost(gap_edge(n2,n1,i, k-1))) {
          F(i,j)=G(k,j)+gap_cost(gap_edge(n2,n1,i, k-1)); TF(i,j)=k-i;
        }
      }
      V(i,j)=G(i,j); TVE(i,j)=0; TVF(i,j)=0;
      for(k=j+1; k<=l2; k++) {
        if(V(i,j)<F(i,k)+gap_cost(gap_edge(n1,n2,j, k-1))) {
          V(i,j)=F(i,k)+gap_cost(gap_edge(n1,n2,j, k-1));
          TVF(i,j)=TF(i,k);
          TVE(i,j)=k-j;
        }
      }
      if(E(i,j)>V(i,j)) { V(i,j)=E(i,j); TVE(i,j)=TE(i,j); TVF(i,j)=0; }
      if(F(i,j)>V(i,j)) { V(i,j)=F(i,j); TVF(i,j)=TF(i,j); TVE(i,j)=0; }

    }
  }

  cout<<V(0,0)<<" pairwise rev\n";

  return V(0,0);
};


#line 1027 "create_gml.lw"
void read_config_file() {
  ifstream is("config");
  LEDA::string s;
  
  while(!is.eof()) {
    is>>s;
    if(s=="FewVariables") is>>fewvariables;
    if(s=="AllNonPruned") is>>allnonpruned;
    if(s=="ReadLowerBound") is>>readoptvalue;
    if(s=="WriteOptValue") is>>writeoptvalue;
    if(s=="EdgeBounce") is>>edgebounce;
    if(s=="GapBounce") is>>gapbounce;
    if(s=="PricingTolerance") is>>pricingtolerance;
    if(s=="InitTolerance") is>>inittolerance;
    if(s=="MaxUpperLower") is>>maxupperlower;
    if(s=="GlobalPart") is>>gp;  
    if(s=="LinearPart") is>>lp;  
    if(s=="SquareRootPart") is>>sp;
    if(s=="EndGapPenality") is>>egp;  
    if(s=="Shift") is>>add; 

  }

  cout<<"\nConfiguration\n===============\n";
  cout<<"\nFewVariables: "<<fewvariables;
  cout<<"\nAllNonPruned: "<<allnonpruned;
  cout<<"\nReadLowerBound: "<<readoptvalue;
  cout<<"\nWriteOptValue: "<<writeoptvalue;
  cout<<"\nEdgeBounce: "<<edgebounce;
  cout<<"\nGapBounce: "<<gapbounce;
  cout<<"\nPricingTolerance: "<<pricingtolerance;
  cout<<"\nInitTolerance: "<<inittolerance<<"\n\n";
}

#line 244 "create_gml.lw"
int main(int argc, char* argv[]) {

  rootbound=-1;
  allgvars=0;
  allavars=0;
  redavars=0;
  prepro.start();

  if(argc<2) {
    cout<<"Usage: Infile(fasta) Outfile a b c (in a+b*l+c*sqrt(l))\n\n";
    return 1;
  }

  if(argc>2) {
    gp=atoi(argv[3]);
    lp=atoi(argv[4]);
    sp=atoi(argv[5]);
    //gp=atoi(argv[6]);
    add=gp+lp+sp;
    cout<<gp<<" "<<lp<<" "<<sp<<" "<<egp<<" "<<add<<" score data\n";
  }

  read_config_file();

  if (readoptvalue) {
    LEDA::string s=argv[1];
    s+=".opt";
    ifstream ifstr(s);
    ifstr>>mylowerbound;
    mylowerbound=mylowerbound*0.9999999;
    cout<<mylowerbound<<" given lower bound\n\n";
  }

  int sequences;

  
#line 357 "create_gml.lw"
ifstream is(argv[1]);

sequences=-1;

LEDA::string s;
while (!is.eof()) {
  s=read_line(is);
  if(s[0]=='>') {
   sequences++;
   Name[sequences]=s.del(0);
   int k=0;
   while((k<Name[sequences].length()) && (Name[sequences][k]!=','))
     k++;
   Name[sequences]=Name[sequences].del(k,Name[sequences].length());
   cout<<Name[sequences]<<" new sequence\n";
  } else if(sequences!=-1)
   strings[sequences]+=s;
}
sequences++;

cout<<"here\n";

for(int i=0; i<sequences; i++) {
  int k=0;
  while (k<strings[i].length()) {
    if((strings[i][k]>='A') && (strings[i][k]<='Z')) k++;
    else {
      if ((strings[i][k]>='a') && (strings[i][k]<='z')) {
        strings[i][k]=strings[i][k]+'A'-'a'; 
        k++;
      } else strings[i]=strings[i].del(k);
    }
  };
  cout<<strings[i]<<endl;
};

for(int j=0; j<sequences; j++) {
  Nodes[j]=new node[strings[j].length()];
  for(int i=0; i<strings[j].length(); i++) {
    Nodes[j][i]=G.new_node(node_information(i, j, strings[j][i]));
  }
  length[j]=strings[j].length();
}

cout<<sequences<<"\n reading done\n";

cout<<"read matrix\n";

ifstream matrix("matrix");
//string 
s="#";
while ((s[0]=='#') || (s.length()==0)) s.read_line(matrix);

int ns=0;
cout<<s<<" : list of characters\n";
for(int i=0; i<s.length(); i++) if (s[i]!=' ') {
  CMap[s[i]]=ns;
  ns++;
}

CM=new double*[ns+1];
for(int i=0; i<ns+1; i++) CM[i]=new double[ns+1];
for(int i=0; i<ns; i++) {
  matrix>>s; cout<<endl;
  for(int j=0; j<ns; j++) {
    matrix>>CM[i][j];
    //CM[i][j]=-CM[i][j];
    cout<<CM[i][j]<<" ";
  }
}

CMap['+']=ns;
for(int i=0; i<ns; i++) { CM[ns][i]=-10000; CM[i][ns]=-10000; }
CM[ns][ns]=1000;

cout<<"setup graph\n";

// ! Order of the edges is important for pairgraph_separator heuristic
node u, v;
forall_nodes(u, G) forall_nodes(v, G) if(G[u].seq()>G[v].seq()) {
  AM(u,v)=G.new_edge(u, v, two_tuple<int, double>(-1, align_cost(G[u].letter(),G[v].letter())));
}


cout<<"problem created\n\n";

//edge_array<double> UPE(G);
edge_array<double> fix(G);
set<int> GE[sequences][sequences];
LEDA::string ts1, ts2, s1t,s2t;
int l1, l2;
edge e;
double pa[sequences][sequences];
double hs=0;
for(int i=0; i<sequences; i++) for(int j=i+1; j<sequences; j++) {
  cout<<i<<" "<<j<<" sequences\n";
  l1=strings[i].length();
  l2=strings[j].length();
  array2<double> V1(-1,l1,-1,l2), V2(-1,l1,-1,l2);
  int p1=0, p2=0, g=0;
  s1t="";
  for(int k=l1-1; k>=0; k--) s1t+=strings[i][k];
  s2t="";
  for(int k=l2-1; k>=0; k--) s2t+=strings[j][k];
  pairwise_align(s1t,s2t,V2, i,j);
  pa[i][j]=pairwise_align(strings[i],strings[j],ts1,ts2,i,j, G.edge_data(), V1);
  for(int k=0; k<l1; k++) for(int l=0; l<l2; l++) {
    e=AM(Nodes[j][l],Nodes[i][k]);
    UPE[e]=V1(k-1,l-1)+V2(l1-k-2,l2-l-2)+G[e].second();
    fix[e]=V1(k,l)+V2(l1-k-2,l2-l-2);
    for(int N=leda_max(k-25,0); N<k; N++) 
      for(int M=k+1; M<leda_min(k+25,strings[i].length()); M++) {
        fix[e]=leda_max(fix[e], V1(N-1,l)+V2(l1-M-2,l2-l-2)+
                        gap_cost(M,N,strings[i].length()));
      }
    for(int N=leda_max(l-25,0); N<l; N++) 
      for(int M=l+1; M<leda_min(l+25,strings[j].length()); M++) {
        fix[e]=leda_max(fix[e], V1(k,N-1)+V2(l1-k-2,l2-M-2)+
                        gap_cost(M,N,strings[j].length()));
      }
if(fix[e]>pa[i][j]+0.001) cout<<"fix better than pairwise upper\n";
    // The best one can do if all letters before k are aligned before l and ...
  }

  for(int k=0; k<l2; k++) {
    for(int l=k; l<l2; l++) {
      double zg=0;
      if(k==0) zg=V2(l1-1,l2-l-2);
      if(l==l2-1) zg=V1(l1-1,k-1);
      if((k>0) && (l<l2-1)) {
        for(int z=-1; z<l1; z++) {
          zg=leda_max(zg, V1(z,k-1)+V2(l1-z-2,l2-l-2));
        }
      }
      gap_edge ge(i,j,k,l);
      UPG.insert(ge, zg+gap_cost(l,k,j));
    }
  }
  for(int k=0; k<l1; k++) {
    for(int l=k; l<l1; l++) {
      double zg=0;
      if(k==0) zg=V2(l1-l-2, l2-1);
      if(l==l1-1) zg=V1(k-1, l2-1);
      if((k>0) && (l<l1-1)) {
        for(int z=-1; z<l2; z++) {
          zg=leda_max(zg, V1(k-1,z)+V2(l1-l-2, l2-z-2));
        }
      }
      gap_edge ge(j,i,k,l);
      UPG.insert(ge, zg+gap_cost(l,k,i));
    }
  }

  hs+=pa[i][j];
  cout<<pa[i][j]<<" pw "<<hs<<" sum\n";
  pa[j][i]=pa[i][j];
  for(int k=0; k<ts1.length(); k++) {
    if((ts1[k]!='-') && (ts2[k]!='-')) {
      if(g>0) {
        gap_edge ge(i,j,p2-g,p2-1); // TODO exchanged i and j (4 times)
        if(g<=gapbounce) LPGap.insert(ge);
      }
      if(g<0) {
        gap_edge ge(j,i,p1+g,p1-1);
        if(-g<=gapbounce) LPGap.insert(ge);
      }
      GE[i][j].insert(1000*p1+p2);
      p1++;p2++; g=0;
    } else {
      if(ts1[k]=='-') {
        g++; p2++;
      } else {
        g--; p1++;
      }
    }
  }
  if(g>0) {
    gap_edge ge(i,j,p2-g,p2-1);
    if(g<=gapbounce) LPGap.insert(ge);
  }
  if(g<0) {
    gap_edge ge(j,i,p1+g,p1-1);
    if(-g<=gapbounce) LPGap.insert(ge);
  }
}

mylowerbound=leda_max(mylowerbound, hs-maxupperlower);

edge_array<double> UPE2(G);
node S,t;
forall_edges(e, G) {
  S=source(e); t=target(e);
  if((hs+UPE[e]-pa[G[S].seq()][G[t].seq()] >= mylowerbound) &&
     (G[S].pos()>0) && (G[t].pos()>0) &&
     (G[S].pos()<strings[G[S].seq()].length()-1) &&
     (G[t].pos()<strings[G[t].seq()].length()-1)) {
    double ub=0;
    for(int i=1; i<sequences; i++) 
      if ((i!=G[S].seq()) && (i!=G[t].seq())) {
        double m=fix[AM_(Nodes[i][strings[i].length()-1],
                         Nodes[G[S].seq()][G[S].pos()])]+
          fix[AM_(Nodes[i][strings[i].length()-1],
                  Nodes[G[t].seq()][G[t].pos()])];
          //gap_cost(G[S].pos()+1,strings[G[S].seq()].length()-1,G[S].seq())+
          //gap_cost(G[t].pos()+1,strings[G[t].seq()].length()-1,G[t].seq());

        for(int k=0; k<strings[i].length()-1; k++) {
          edge e1=AM_(Nodes[i][k],Nodes[G[S].seq()][G[S].pos()]),
            e2=AM_(Nodes[i][k],Nodes[G[t].seq()][G[t].pos()]);
          m=leda_max(m,fix[e1]+fix[e2]);
        }
        ub+=leda_min(0.0, m-pa[i][G[S].seq()]-pa[i][G[t].seq()]); // TODO required?
        //cout<< m-pa[i][G[S].seq()]-pa[i][G[t].seq()]<<" save over pairwise\n";
      }
    UPE2[e]=hs+UPE[e]-pa[G[S].seq()][G[t].seq()]+ub;
    //cout<<hs+UPE[e]-pa[G[S].seq()][G[t].seq()]<<" "<<UPE2[e]<<" old and new bound\n";
  } else UPE2[e]=hs+UPE[e]-pa[G[S].seq()][G[t].seq()];
}


forall_edges(e, G) {
  UPE[e]=hs+UPE[e]-pa[G[source(e)].seq()][G[target(e)].seq()];
  //UPE2[e]=UPE[e]; // TODO remove
}

cout<<hs<<" Obere Schranke\n";

cout<<"init LPedge\n";
LPEdge.init(G, !fewvariables);

l1=0; l2=0; int ltest=0;
forall_edges(e, G) {
  if(UPE2[e] < mylowerbound ) { // TODO UPE2 !
    l1++; 
    if(UPE[e] >= mylowerbound) {
      ltest++;
      //cout<<G[source(e)].seq()<<" "<<G[source(e)].pos()<<" "
      //  <<G[target(e)].seq()<<" "<<G[target(e)].pos()<<" edge\n";
    }
  } else { 
    l2++;
    if ((allnonpruned) || (UPE[e] >= hs*inittolerance)) LPEdge[e]=true;
  }
}
cout<<l1<<" pruned edges "<<l2<<" unpruned edges\n";
cout<<ltest<<" newly prunded\n";

forall_edges(e,G) if (!LPEdge[e]) G.hide_edge(e);

l1=0; l2=0;
if (allnonpruned) { // Take all needed gaps
  seq_item si;
  forall_items(si, UPG) {
    gap_edge ge=UPG.key(si);
    UPG.change_inf(si, hs+UPG.inf(si)-pa[leda_min(ge.start_seq(),ge.end_seq())]
                   [leda_max(ge.start_seq(),ge.end_seq())]);
    if((!LPGap.member(ge)) && (UPG.inf(si) >= mylowerbound )) {
      l2++;
      if(ge.end_pos()-ge.start_pos()<gapbounce) LPGap.insert(ge);
    } else l1++;
  }
} else {
  seq_item si;
  forall_items(si, UPG) { // Take all gaps close to pairwise optimal
    gap_edge ge=UPG.key(si);
    UPG.change_inf(si, hs+UPG.inf(si)-pa[leda_min(ge.start_seq(),ge.end_seq())]
                   [leda_max(ge.start_seq(),ge.end_seq())]);
    if((!LPGap.member(ge)) && (UPG.inf(si) >= hs*inittolerance)) {
      l2++;
      if(ge.end_pos()-ge.start_pos()<gapbounce) LPGap.insert(ge);
    } else l1++;
  }
}

G.restore_all_edges();
cout<<l1<<" pruned gaps "<<l2<<" unpruned gaps\n";

#line 279 "create_gml.lw"
                  ;

  
#line 289 "create_gml.lw"
//edge e;
if (!allnonpruned) {
  forall_edges(e, G) {
    if((GE[G[target(e)].seq()][G[source(e)].seq()].member(1000*G[target(e)].pos()
                                                          +G[source(e)].pos())) 
       || (abs(G[source(e)].pos()-G[target(e)].pos())<=edgebounce)) 
    LPEdge[e]=true;
  }

  int i, j, k, l;
  for(i=0; i<sequences; i++) {
    for(j=0; j<sequences; j++) if (i!=j) {
      for(k=0; k<length[j]; k++) {
        for(l=k; l<length[j]; l++) {
          gap_edge ge(i,j,k,l);
          if (l-k<gapbounce)
            LPGap.insert(ge);
        }
      }
    }
  }
}

forall_edges(e, G) if (!LPEdge[e]) G.del_edge(e);

/*
gap_edge ge;
forall(ge, LPGap) {
  G.new_edge(Nodes[ge.end_seq()][ge.start_pos()],
             Nodes[ge.end_seq()][ge.end_pos()], 
             two_tuple<int, double>(ge.start_seq(),
               gap_cost(ge.start_pos(), ge.end_pos(), ge.start_seq())));
}

for(int i=0; i<sequences; i++) {
  for(int j=1; j<length[i]; j++) {
    G.new_edge(Nodes[i][j-1], Nodes[i][j], 
               two_tuple<int, double>(-2, 0));
  }
};
*/

ofstream os1(argv[2]);
os1 << sequences <<" " << G.number_of_edges()+LPGap.size() << endl;

int SI[sequences];
SI[0]=0;
for(int i=0; i<sequences; i++) {
  os1 << length[i] <<" "<<Name[i] <<" "<<strings[i] <<endl;
  if(i<sequences-1) SI[i+1]=SI[i]+length[i];
}

forall_edges(e, G) {
  os1 << SI[G[source(e)].seq()]+G[source(e)].pos() <<" "
      << SI[G[target(e)].seq()]+G[target(e)].pos() <<" -1 "
      << G[e].second() << endl;
} 

gap_edge ge;
forall(ge, LPGap) {
  os1 << SI[ge.end_seq()]+ge.start_pos() <<" "
      << SI[ge.end_seq()]+ge.end_pos() <<" "
      << ge.start_seq() <<" " 
      << gap_cost(ge.start_pos(), ge.end_pos(), ge.start_seq()) <<endl;
}

#line 283 "create_gml.lw"
  return 0;
}

---