#include <eutils/emain.h>
#include <eutils/ebasichashmap.h>
#include <eutils/eblockarray.h>
#include <eutils/estrhashof.h>
#include <eutils/einthashof.h>
#include <eutils/efile.h>
#include <eutils/etimer.h>
#include <eutils/eheap.h>
#include <eutils/ethread.h>

#include <map>
using namespace std;

#include "ekmerhashmap.h"
//#include "eseqali.h"

const unsigned long safe_shift[32]={0x0ul,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful,0xfffffffffffffffful};

/*
unsigned char seq_match_table[1u<<16u];

void initMatchTable()
{
  unsigned char tmp[4];
  tmp[0x00u]=1u;
  tmp[0x01u]=0u;
  tmp[0x02u]=0u;
  tmp[0x03u]=0u;
  
  for (uint32_t i=0u; i<(1u<<16u); ++i){
    seq_match_table[i]=0u;
    for (unsigned int k=0u; k<16u; k+=2u)
      seq_match_table[i]+=tmp[0x03u&(i>>k)];
  }
}
*/

unsigned char seq_comp_table[1u<<16u];

uint32_t nuc[]={'a','t','g','c','u'};
uint32_t compnuc[]={0x0u,0x1u,0x2u,0x3u,0x1u};

void initCompressionTable()
{
  for (uint32_t i=0; i<(1u<<16u); ++i)
    seq_comp_table[i]=0x0u;

  for (uint32_t i=0x0u; i<=0xffu; ++i){
    for (int j=0; j<5; ++j)
      seq_comp_table[(i<<8u) | nuc[j]]|=compnuc[j];
  }
  for (int i=0; i<5; ++i){
    for (uint32_t j=0x0u; j<=0xffu; ++j)
      seq_comp_table[(nuc[i]<<8u) | j]|=(compnuc[i]<<2u);
  }
}

class eseq
{
 public:
  estr useq;
  estr seq;
  int miss;
  int gaps;
  int seqlen;
  float quality;

  elongarray uniqind;
  eintarray uniqpos;
//  einthashof<int> seqs;
  eseq();
  eseq(estr& seq);
};

eseq::eseq(): seqlen(0),miss(0),gaps(0) {}

eseq::eseq(estr& ucseq): miss(0),gaps(0)
{
  int i;
  uint32_t tmp;
  useq.reserve(int((ucseq.len()+3)/4)*4);
  // align to 64bits + 64 empty bits
  int clen=int((ucseq.len()+7)/8+1)*8;
  seq.reserve(clen);
  ucseq.reserve(int((ucseq.len()+3)/4)*4);
  useq=ucseq;
  useq.lowercase();
  i=useq.find("nn");
  if (i!=-1) useq.del(i);
  for (i=useq.len()-1; (i>=0 && useq[i]=='n') || (i>0 && useq[i-1]=='n'); --i);
  if (i<0) { useq.clear(); seqlen=0; return; } // trim bad quality ends
  useq.del(i+1);

  quality=0.0;
  for (i=0; i<useq.len(); ++i) { if (useq[i]=='n') ++quality; }
  quality=1.0-quality/useq.len();

  useq.del(i);
  useq.replace("u","t");
  for (i=0; i<ucseq.len()-4; i+=4){
    tmp=*reinterpret_cast<uint32_t*>(&ucseq._str[i]);
    seq._str[i/4]=seq_comp_table[tmp&0xffffu]|(seq_comp_table[(tmp>>16)&0xffffu]<<4);
  }
  switch (ucseq.len()%4){
    case 3: 
      tmp=*(uint32_t*)(&ucseq._str[i]);
      seq._str[i/4]=seq_comp_table[tmp&0xffffu]|(seq_comp_table[(tmp>>16)&0x00ffu]<<4);
     break;
    case 2: 
      tmp=*(uint32_t*)(&ucseq._str[i]);
      seq._str[i/4]=seq_comp_table[tmp&0xffffu];
     break;
    case 1: 
      tmp=*(uint32_t*)(&ucseq._str[i]);
      seq._str[i/4]=seq_comp_table[tmp&0x00ffu];
     break;
  }
  for (i=i/4+1; i<clen; ++i)
    seq._str[i]=0x00;
  seq._strlen=clen;
  seqlen=ucseq.len();
}

uint32_t chr2kmer(const char *str)
{
  uint32_t *tmp=(uint32_t*)str;
  return(seq_comp_table[tmp[0]&0xffffu]|(seq_comp_table[(tmp[0]>>16u)&0xffffu]<<4u)|(seq_comp_table[tmp[1]&0xffffu]<<8u)|(seq_comp_table[(tmp[1]>>16u)&0xffffu]<<12u));
}

ostream& operator<<(ostream& stream,const eseq& seq)
{
  for (int i=0; i<seq.seq.len(); ++i){
    unsigned char t=seq.seq[i];
    for (int k=0; k<4; ++k,t>>=2u){
      switch (0x03u&t){
        case 0x00u: stream << "a"; break;
        case 0x01u: stream << "u"; break;
        case 0x02u: stream << "g"; break;
        case 0x03u: stream << "c"; break;
      }
    }
  }
  return(stream);
}

/*
int cmp_seq(const eseq& s1,const eseq& s2,int delta)
{
  uint32_t v,v2;
  int c=0u;
  int idelta;
  if (delta<0){
    delta=-delta;
    idelta=delta/4;
    v=uint32_t(s1.seq._str[0])|(uint32_t(s1.seq._str[1])<<8u);
    v2=uint32_t(s2.seq._str[idelta])|(uint32_t(s2.seq._str[idelta+1])<<8u);
    for (int j=0; j<s1.seq._strlen-4 && j+idelta<s2.seq._strlen-4; ++j){
      v=(v&0x0000ffffu)|(uint32_t(s1.seq._str[j+2])<<16u)|(uint32_t(s1.seq._str[j+3])<<24u);
      v2=(v2&0x0000ffffu)|(uint32_t(s2.seq._str[j+2+idelta])<<16u)|(uint32_t(s2.seq._str[j+3+idelta])<<24u);
      v2>>=2u*(delta%4);
      c+=seq_match_table[(v^v2)&0x0ffffu];
      v>>=16u,v2>>=16u-2u*(delta%4);
    }
  }else{
    idelta=delta/4;
    v=uint32_t(s1.seq._str[idelta])|(uint32_t(s1.seq._str[idelta+1])<<8u);
    v2=uint32_t(s2.seq._str[0])|(uint32_t(s2.seq._str[1])<<8u);
    for (int j=0; j+idelta<s1.seq._strlen-4 && j<s2.seq._strlen-4; ++j){
      v=(v&0x0000ffffu)|(uint32_t(s1.seq._str[j+2+idelta])<<16u)|(uint32_t(s1.seq._str[j+3+idelta])<<24u);
      v2=(v2&0x0000ffffu)|(uint32_t(s2.seq._str[j+2])<<16u)|(uint32_t(s2.seq._str[j+3])<<24u);
      v>>=2u*(delta%4);
      c+=seq_match_table[(v^v2)&0x0ffffu];
      v>>=16u-2u*(delta%4),v2>>=16u;
    }
  }
  return(c);
}
*/

//ebasicstrhashof<ebasicarray<ekmer> > kmer_hash;

inline void xy2estr(int x,int y,estr& str)
{
  str.clear();
  if (x<y){
    serialint(x,str);
    serialint(y,str);
  }else{
    serialint(y,str);
    serialint(x,str);
  }
}


#define IKMERSIZE 8ul
#define IKMERBITS (IKMERSIZE*2ul)
#define IKMERMAX (1ul<<IKMERBITS)
#define IKMERMASK (IKMERMAX-1ul)

extern unsigned char seqident_lt[IKMERMAX];

//unsigned char seqident_match[IKMERMAX];
//char seqident_delta[IKMERMAX];

/*
void initSeqIdent()
{
  for (uint32_t i=0; i<IKMERMAX; ++i){
    for (uint32_t j=0; j<IKMERMAX; ++j){
      seqident_match[i^j]=0u;
      uint32_t ti=i; uint32_t tj=j;
      unsigned char firstmiss=IKMERSIZE;
      for (int k=0; k<IKMERSIZE; ++k,tj>>=2u,ti>>=2u){
        if ((0x03u&ti)==(0x03u&tj))
          ++seqident_match[i^j];
        else if (firstmiss==IKMERSIZE) firstmiss=k;
      }
      seqident_match[i^j]|=(firstmiss<<4u);
    }
  }
}
*/
void seqident_aligned(const eseq& s1,const eseq& s2,int ip1,int ip2,int len,int& miss,int &gaps)
{
  int k;
  int p1=ip1;
  int p2=ip2;
  int ep1=ip1+len;
  int ep2=ip2+len;
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s1.seq._str);
  unsigned long *pstr2=reinterpret_cast<unsigned long*>(s2.seq._str);
  unsigned long v1,v2;

  unsigned char tmp,tmp2;
  unsigned char tmppos;
  char tmpgap;
  if (ep1>s1.seqlen) ep1=s1.seqlen;
  if (ep2>s2.seqlen) ep2=s2.seqlen;
  for (; p1<ep1-IKMERSIZE && p2<ep2-IKMERSIZE; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
    for (k=0; k<32u-IKMERSIZE; k+=IKMERSIZE,v1>>=IKMERSIZE*2u,v2>>=IKMERSIZE*2u){
      tmp=seqident_lt[(v1^v2)&IKMERMASK]; // first 4 bits are position of first mismatch, last 4 bits are matched nucleotides
      miss+=int(IKMERSIZE)-int(tmp&0xfu);
    }
  }
  if (p1<ep1 && p2<ep2){
    ldieif(MIN(ep1-p1,ep2-p2) > IKMERSIZE,"wrong");
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    tmp=seqident_lt[(v1^v2)&(IKMERMASK>>(IKMERBITS-MIN(ep1-p1,ep2-p2)*2u))]; // first 4 bits are position of first mismatch, last 4 bits are matched nucleotides
    miss+=int(IKMERSIZE)-int(tmp&0xfu);
  }
}

void seqident(const eseq& s1,const eseq& s2,int ip1,int ip2,int len,int& miss,int &gaps)
{
  int p1=ip1;
  int p2=ip2;
  int ep1=ip1+len;
  int ep2=ip2+len;
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s1.seq._str);
  unsigned long *pstr2=reinterpret_cast<unsigned long*>(s2.seq._str);
  unsigned long v1,v2;

  unsigned char tmp,tmp2;
  unsigned char tmppos;
  char tmpgap;
  if (ep1>s1.seqlen) ep1=s1.seqlen;
  if (ep2>s2.seqlen) ep2=s2.seqlen;
  for (; p1<ep1-IKMERSIZE && p2<ep2-IKMERSIZE; p1+=IKMERSIZE,p2+=IKMERSIZE){
    v1=pstr1[p1/32u]>>(2u*(p1%32u))|pstr1[p1/32u+1u]<<(64u-2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u))|pstr2[p2/32u+1u]<<(64u-2u*(p2%32u));
    tmp=seqident_lt[(v1^v2)&IKMERMASK]; // first 4 bits are position of first mismatch, last 4 bits are matched nucleotides
    if ((tmp&0xfu)!=IKMERSIZE) {
      tmppos=(tmp>>4u)&0xfu;
      ldieif(tmppos>=IKMERSIZE,"something wrong");
      v1>>=tmppos*2u;
      v2>>=tmppos*2u;
      tmpgap=1;
      tmp2=seqident_lt[((v1>>2u)^v2)&IKMERMASK];
      if ((tmp2&0xfu) != IKMERSIZE){
        tmpgap=2;
        tmp2=seqident_lt[((v1>>4u)^v2)&IKMERMASK];
      }
      if ((tmp2&0xfu) != IKMERSIZE){
        tmpgap=-1;
        tmp2=seqident_lt[(v1^(v2>>2u))&IKMERMASK];
      }
      if ((tmp2&0xfu) != IKMERSIZE){
        tmpgap=-2;
        tmp2=seqident_lt[(v1^(v2>>4u))&IKMERMASK];
      }
      if ((tmp2&0xfu) == IKMERSIZE){
        if (tmpgap<0) p2+=-tmpgap;
        else p1+=tmpgap;
        ++gaps;
      }else{
        miss+=int(IKMERSIZE)-int(tmp&0xfu);
      }
    }
  }
  if (p1<ep1 && p2<ep2){
    ldieif(MIN(ep1-p1,ep2-p2) > IKMERSIZE,"wrong");
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    tmp=seqident_lt[(v1^v2)&(IKMERMASK>>(IKMERBITS-MIN(ep1-p1,ep2-p2)*2u))]; // first 4 bits are position of first mismatch, last 4 bits are matched nucleotides
    miss+=int(IKMERSIZE)-int(tmp&0xfu);
  }
}

long mt_nj;
long mt_i;
long mt_j;
long mt_k;
long mt_m;
int nthreads=8;
int maxmiss=2;
emutex mt_mut;

estrarrayof<eseq> seqs;
earray<earray<elongarray> > dupslist;
elongarray seqsind;
elongarray seqscount;
int kmersize=50;

void taskSeqident()
{
//  cout << "Starting task" << endl;
  while(1) {
    mt_mut.lock();
    long i=mt_i;
    long m=mt_m;
    long nj=mt_nj;
    if (mt_k==dupslist[i][m].size()) { mt_mut.unlock(); return; }
  
    long bk=mt_k;
    long ek=bk + (dupslist[i][m].size()-bk)/nthreads/2;
    if (ek<bk+100) ek=bk+100;
    if (ek>dupslist[i][m].size()) ek=dupslist[i][m].size();
    mt_k=ek;
//    cout << "p: " << bk << " " << ek << endl;
    mt_mut.unlock();
    
    eseq &s(seqs.values(nj));
    for (long k=bk; k<ek; ++k){
      eseq &s2(seqs.values(dupslist[i][m][k]));
      if (seqscount[dupslist[i][m][k]]==0 || s.seqlen < s2.seqlen) continue;
      int miss=0;
      for (long l=0; miss<=maxmiss && l<s.uniqind.size() && l<s2.uniqind.size(); ++l){
        if (s.uniqind[l]!=s2.uniqind[l])
          ++miss;
      }
      if (miss<=maxmiss){
        int gaps=0;
        if (miss>0){
          miss=0;
          if (s.uniqind.size() && s2.uniqind.size() && s.uniqind[0]!=s2.uniqind[0])
            seqident_aligned(s,s2,0,0,kmersize,miss,gaps);
          for (int l=1; miss+gaps<=maxmiss && l<s.uniqind.size() && l<s2.uniqind.size(); ++l){
            if (s.uniqind[l]!=s2.uniqind[l]) {
              seqident_aligned(s,s2,s.uniqpos[l-1]+kmersize,s2.uniqpos[l-1]+kmersize,kmersize,miss,gaps);
            }
          }
        }
        if (miss+gaps<=maxmiss){
          s2.miss=miss; s2.gaps=gaps;
          seqsind[dupslist[i][m][k]]=nj;
          seqscount[nj]+=seqscount[dupslist[i][m][k]];
          seqscount[dupslist[i][m][k]]=0;
          break;
        }
      }
    }
  }
}




int emain()
{
  epregister(maxmiss);
  epregister(nthreads);
  int lcutoff=50;
  epregister(lcutoff);
  float qcutoff=0.98;
  epregister(qcutoff);

  bool aligned=false;
  epregister(aligned);

  int maxshift=5;
  epregister(maxshift);

  epregister(kmersize);

  estr badfile;
  epregister(badfile);
  estr of;
  epregister(of);

  int hashpos=0;
  int hashlen=8;
  epregister(hashpos);
  epregister(hashlen);
  eparseArgs();
  ldieif(getParser().args.size()<2,"syntax: denoiser <input>");
  if (aligned){
    maxshift=0;
  }

  initCompressionTable();
//  initMatchTable();
//  initSeqIdent();

  efile fbad;
  if (badfile.len())
    fbad.open(badfile,"w");
  
  efile f;
  estr line;
  estrarray args;

  estr str2id,str2seq;
  int maxseqlen=0;

  if (getParser().args[1]=="-"){
    f.open(stdin);
    cout << "# stdin reading sequences" << endl;
  }else
    f.open(getParser().args[1],"r");

  int totalseqs=0;
  estrarray parts;
  cout << "# reading sequences" << endl;
  seqs.reserve(60000000);
  f.readln(line);
  while (!f.eof()){
    ldieif(line.len()==0 || line[0]!='>',"line missing >: "+line);
    str2id=line;
    str2seq.clear();
    while (f.readln(line) && line.len() && line[0]!='>') { str2seq+=line; }
    str2seq.lowercase();
    eseq tmps(str2seq);
    parts=str2id.explode(" ");
    str2id=parts[0].substr(1);
    if (tmps.quality>=qcutoff && tmps.seqlen>=lcutoff)
      seqs.add(str2id,tmps);
    else if (badfile.len())
      fbad.write(str2id+" "+tmps.quality+" "+tmps.seqlen+" "+str2seq.len()+"\n");
    if (tmps.seqlen>maxseqlen) maxseqlen=tmps.seqlen;
    ++totalseqs;
  }
 
  cout << "# sequences: " << totalseqs << endl;
  cout << "# accepted: " << seqs.size() << endl;
  cout << "# removed: " << totalseqs-seqs.size() << endl;
  cout << "# maxlen: " << maxseqlen << endl;
  
  if(seqs.size()==0) { lerror("empty database"); return(0); }

  etimer t1;
  t1.reset();


  int maxdupsize=0;
  ebasicstrhashof<int> duphash;
  ebasicstrhashof<int>::iter it;
//  eintarray uniqind;
  ebasicstrhashof<int>::iter sit;
  duphash.reserve(seqs.size());

  // look for full sequence duplicates
  duphash.clear();

  seqsind.reserve(seqs.size());
//  for (int i=0; i<seqsind.size(); ++i) seqsind[i]=i;
  seqscount.init(seqs.size(),1);

  elongarray uniqind;
  uniqind.reserve(seqs.size());

  duphash.reserve(seqs.size());

  long i;
  for (i=0; i<seqs.size(); ++i){
    eseq &seq(seqs.values(i));
    if ((i%10000)==0)
      fprintf(stderr,"\r%li/%li",(long)i,(long)seqs.size());
    it=duphash.get(seq.useq);
    if (it==duphash.end()){
      duphash.add(seq.useq,i);
      uniqind.add(i);
      seqsind.add(i);
    } else {
      seqsind.add(it.value());
      seqscount[i]=0;
      ++seqscount[it.value()];
    }
  }
  fprintf(stderr,"\r%li/%li %li\n",(long)i,(long)seqs.size(),(long)duphash.size());
  cout << "# unique: " << duphash.size() << endl;

  elongarray seqdupcount;
  seqdupcount.reserve(uniqind.size());
  for (long i=0; i<uniqind.size(); ++i)
    seqdupcount.add(seqscount[i]); // add duplicate counts for every sequence

  elongarray si(lheapsort(seqdupcount));
  si.reverse();
  uniqind=uniqind[si];

  // look for shared kmers of *kmersize* length in the first kmer
  int k;
  for (k=0; k<kmersize; k+=kmersize){
    duphash.clear();
    dupslist.add(earray<elongarray>());
    earray<elongarray> &dlist(dupslist[dupslist.size()-1]);
    for (i=0; i<uniqind.size(); ++i){
      eseq &seq(seqs.values(uniqind[i]));
      if ((i%10000)==0)
        fprintf(stderr,"\r%i %li/%li",k,(long)i,(long)seqs.size());
      if (seq.useq.len()<k) continue;
      it=duphash.get(seq.useq.substr(k,kmersize));
      if (it==duphash.end()) {
        seq.uniqind.add(dlist.size());
        seq.uniqpos.add(k);
        duphash.add(seq.useq.substr(k,kmersize),dlist.size());
        dlist.add(elongarray(uniqind[i]));
      } else {
        elongarray &dups(dlist[it.value()]);
        seq.uniqind.add(it.value());
        seq.uniqpos.add(k);
        dups.add(uniqind[i]);
      }
    }
    fprintf(stderr,"\r%i %li/%li %li\n",k,(long)i,(long)uniqind.size(),(long)dlist.size());
  }

/*
  eintarray seqdupcount;
  seqdupcount.reserve(seqs.size());
  for (int i=0; i<uniqind.size(); ++i)
    seqdupcount.add(dupslist[0][seqs[uniqind[i]].uniqind[0]].size()); // add duplicate counts for every sequence

  eintarray si(iheapsort(seqdupcount));
*/

  // look for shared kmers of 50bp length every 50bp, start with more abundant (based on first 50bp) sequences first
  for (; k<maxseqlen; k+=kmersize){
    duphash.clear();
    dupslist.add(earray<elongarray>());
    earray<elongarray> &dlist(dupslist[dupslist.size()-1]);
    for (i=0; i<uniqind.size(); ++i){
      eseq &seq(seqs.values(uniqind[i]));
      if ((i%10000)==0)
        fprintf(stderr,"\r%i %li/%li",k,(long)i,(long)uniqind.size());
      int pos=seq.uniqpos[seq.uniqpos.size()-1]+kmersize;
      int z;
      if (seq.useq.len()<=pos+20) continue; // require at least 20 bp
      it=duphash.end();
      for (z=-maxshift; z<=maxshift && k+z>=0 && it==duphash.end(); ++z) it=duphash.get(seq.useq.substr(pos+z,kmersize));
      if (it==duphash.end()) {
        seq.uniqind.add(dlist.size());
        seq.uniqpos.add(pos);
        duphash.add(seq.useq.substr(pos,kmersize),dlist.size());
        dlist.add(elongarray(uniqind[i]));
      } else {
        elongarray &dups(dlist[it.value()]);
        seq.uniqind.add(it.value());
        seq.uniqpos.add(pos+z);
        dups.add(uniqind[i]);
      }
    }
    fprintf(stderr,"\r%i %li/%li %li\n",k,(long)i,(long)uniqind.size(),(long)dlist.size());
  }

  ethreads t;

  elongarray checkedseqs;
  for (long i=0; i<maxmiss+1 && i<dupslist.size(); ++i){
    checkedseqs.init(seqs.size(),-1l);
    long m;
    for (m=0; m<dupslist[i].size(); ++m){
      if (m%100==0)
        fprintf(stderr,"\r%li %li/%li",i,(long)m,(long)dupslist[i].size());
      for (long j=0; j<dupslist[i][m].size(); ++j){
        if (dupslist[i][m].size()>=10000 && j%10000==0)
          fprintf(stderr,"\r%li %li/%li %li/%li",i,(long)m,(long)dupslist[i].size(),j,(long)dupslist[i][m].size());
//        if (seqcount[dupslist[i][m][j]]==0) continue;
        long nj;
        for (nj=seqsind[dupslist[i][m][j]]; seqsind[nj]!=nj; nj=seqsind[nj]);
        seqsind[dupslist[i][m][j]]=nj;
        if (checkedseqs[nj]==m) continue; // this ensure we do not keep comparing sequences to already compared "references"
        checkedseqs[nj]=m;
        eseq &s(seqs.values(nj));
        if (dupslist[i][m].size()-j-1<100){
//        if (0){
          for (long k=j+1; k<dupslist[i][m].size(); ++k){
            eseq &s2(seqs.values(dupslist[i][m][k]));
            if (seqscount[dupslist[i][m][k]]==0 || s.seqlen < s2.seqlen) continue;
            int miss=0;
            for (long l=0; miss<=maxmiss && l<s.uniqind.size() && l<s2.uniqind.size(); ++l){
              if (s.uniqind[l]!=s2.uniqind[l])
                ++miss;
            }
  //          cout << "miss: " << miss << endl;
            if (miss<=maxmiss){
              int gaps=0;
              if (miss>0){
                miss=0;
                if (s.uniqind.size() && s2.uniqind.size() && s.uniqind[0]!=s2.uniqind[0])
                  seqident_aligned(s,s2,0,0,kmersize,miss,gaps);
                for (int l=1; miss+gaps<=maxmiss && l<s.uniqind.size() && l<s2.uniqind.size(); ++l){
                  if (s.uniqind[l]!=s2.uniqind[l]) {
                    seqident_aligned(s,s2,s.uniqpos[l-1]+kmersize,s2.uniqpos[l-1]+kmersize,kmersize,miss,gaps);
                  }
                }
              }
              if (miss+gaps<=maxmiss){
                s2.miss=miss; s2.gaps=gaps;
                seqsind[dupslist[i][m][k]]=nj;
                seqscount[nj]+=seqscount[dupslist[i][m][k]];
                seqscount[dupslist[i][m][k]]=0;
  /*
                if (miss+gaps>0){
                  cout << s.useq << " " << miss << " " << gaps << endl;
                  cout << s2.useq << " " << miss << " " << gaps << endl << endl;
                }
  */
                break;
              }
            }
          }
        }else{ // duplicate size larger than 100, go multithreaded
          mt_nj=nj;
          mt_i=i;
          mt_j=j;
          mt_k=j+1;
          mt_m=m;
         
//          fprintf(stderr,"\r%li %li/%li %li/%li p",i,(long)m,(long)dupslist[i].size(),j,(long)dupslist[i][m].size());

          t.run(taskSeqident,evararray(),nthreads);
//          fprintf(stderr,"\r%li %li/%li %li/%li wait",i,(long)m,(long)dupslist[i].size(),j,(long)dupslist[i][m].size());
          t.wait();
//          fprintf(stderr,"\r%li %li/%li %li/%li p done\n",i,(long)m,(long)dupslist[i].size(),j,(long)dupslist[i][m].size());
        }
      }
    }
    fprintf(stderr,"\r%li %li/%li\n",(long)i,(long)m,(long)dupslist[i].size());
  }

  elongarray uniqseqs;
  uniqseqs.init(seqs.size(),-1);
  long uniqseqscount=0;
  earray<elongarray> seqduplist;
  
  for (long i=0; i<seqs.size(); ++i){
    long n;
    for (n=seqsind[i]; seqsind[n]!=n; n=seqsind[n]);
    seqsind[i]=n;
    if (uniqseqs[n]==-1){
      uniqseqs[n]=uniqseqscount++;
      seqduplist.add(elongarray(n));
    }
    if (i==n) continue; // do not add two times the representative
    seqduplist[uniqseqs[n]].add(i);
  }
  cout << "# final: " << uniqseqscount << endl;

//  cout << "# finished computing in: " << t1.lap()/1000.0 << " secs" << endl;

  for (long i=0; i<seqduplist.size(); ++i){
//    eseq &sr(seqs.values(seqduplist[i][0]));
    cout << seqs.keys(seqduplist[i][0]) << "\t" << (seqscount[seqduplist[i][0]]-1);
    for (long j=1; j<seqduplist[i].size(); ++j){
//      eseq &s(seqs.values(seqduplist[i][j]));
      cout << "\t" << seqs.keys(seqduplist[i][j]);
    }
    cout << endl;
  }

  return(0);
}