#include <eutils/emain.h>
#include <eutils/ebasichashmap.h>
#include <eutils/estrhashof.h>
#include <eutils/einthashof.h>
#include <eutils/efile.h>
#include <eutils/etimer.h>
#include <eutils/eheap.h>
#include <eutils/ernd.h>
#include <eutils/ethread.h>
#include <eutils/esystem.h>
#include <eutils/eoption.h>
#include <deque>
#include <algorithm>

#define LDEBUG(key,cmd) LDEBUG2(key,cmd)
#define LDEBUG2(key,cmd) LDEBUG_##key(cmd)
#define LDEBUG_(cmd) {}
#define LDEBUG_TRUE(cmd) cmd

#define D_SEQALIGNMENT
#define D_SEQIDENT
#define D_SEQSEARCH

//#define D_SEQIDENT TRUE
//#define D_SEQALIGNMENT TRUE
//#define D_SEQSEARCH TRUE


#include <math.h>

#include <map>
using namespace std;

#include "mapclust-config.h"
#include "ekmerhashmap.h"
#include "eseqali.h"

#include "kmerseqtables-data.h"

#define KMERSIZE 10ul
#define KMERBITS (KMERSIZE*2ul)
#define MAXSIZE (1ul<<KMERBITS)
#define KMERMAX (1ul<<KMERBITS)
#define KMERMASK (MAXSIZE-1ul)

#define KMERSIZE2 8ul
#define KMERBITS2 (KMERSIZE2*2ul)
#define MAXSIZE2 (1ul<<KMERBITS2)
#define KMERMAX2 (1ul<<KMERBITS2)
#define KMERMASK2 (MAXSIZE2-1ul)



const uint32_t nuc[]={'a','t','g','c','u'};
const uint32_t compnuc[]={0x0u,0x1u,0x2u,0x3u,0x1u};
const uint32_t consnuc[]={0x0001u,0x0010u,0x0100u,0x1000u};
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 int akmers[MAXSIZE]; // 6bases

int topotus=60;
int tophits=20;
int otulim=50;
float ccutoff=0.96;

/*
const float matchcost=2.0;
const float misscost=-1.0;
const float gapcost=-1.0;
const float gapopen=-5.0;
*/

//const ealignscore as(1.0,2.0,10.0,1.0);
//const ealignscore as(2.0,1.0,5.0,1.0,30.0);
const ealignscore as(1.0,1.0,5.0,1.0,20.0);
//const ealignscore as(2.0,1.0,5.0,1.0,30.0);

const int NCOUNT_MAXLEN=1600; // IMPORTANT: must be disivible by 8 to be 64bit aligned

float ti=0.0,ts=0.0,ti2=0.0,ts2=0.0,ta=0.0;
float tdp=0.0,tdp1=0.0,tdp2=0.0,tdpfl=0.0,tdpmd=0.0;
etimer t1,t2;

const uint32_t BMASK31=(1u<<31u)-1u;
const uint32_t BMASK16=(1u<<16u)-1u;
const uint32_t BMASK11=(1u<<11u)-1u;
const uint32_t BMASK10=(1u<<10u)-1u;


class etax
{
 public:
  estr name;
  earray<estr> levels;
  earray<estrarray> names;
  earray<estrhashof<int> > ind;
  ebasicarray<eseqtax*> seqs;
  efloatarray cutoff;
  efloatarray cutoffcoef;
};

class esearchws
{
 public:
//  eintarray best;
  ealignws alignws;
  eintarray otukmerpos;
//  uint64_t *kmerbitmask;
//  uint64_t *bitmask;
  ebasicarray<uint32_t> idcount;
  ebasicarray<uint32_t> idlist;
//  ebasicarray<uint32_t> idcount2;
  int maskid;
  eintarray kmermask;
  int offset;
  eintarray kmerpos;
  eintarray kmerposrev;
  int offset2;
  eintarray kmerpos2;
  ebasicarray<eintarray> taxcounts;
  eintarray tmpkmers;

  void init(int mxsize);
};


class ediag
{
 public:
  long i;
  long j;
  long i2;
  long j2;
  float V;
  ediag *bestseg;
  ediag(): i(-1),j(-1),i2(-1),j2(-1),V(-1),bestseg(0x00) {}
//  ediag(int _i,int _j,int _i2,int _j2): i(_i),j(_j),i2(_i2),j2(_j2),V(-1),bestseg(0x00) {}
  ediag(long _i,long _delta,long _len): i(_i),j(_i+_len),i2(_i+_delta),j2(_i+_delta+_len),V(-1),bestseg(0x00) {}
};

/*
class ekmer
{
 public:
  int i;
//  eintarray pos;
  int pos;
  ekmer(): i(-1) {}
  ekmer(int _i,int _pos): i(_i),pos(_pos) {}
//  ekmer(int _i,eintarray& _pos): i(_i),pos(_pos) {}
};
*/

class emaxcount {
 public:
  int i;
  int maxcount;
  eintarray kwordaln;
  emaxcount(): i(-1), maxcount(0) {}
};

char oct2hex(unsigned char o)
{
  switch (o){
    case 0x00: return('0');
    case 0x01: return('1');
    case 0x02: return('2');
    case 0x03: return('3');
    case 0x04: return('4');
    case 0x05: return('5');
    case 0x06: return('6');
    case 0x07: return('7');
    case 0x08: return('8');
    case 0x09: return('9');
    case 0x0a: return('a');
    case 0x0b: return('b');
    case 0x0c: return('c');
    case 0x0d: return('d');
    case 0x0e: return('e');
    case 0x0f: return('f');
  }
  ldie("wrong octet value");
  return('-');
}


char* uint2hex(char *tmpstr,unsigned int c)
{
  tmpstr[0u]='0';
  tmpstr[1u]='x';
  for (unsigned int k=0; k<sizeof(c)*8u/4u; ++k)
    tmpstr[2u+k]=oct2hex((c>>(sizeof(c)*8u-4u-k*4u))&0xf);
  tmpstr[2u+sizeof(c)*8u/4u]='u';
  tmpstr[2u+sizeof(c)*8u/4u+1u]=0x00;
  return(tmpstr);
}


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));
}

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



estr long2str(unsigned long kmer)
{
  estr tmpstr;
  for (int i=0; i<64u; i+=2u,kmer>>=2u){
    switch (0x03u&kmer){
      case 0x00u: tmpstr+='A'; break;
      case 0x01u: tmpstr+='T'; break;
      case 0x02u: tmpstr+='G'; break;
      case 0x03u: tmpstr+='C'; break;
    }
  }
  return(tmpstr);
}
estr kmer2str(uint32_t kmer)
{
  estr tmpstr;
  for (int i=0; i<IKMERSIZE; ++i,kmer>>=2u){
    switch (0x03u&kmer){
      case 0x00u: tmpstr+='A'; break;
      case 0x01u: tmpstr+='T'; break;
      case 0x02u: tmpstr+='G'; break;
      case 0x03u: tmpstr+='C'; break;
    }
  }
  return(tmpstr);
}

class eseqdb
{
 public:
  estrarrayof<eseq> seqs;
  estrhashof<int> seqind;
  eintarray seqotu2;
  earray<eintarray> otus2;
  eintarray seqotu;
  earray<eintarray> otus;
  ebasicarray<deque<uint32_t> > otukmers;
  ebasicarray<char> otukmermask;
  ebasicarray<etax> taxa;
};

unsigned long seqkmer(const eseq& s,long p1)
{
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s.seq._str);
  return(((pstr1[p1/32u]>>(2u*(p1%32u)))|((pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u]))&KMERMASK);
}

unsigned long seqrevkmer(const eseq& s,long p1)
{
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s.seq._str);
  p1=s.seqlen-p1;
  return(kmer_rev_lt[((pstr1[p1/32u]>>(2u*(p1%32u)))|((pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u]))&KMERMASK]);
}

void print_tncount(unsigned char* tncounts,int p,int e)
{
  uint64_t vid;

  for (int i=p; i<e; ++i){
    vid=tncounts[i];
    switch (vid&0x3u){ 
      case 0x0u: cout << "i: " << i << " -" << endl; break;
      case 0x1u: cout << "i: " << i << " m" << endl; break;
      case 0x2u: cout << "i: " << i << " x" << endl; break;
    }
  }
}

unsigned char d1_lt[1u<<16u];
unsigned char d2_lt[1u<<16u];
unsigned char id_lt[1u<<16u];


void initdlt()
{
  for (uint32_t i=0u; i<(1u<<16u); ++i){
    id_lt[i]=0u;
    d1_lt[i]=8u;
    d2_lt[i]=8u;
  }
  id_lt[0u]=1u;
  d2_lt[0u]=0u;
  for (uint32_t i=0u; i<8u; ++i)
    d1_lt[i]=i;
}

void seqident_seg_left(const eseq& s1,int p1,const eseq& s2,int p2,ealigndata& adata,const ealignscore& as){
  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[5];
  int k;
  int itmp;

  if (p1==0 || p2==0) return;

//  cout << "seg_left id: " << p1 << " :: " << p2 << " len1: " << p1 << " len2: " << p2 << endl;
  for (; p1>=32 && p2>=32; p1+=k-32u+IKMERSIZE,p2+=k-32u+IKMERSIZE){
    v1=pstr1[p1/32u-1u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u-1u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//    cout << p1 << " " << long2str(v1) << endl;
//    cout << p2 << " " << long2str(v2) << endl;
    for (k=32u-IKMERSIZE; k>=0; k-=IKMERSIZE,v1<<=IKMERSIZE*2u,v2<<=IKMERSIZE*2u){
//      cout << p1 << "+" << k << " " << kmer2str((v1>>((32u-IKMERSIZE)*2u))&IKMERMASK) << endl; 
//      cout << p2 << "+" << k << " " << kmer2str((v2>>((32u-IKMERSIZE)*2u))&IKMERMASK) << endl; 
      tmp[0]=seq_ident_table[((v1^v2)>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      tmp[1]=seq_ident_table[((v1^(v2>>2ul))>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      tmp[2]=seq_ident_table[((v1^(v2>>4ul))>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      tmp[3]=seq_ident_table[((v2^(v1>>2ul))>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      tmp[4]=seq_ident_table[((v2^(v1>>4ul))>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      
      itmp=0;
      for (int i=1; i<5; ++i){
        if (tmp[i]>tmp[itmp]) itmp=i;
      }
      if (itmp!=0 && tmp[itmp]>tmp[0]+IKMERSIZE/4){
//        matches+=tmp[0]; mismatches+=IKMERSIZE-tmp[0];
        switch(itmp){
          case 1: p1-=1; adata.gaps+=1; break;
          case 2: p1-=2; adata.gaps+=2; break;
          case 3: p2-=1; adata.gaps+=1; break;
          case 4: p2-=2; adata.gaps+=2; break;
        }
        break;
      }else{
        adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
        adata._score+=tmp[0]*as.match-(IKMERSIZE-tmp[0])*as.mismatch;
      }
    }
  }
//  cout << "seg_left rest id: " << p1 << " :: " << p2 << " len1: " << p1 << " len2: " << p2 << endl;
  for (;p1>=int(IKMERSIZE) && p2>=int(IKMERSIZE); p1+=k+int(IKMERSIZE)-32,p2+=k+int(IKMERSIZE)-32){ // needed to handle gap shifts gracefully
//    cout << p1 << endl;
//    cout << p2 << endl;
    v1=(p1>=32?pstr1[p1/32u-1u]>>(2u*(p1%32u)):0ul);
    v2=(p2>=32?pstr2[p2/32u-1u]>>(2u*(p2%32u)):0ul);
    v1|=(pstr1[p1/32u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//    cout << p1 << " " << long2str(v1) << endl;
//    cout << p2 << " " << long2str(v2) << endl;
    for (k=32u-IKMERSIZE; k>=0 && p1+k>=int(IKMERSIZE) && p2+k>=int(IKMERSIZE); k-=IKMERSIZE,v1<<=IKMERSIZE*2u,v2<<=IKMERSIZE*2u){
//      cout << p1 << "+" << k << " " << kmer2str((v1>>((32u-IKMERSIZE)*2u))&IKMERMASK) << endl; 
//      cout << p2 << "+" << k << " " << kmer2str((v2>>((32u-IKMERSIZE)*2u))&IKMERMASK) << endl; 
      tmp[0]=seq_ident_table[((v1^v2)>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      tmp[1]=seq_ident_table[((v1^(v2>>2ul))>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      tmp[2]=seq_ident_table[((v1^(v2>>4ul))>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      tmp[3]=seq_ident_table[((v2^(v1>>2ul))>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      tmp[4]=seq_ident_table[((v2^(v1>>4ul))>>((32u-IKMERSIZE)*2u))&IKMERMASK];
      
      itmp=0;
      for (int i=1; i<5; ++i){
        if (tmp[i]>tmp[itmp]) itmp=i;
      }
      if (itmp!=0 && tmp[itmp]>tmp[0]+IKMERSIZE/4){
        switch(itmp){
          case 1: p1-=1; adata.gaps+=1; break;
          case 2: p1-=2; adata.gaps+=2; break;
          case 3: p2-=1; adata.gaps+=1; break;
          case 4: p2-=2; adata.gaps+=2; break;
        }
        break;
      }else{
        adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
        adata._score+=tmp[0]*as.match-(IKMERSIZE-tmp[0])*as.mismatch;
      }
    }
  }

  int l=(p1<p2?p1:p2);
//  cout << "missing: " << p1 << " l: " << l << endl;
  if (l<=0) return;
  p1-=l;
  p2-=l;
  v1=pstr1[p1/32u]>>(2u*(p1%32u));
  v2=pstr2[p2/32u]>>(2u*(p2%32u));
  v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
  v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;

  tmp[0]=seq_ident_table[(v1^v2)&((0x1ul<<(l*2ul))-1ul)];
//  cout << long2str(v1&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
//  cout << long2str(v2&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
  adata.matches+=tmp[0]-(IKMERSIZE-l);
  adata.mismatches+=IKMERSIZE-tmp[0];
  adata._score+=(tmp[0]-(IKMERSIZE-l))*as.match-(IKMERSIZE-tmp[0])*as.mismatch;

/*
  tmp[0]=seq_ident_table[((v1^v2)>>((32u-IKMERSIZE)*2u))&~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul)&IKMERMASK]; // put mask at end of 64bit int, i.e.: FF000000

  cout << long2str((v1>>((32u-IKMERSIZE)*2u))&(~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul))&IKMERMASK) << " " << long2str((~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul))&IKMERMASK) << endl;
  cout << long2str((v2>>((32u-IKMERSIZE)*2u))&~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul)&IKMERMASK) << " " << long2str(~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul)&IKMERMASK) << endl;
  matches+=tmp[0]-(IKMERSIZE-l);
  mismatches+=IKMERSIZE-tmp[0];
*/
}


bool galign=true;


inline void sumcounts(unsigned char* tncount,int p1,uint64_t tmpnuc)
{
  *reinterpret_cast<uint64_t*>(&tncount[p1])|=tmpnuc;
/*
  for (int i=0; i<8; ++i,tmpnuc>>=8u)
    tncount[p1+i]|=(tmpnuc&0xff);
*/
//  cout << uint2hex(tmpstr,*p) << " " << uint2hex(tmpstr2,tmpnuc) << " " << uint2hex(tmpstr3,tncount[p1]) << endl;
}

void seqncount(const eseq& s1,long p1,long e1,const eseq& s2,long p2,long e2,ealigndata& adata,const ealignscore& as)
{
  int k;
  adata.matches+=e1-p1;
  adata._score+=(e1-p1)*as.match;
  adata.profile.add(AT_MATCH,e1-p1);
  if (adata.aln==0x00) return;

/*
  uint64_t *pstr1=reinterpret_cast<uint64_t*>(s1.seq._str);
  uint64_t *pstr2=reinterpret_cast<uint64_t*>(s2.seq._str);
  uint64_t v1,v2;
*/

//  cout << "id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p2+IKMERSIZE<=e2; p2+=IKMERSIZE,p1+=IKMERSIZE){
//    cout << "sumcounts: " << p1 << endl;
//    v1=(pstr1[p1/32u]>>(2u*(p1%32u)))|((pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u]);
//    v2=(pstr2[p2/32u]>>(2u*(p2%32u)))|((pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u]);
    sumcounts(adata.aln,p1,seq_alignment_lt[0u]);
//    tmpas1+=kmer2str(v1);
//    tmpas2+=kmer2str(v2);
  }

  int l=e2-p2;
  if (l<=0) return;
//  if (l<=0) goto fend;
//  cout << "missing: p1: " << p1 << " p2: " << p2 << " l: " << l << endl;
  
  sumcounts(adata.aln,p1,seq_alignment_lt[0u]&((0x1ul<<(l*8u))-1ul));
//  v1=(pstr1[p1/32u]>>(2u*(p1%32u)))|((pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u]);
//  v2=(pstr2[p2/32u]>>(2u*(p2%32u)))|((pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u]);
//  tmpas1+=kmer2str(v1).substr(0,l);
//  tmpas2+=kmer2str(v2).substr(0,l);
// fend:
//  pas1=tmpas1+pas1;
//  pas2=tmpas2+pas2;
}



void seqident_seg_left_local(const eseq& s1,long p1,const eseq& s2,long p2,ealigndata& adata,const ealignscore& as){
  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[5];
  int k;
  int itmp;

  if (p1==0 || p2==0) {
    adata.s1=p1;
    adata.s2=p2;
    return;
  }

//  cout << "seg_left id: " << p1 << " :: " << p2 << " len1: " << p1 << " len2: " << p2 << endl;
  for (; p1>=32 && p2>=32; p1-=k,p2-=k){
    v1=pstr1[p1/32u-1u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u-1u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//    cout << p1 << " " << long2str(v1) << endl;
//    cout << p2 << " " << long2str(v2) << endl;
    for (k=IKMERSIZE; k<32; k+=IKMERSIZE,v1<<=IKMERSIZE*2u,v2<<=IKMERSIZE*2u){
//      cout << p1 << "+" << k << " " << kmer2str((v1>>((32u-IKMERSIZE)*2u))&IKMERMASK) << endl; 
//      cout << p2 << "+" << k << " " << kmer2str((v2>>((32u-IKMERSIZE)*2u))&IKMERMASK) << endl; 
      tmp[0]=seq_ident_table[((v1^v2)>>(64u-IKMERSIZE*2u))&IKMERMASK];
     
//      if (!galign && tmp[0]<IKMERSIZE/2){
      if (tmp[0]<IKMERSIZE/2){
        adata.s1=p1-k;
        adata.s2=p2-k;
        return;
      }
      adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
//      sumcounts(tncount,p1-k,seq_alignment_lt[((v1^v2)>>(64u-IKMERSIZE*2u))&IKMERMASK]);
    }
  }
//  cout << "seg_left rest id: " << p1 << " :: " << p2 << " len1: " << p1 << " len2: " << p2 << endl;
//  cout << "2nd last p1: " << p1 << " p2: " << p2 << endl;
  for (;p1>=int(IKMERSIZE) && p2>=int(IKMERSIZE); p1-=k,p2-=k){ // needed to handle gap shifts gracefully
//    cout << "- 2nd last p1: " << p1 << " p2: " << p2 << endl;
//    cout << p1 << endl;
//    cout << p2 << endl;
    v1=(p1>=32?pstr1[p1/32u-1u]>>(2u*(p1%32u)):0ul);
    v2=(p2>=32?pstr2[p2/32u-1u]>>(2u*(p2%32u)):0ul);
    v1|=(pstr1[p1/32u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//    cout << p1 << " " << long2str(v1) << endl;
//    cout << p2 << " " << long2str(v2) << endl;
    for (k=IKMERSIZE; p1>=k && p2>=k && k<=32; k+=IKMERSIZE,v1<<=IKMERSIZE*2u,v2<<=IKMERSIZE*2u){
//      cout << p1 << "+" << k << " " << kmer2str((v1>>((32u-IKMERSIZE)*2u))&IKMERMASK) << endl; 
//      cout << p2 << "+" << k << " " << kmer2str((v2>>((32u-IKMERSIZE)*2u))&IKMERMASK) << endl; 
      tmp[0]=seq_ident_table[((v1^v2)>>(64u-IKMERSIZE*2u))&IKMERMASK];
      
//      if (!galign && tmp[0]<=IKMERSIZE/2){
      if (tmp[0]<=IKMERSIZE/2){
        adata.s1=p1-k;
        adata.s2=p2-k;
//        cout << "break p1: " << adata.s1 << " p2: " << adata.s2 << endl;
        return;
      }
      adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
//      sumcounts(tncount,p1-k,seq_alignment_lt[((v1^v2)>>(64u-IKMERSIZE*2u))&IKMERMASK]);
    }
  }

  long l=(p1<p2?p1:p2);
//  ldieif(p1<0 || p2<0,"negative sequence positions, p1: " + estr(p1) +" p2: "+p2);
//  cout << "p1: " << p1 << " p2: " << p2 << " l: " << l << endl;
  if (l<=0) { 
    adata.s1=p1-l;
    adata.s2=p2-l;
    return;
  }

  p1+=IKMERSIZE-l; p2+=IKMERSIZE-l; // adjust p1 and p2, such that we get a full kmer which includes the left end, we will mask the right side afterwards

  v1=(p1>=32?pstr1[p1/32u-1u]>>(2u*(p1%32u)):0ul);
  v2=(p2>=32?pstr2[p2/32u-1u]>>(2u*(p2%32u)):0ul);
  v1|=(pstr1[p1/32u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
  v2|=(pstr2[p2/32u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;

  tmp[0]=seq_ident_table[((v1^v2)>>(62u-IKMERSIZE*2u)) & ((0x1ul<<(l*2ul))-1ul) & IKMERMASK];
//  if (!galign && tmp[0]<l/2){
  if (tmp[0]<l/2){
    adata.s1=p1-IKMERSIZE+l;
    adata.s2=p2-IKMERSIZE+l;
    return;
  }

//  cout << long2str(v1&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
//  cout << long2str(v2&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
  adata.matches+=tmp[0]-(IKMERSIZE-l);
  adata.mismatches+=IKMERSIZE-tmp[0];
  adata.s1=p1-IKMERSIZE;
  adata.s2=p2-IKMERSIZE;
//  sumcounts(tncount,p1-IKMERSIZE,seq_alignment_lt[((v1^v2)>>(64u-IKMERSIZE*2u)) & IKMERMASK]&((0x1ul<<(l*8ul))-1ul));
//  cout << "final p1: " << p1 << " p2: " << p2 << " l: " << l << endl;

/*
  tmp[0]=seq_ident_table[((v1^v2)>>((32u-IKMERSIZE)*2u))&~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul)&IKMERMASK]; // put mask at end of 64bit int, i.e.: FF000000

  cout << long2str((v1>>((32u-IKMERSIZE)*2u))&(~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul))&IKMERMASK) << " " << long2str((~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul))&IKMERMASK) << endl;
  cout << long2str((v2>>((32u-IKMERSIZE)*2u))&~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul)&IKMERMASK) << " " << long2str(~((0x1ul<<(IKMERSIZE*2u-l*2ul))-1ul)&IKMERMASK) << endl;
  matches+=tmp[0]-(IKMERSIZE-l);
  mismatches+=IKMERSIZE-tmp[0];
*/
}

void seqident_seg_local_right(const eseq& s1,long p1,long e1,const eseq& s2,long p2,long e2,ealigndata& adata,ealignscore& as){
  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[5];
  int k;
  int itmp;

//  cout << "id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1<e1-32 && p2<e2-32; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%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[0]=seq_ident_table[(v1^v2)&IKMERMASK];
//      if (!galign && tmp[0]<IKMERSIZE/2){
      if (tmp[0]<IKMERSIZE/2){
        adata.e1=p1+k;
        adata.e2=p2+k;
        return;
      }
      adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
//      sumcounts(tncount,p1+k,seq_alignment_lt[(v1^v2)&IKMERMASK]);
    }
  }

//  cout << "rest id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1+IKMERSIZE<e1 && p2+IKMERSIZE<e2; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//    cout << p1 << " " << long2str(v1) << endl;
//    cout << p2 << " " << long2str(v2) << endl;
    for (k=0; p1+k+IKMERSIZE<e1 && p2+k+IKMERSIZE<e2; k+=IKMERSIZE,v1>>=IKMERSIZE*2u,v2>>=IKMERSIZE*2u){
//      cout << p1 << "+" << k << " " << kmer2str(v1) << endl;
//      cout << p2 << "+" << k << " " << kmer2str(v2) << endl;
      tmp[0]=seq_ident_table[(v1^v2)&IKMERMASK];
     
//      if (!galign && tmp[0]<IKMERSIZE/2){
      if (tmp[0]<IKMERSIZE/2){
        adata.e1=p1+k;
        adata.e2=p2+k;
        return;
      }
      adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
//      sumcounts(tncount,p1+k,seq_alignment_lt[(v1^v2)&IKMERMASK]);
    }
  }
  long l=(e1-p1<e2-p2?e1-p1:e2-p2);
//  cout << "missing: " << p1 << " l: " << l << endl;
//  ldieif(l<0,"error past end of sequence");
  if (l<=0) {
    adata.e1=p1+l;
    adata.e2=p2+l;
    return;
  }
  v1=pstr1[p1/32u]>>(2u*(p1%32u));
  v2=pstr2[p2/32u]>>(2u*(p2%32u));
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;
  v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
  v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;
  tmp[0]=seq_ident_table[(v1^v2)&((0x1ul<<(l*2ul))-1ul)];
//  cout << long2str(v1&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
//  cout << long2str(v2&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;

//  if (!galign && tmp[0]<(IKMERSIZE-l)/2) {
  if (tmp[0]<(IKMERSIZE-l)/2) {
    adata.e1=p1;
    adata.e2=p2;
    return;
  }

  adata.matches+=tmp[0]-(IKMERSIZE-l);
  adata.mismatches+=IKMERSIZE-tmp[0];
//  sumcounts(tncount,p1,seq_alignment_lt[(v1^v2)&IKMERMASK]&((0x1ul<<(l*8u))-1ul));
  adata.e1=p1+l;
  adata.e2=p2+l;
}


void seqident_seg_nogaps(const eseq& s1,long p1,long e1,const eseq& s2,long p2,long e2,ealigndata& adata,const ealignscore& as){
  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;
  int k;
  int itmp;

//  cout << "id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1<e1-32 && p2<e2-32; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%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=seq_ident_table[(v1^v2)&IKMERMASK];
      adata.matches+=tmp;
      adata.mismatches+=IKMERSIZE-tmp;
      adata.profile.add(AT_MATCH,tmp);
      adata.profile.add(AT_MISS,IKMERSIZE-tmp);
      adata._score+=tmp*as.match-(IKMERSIZE-tmp)*as.mismatch;
    }
  }

  for (; p1+IKMERSIZE<e1 && p2+IKMERSIZE<e2; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
    for (k=0; p1+k+IKMERSIZE<e1 && p2+k+IKMERSIZE<e2; k+=IKMERSIZE,v1>>=IKMERSIZE*2u,v2>>=IKMERSIZE*2u){
      tmp=seq_ident_table[(v1^v2)&IKMERMASK];
      adata.matches+=tmp;
      adata.mismatches+=IKMERSIZE-tmp;
      adata.profile.add(AT_MATCH,tmp);
      adata.profile.add(AT_MISS,IKMERSIZE-tmp);
      adata._score+=tmp*as.match-(IKMERSIZE-tmp)*as.mismatch;
    }
  }

  long l=(e1-p1<e2-p2?e1-p1:e2-p2);
//  long ml=(e1-p1>e2-p2?e1-p1:e2-p2);
  if (l<=0) return;
  v1=pstr1[p1/32u]>>(2u*(p1%32u));
  v2=pstr2[p2/32u]>>(2u*(p2%32u));
  v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
  v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
  tmp=seq_ident_table[(v1^v2)&((0x1ul<<(l*2ul))-1ul)];
  adata.matches+=tmp-(IKMERSIZE-l);
  adata.mismatches+=IKMERSIZE-tmp;
  adata._score+=(tmp-(IKMERSIZE-l))*as.match-(IKMERSIZE-tmp)*as.mismatch;
  adata.profile.add(AT_MATCH,tmp-(IKMERSIZE-l));
  adata.profile.add(AT_MISS,IKMERSIZE-tmp);
}


void seqident_seg(const eseq& s1,long p1,long e1,const eseq& s2,long p2,long e2,ealigndata& adata,const ealignscore& as){
  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[5];
  int k;
  int itmp;

//  cout << "id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1<e1-32 && p2<e2-32; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%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[0]=seq_ident_table[(v1^v2)&IKMERMASK];
      tmp[1]=seq_ident_table[(v1^(v2>>2ul))&IKMERMASK];
      tmp[2]=seq_ident_table[(v1^(v2>>4ul))&IKMERMASK];
      tmp[3]=seq_ident_table[(v2^(v1>>2ul))&IKMERMASK];
      tmp[4]=seq_ident_table[(v2^(v1>>4ul))&IKMERMASK];
      
      itmp=0;
      for (int i=1; i<5; ++i){
        if (tmp[i]>tmp[itmp]) itmp=i;
      }
      if (itmp!=0 && tmp[itmp]>tmp[0]+IKMERSIZE/4){
        switch(itmp){
          case 1: p2+=1; adata.gaps+=1; break;
          case 2: p2+=2; adata.gaps+=2; break;
          case 3: p1+=1; adata.gaps+=1; break;
          case 4: p1+=2; adata.gaps+=2; break;
        }
        break;
      }else{
        adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
        adata._score+=tmp[0]*as.match-(IKMERSIZE-tmp[0])*as.mismatch;
      }
    }
  }

//  cout << "rest id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1+IKMERSIZE<e1 && p2+IKMERSIZE<e2; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//    cout << p1 << " " << long2str(v1) << endl;
//    cout << p2 << " " << long2str(v2) << endl;
    for (k=0; p1+k+IKMERSIZE<e1 && p2+k+IKMERSIZE<e2; k+=IKMERSIZE,v1>>=IKMERSIZE*2u,v2>>=IKMERSIZE*2u){
//      cout << p1 << "+" << k << " " << kmer2str(v1) << endl;
//      cout << p2 << "+" << k << " " << kmer2str(v2) << endl;
      tmp[0]=seq_ident_table[(v1^v2)&IKMERMASK];
      tmp[1]=seq_ident_table[(v1^(v2>>2ul))&IKMERMASK];
      tmp[2]=seq_ident_table[(v1^(v2>>4ul))&IKMERMASK];
      tmp[3]=seq_ident_table[(v2^(v1>>2ul))&IKMERMASK];
      tmp[4]=seq_ident_table[(v2^(v1>>4ul))&IKMERMASK];
      
      itmp=0;
      for (int i=1; i<5; ++i){
        if (tmp[i]>tmp[itmp]) itmp=i;
      }
      if (itmp!=0 && tmp[itmp]>tmp[0]+IKMERSIZE/4){
        switch(itmp){
          case 1: p2+=1; adata.gaps+=1; break;
          case 2: p2+=2; adata.gaps+=2; break;
          case 3: p1+=1; adata.gaps+=1; break;
          case 4: p1+=2; adata.gaps+=2; break;
        }
        break;
      }else{
        adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
        adata._score+=tmp[0]*as.match-(IKMERSIZE-tmp[0])*as.mismatch;
      }
    }
  }
  long l=(e1-p1<e2-p2?e1-p1:e2-p2);
  long ml=(e1-p1>e2-p2?e1-p1:e2-p2);
  if (e1!=s1.seqlen && e2!=s2.seqlen)
    adata.gaps+=ml-l;
//  cout << "missing: " << p1 << " l: " << l << endl;
  if (l<=0) return;
  v1=pstr1[p1/32u]>>(2u*(p1%32u));
  v2=pstr2[p2/32u]>>(2u*(p2%32u));
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;
  v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
  v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;
  tmp[0]=seq_ident_table[(v1^v2)&((0x1ul<<(l*2ul))-1ul)];
//  cout << long2str(v1&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
//  cout << long2str(v2&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
  adata.matches+=tmp[0]-(IKMERSIZE-l);
  adata.mismatches+=IKMERSIZE-tmp[0];
  adata._score+=(tmp[0]-(IKMERSIZE-l))*as.match-(IKMERSIZE-tmp[0])*as.mismatch;
}



void seqident_seg(const eseq& s1,long p1,long e1,const eseq& s2,long p2,long e2,ealigndata& adata,unsigned char *tncount){
  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[5];
  int k;
  int itmp;

//  cout << "id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1<e1-32 && p2<e2-32; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%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[0]=seq_ident_table[(v1^v2)&IKMERMASK];
      tmp[1]=seq_ident_table[(v1^(v2>>2ul))&IKMERMASK];
      tmp[2]=seq_ident_table[(v1^(v2>>4ul))&IKMERMASK];
      tmp[3]=seq_ident_table[(v2^(v1>>2ul))&IKMERMASK];
      tmp[4]=seq_ident_table[(v2^(v1>>4ul))&IKMERMASK];
      
      itmp=0;
      for (int i=1; i<5; ++i){
        if (tmp[i]>tmp[itmp]) itmp=i;
      }
      if (itmp!=0 && tmp[itmp]>tmp[0]+IKMERSIZE/4){
        switch(itmp){
          case 1: p2+=1; adata.gaps+=1; break;
          case 2: p2+=2; adata.gaps+=2; break;
          case 3: p1+=1; adata.gaps+=1; break;
          case 4: p1+=2; adata.gaps+=2; break;
        }
        break;
      }else{
        adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
        sumcounts(tncount,p1+k,seq_alignment_lt[(v1^v2)&IKMERMASK]);
      }
    }
  }

//  cout << "rest id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1+IKMERSIZE<e1 && p2+IKMERSIZE<e2; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//    cout << p1 << " " << long2str(v1) << endl;
//    cout << p2 << " " << long2str(v2) << endl;
    for (k=0; p1+k+IKMERSIZE<e1 && p2+k+IKMERSIZE<e2; k+=IKMERSIZE,v1>>=IKMERSIZE*2u,v2>>=IKMERSIZE*2u){
//      cout << p1 << "+" << k << " " << kmer2str(v1) << endl;
//      cout << p2 << "+" << k << " " << kmer2str(v2) << endl;
      tmp[0]=seq_ident_table[(v1^v2)&IKMERMASK];
      tmp[1]=seq_ident_table[(v1^(v2>>2ul))&IKMERMASK];
      tmp[2]=seq_ident_table[(v1^(v2>>4ul))&IKMERMASK];
      tmp[3]=seq_ident_table[(v2^(v1>>2ul))&IKMERMASK];
      tmp[4]=seq_ident_table[(v2^(v1>>4ul))&IKMERMASK];
      
      itmp=0;
      for (int i=1; i<5; ++i){
        if (tmp[i]>tmp[itmp]) itmp=i;
      }
      if (itmp!=0 && tmp[itmp]>tmp[0]+IKMERSIZE/4){
        switch(itmp){
          case 1: p2+=1; adata.gaps+=1; break;
          case 2: p2+=2; adata.gaps+=2; break;
          case 3: p1+=1; adata.gaps+=1; break;
          case 4: p1+=2; adata.gaps+=2; break;
        }
        break;
      }else{
        adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
        sumcounts(tncount,p1+k,seq_alignment_lt[(v1^v2)&IKMERMASK]);
      }
    }
  }
  long l=(e1-p1<e2-p2?e1-p1:e2-p2);
  long ml=(e1-p1>e2-p2?e1-p1:e2-p2);
  if (e1!=s1.seqlen && e2!=s2.seqlen)
    adata.gaps+=ml-l;
//  cout << "missing: " << p1 << " l: " << l << endl;
  if (l<=0) return;
  v1=pstr1[p1/32u]>>(2u*(p1%32u));
  v2=pstr2[p2/32u]>>(2u*(p2%32u));
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;
  v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
  v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;
  tmp[0]=seq_ident_table[(v1^v2)&((0x1ul<<(l*2ul))-1ul)];
//  cout << long2str(v1&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
//  cout << long2str(v2&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
  adata.matches+=tmp[0]-(IKMERSIZE-l);
  adata.mismatches+=IKMERSIZE-tmp[0];
  sumcounts(tncount,p1,seq_alignment_lt[(v1^v2)&IKMERMASK]&((0x1ul<<(l*8u))-1ul));
}

void seqident_rev_seg(const eseq& s1,long p1,long e1,const eseq& s2,long p2,long e2,ealigndata& adata){
  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[5];
  int k;
  int itmp;

//  cout << "id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1<e1-32 && p2<e2-32; p1+=k,p2+=k){
    long rp2=s2.seqlen-p2-32;
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[rp2/32u]>>(2u*(rp2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[rp2/32u+1u]<<(64u-2u*(rp2%32u)))&safe_shift[rp2%32u];
    for (k=0; k<32u-IKMERSIZE; k+=IKMERSIZE,v1>>=IKMERSIZE*2u,v2<<=IKMERSIZE*2u){
      tmp[0]=seq_ident_table[(v1^kmer_rev_lt[(v2>>((32u-IKMERSIZE)*2u))&IKMERMASK])&IKMERMASK];
      tmp[1]=seq_ident_table[(v1^kmer_rev_lt[(v2>>((32u-IKMERSIZE)*2u-2u))&IKMERMASK])&IKMERMASK];
      tmp[2]=seq_ident_table[(v1^kmer_rev_lt[(v2>>((32u-IKMERSIZE)*2u-4u))&IKMERMASK])&IKMERMASK];
      tmp[3]=seq_ident_table[(kmer_rev_lt[(v2>>((32u-IKMERSIZE)*2u))&IKMERMASK]^(v1>>2ul))&IKMERMASK];
      tmp[4]=seq_ident_table[(kmer_rev_lt[(v2>>((32u-IKMERSIZE)*2u))&IKMERMASK]^(v1>>4ul))&IKMERMASK];
      
      itmp=0;
      for (int i=1; i<5; ++i){
        if (tmp[i]>tmp[itmp]) itmp=i;
      }
      if (itmp!=0 && tmp[itmp]>tmp[0]+IKMERSIZE/4){
        switch(itmp){
          case 1: p2+=1; adata.gaps+=1; break;
          case 2: p2+=2; adata.gaps+=2; break;
          case 3: p1+=1; adata.gaps+=1; break;
          case 4: p1+=2; adata.gaps+=2; break;
        }
        break;
      }else{
        adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
      }
    }
  }

//  cout << "rest id: " << p1 << "," << e1 << " :: " << p2 << "," <<e2 << " len1: " << e1-p1 << " len2: " << e2-p2 << endl;
  for (; p1+IKMERSIZE<e1 && p2+IKMERSIZE<e2; p1+=k,p2+=k){
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//    cout << p1 << " " << long2str(v1) << endl;
//    cout << p2 << " " << long2str(v2) << endl;
    for (k=0; p1+k+IKMERSIZE<e1 && p2+k+IKMERSIZE<e2; k+=IKMERSIZE,v1>>=IKMERSIZE*2u,v2>>=IKMERSIZE*2u){
//      cout << p1 << "+" << k << " " << kmer2str(v1) << endl;
//      cout << p2 << "+" << k << " " << kmer2str(v2) << endl;
      tmp[0]=seq_ident_table[(v1^v2)&IKMERMASK];
      tmp[1]=seq_ident_table[(v1^(v2>>2ul))&IKMERMASK];
      tmp[2]=seq_ident_table[(v1^(v2>>4ul))&IKMERMASK];
      tmp[3]=seq_ident_table[(v2^(v1>>2ul))&IKMERMASK];
      tmp[4]=seq_ident_table[(v2^(v1>>4ul))&IKMERMASK];
      
      itmp=0;
      for (int i=1; i<5; ++i){
        if (tmp[i]>tmp[itmp]) itmp=i;
      }
      if (itmp!=0 && tmp[itmp]>tmp[0]+IKMERSIZE/4){
        switch(itmp){
          case 1: p2+=1; adata.gaps+=1; break;
          case 2: p2+=2; adata.gaps+=2; break;
          case 3: p1+=1; adata.gaps+=1; break;
          case 4: p1+=2; adata.gaps+=2; break;
        }
        break;
      }else{
        adata.matches+=tmp[0]; adata.mismatches+=IKMERSIZE-tmp[0];
      }
    }
  }
  long l=(e1-p1<e2-p2?e1-p1:e2-p2);
//  cout << "missing: " << p1 << " l: " << l << endl;
  if (l<=0) return;
  v1=pstr1[p1/32u]>>(2u*(p1%32u));
  v2=pstr2[p2/32u]>>(2u*(p2%32u));
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;
  v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
  v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
//  cout << p1 << " " << long2str(v1) << endl;
//  cout << p2 << " " << long2str(v2) << endl;
  tmp[0]=seq_ident_table[(v1^v2)&((0x1ul<<(l*2ul))-1ul)];
//  cout << long2str(v1&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
//  cout << long2str(v2&((0x1ul<<(l*2u))-1ul)) << " " << long2str(((0x1ul<<(l*2u))-1ul)) << endl;
  adata.matches+=tmp[0]-(IKMERSIZE-l);
  adata.mismatches+=IKMERSIZE-tmp[0];
}



void seqident_rev(const eseq& s1,long offset1,eintarray& kmerpos1,const eseq& s2,long offset2,eintarray& kmerpos2,ealigndata& adata)
{
  long p1,p2;
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s1.seq._str);
  unsigned long *pstr2=reinterpret_cast<unsigned long*>(s2.seq._str);
  unsigned long v1,v2;
  long ndelta=-s1.seqlen-s2.seqlen;
  long lastkmerlen,lastkmerpos,lastndelta;

//  for (int i=0; i<11; ++i)
//    match[i]=0;
  adata.matches=0; adata.mismatches=0; adata.gaps=0;

  char tmp[5];
  char itmp;
  int k;

  ebasicarray<ediag> diags;

  if (s1.seqlen<s2.seqlen){
    lastkmerpos=-2*int(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1.seqlen-s2.seqlen;
    for (p1=0; p1<s1.seqlen-32; p1+=k){
      v1=pstr1[p1/32u]>>(2u*(p1%32u));
      v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
      for (k=0; k<32u-KMERSIZE; ++k,v1>>=2u){
        int kpos2=kmerpos2[kmer_rev_lt[v1&KMERMASK]]-offset2;
        if (kpos2>-1){
          kpos2=s2.seqlen-kpos2;
          int d=p1+k-lastkmerpos;
          if (d<=KMERSIZE && ((lastndelta==kpos2-p1-k) || (v1&KMERMASK)==seqrevkmer(s2,p1+k+lastndelta))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p1: " << lastkmerpos-lastkmerlen+KMERSIZE<< "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=kpos2-p1-k;
          }
          lastkmerpos=p1+k;
        }
      }
    }
    if (lastkmerlen){
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
//      cout << "end p1: " << lastkmerpos-lastkmerlen+KMERSIZE << "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
    }
  }
/*
  else{
    lastkmerpos=-2*int(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1.seqlen-s2.seqlen;
    for (p2=0; p2<s2.seqlen-32 && ndelta==-s1.seqlen-s2.seqlen; p2+=k){
      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-KMERSIZE; ++k,v2>>=2u){
        int kpos1=kmerpos1[kmer_rev_lt[v2&KMERMASK]]-offset1;
        if (kpos1>-1){
          int d=p2+k-lastkmerpos;
          if (d<2*KMERSIZE && ((lastndelta==kpos1-p2-k) || kmer_rev_lt[v2&KMERMASK]==seqkmer(s1,p2+k-lastndelta))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=p2+k-kpos1;
          }
          lastkmerpos=p2+k;
        }
      }
    }
    if (lastkmerlen){
//      cout << "end p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
    }
  }
*/

  // sparse dynamic programming to find longest chain of segments
  multimap<int,ediag*> segments;
  for (int j=0; j<diags.size(); ++j){
    ediag& diag(diags.at(j));
    diag.V=diag.j-diag.i;
    segments.insert(pair<int,ediag*>(diag.i,&diag));
    segments.insert(pair<int,ediag*>(diag.j,&diag));
  }

  multimap<int,ediag*> Lsegments;
  multimap<int,ediag*>::iterator sit;
  multimap<int,ediag*>::iterator slb;
  ediag *bestseg=0x00;
  for (sit=segments.begin(); sit!=segments.end(); ++sit){
//      cout << "x: " << sit->first << " start: " << (sit->first==sit->second->i?"1":"0") << " i: " << sit->second->i << " i2: " << sit->second->i2 << " Lsegs: " << Lsegments.size() << endl;
    if (sit->first==sit->second->i) { // start of rectangle
      slb=Lsegments.upper_bound(sit->second->i2-1);
      if (Lsegments.size()==0 || slb==Lsegments.begin()) continue;
      --slb;
      sit->second->V=sit->second->j-sit->second->i+slb->second->V;
      sit->second->bestseg=slb->second;
//        cout << "linking to: " << slb->second->i << endl;
    } else {  // end of rectangle
      slb=Lsegments.lower_bound(sit->second->j2);
      if (slb==Lsegments.end() || sit->second->V > slb->second->V){
        while (slb!=Lsegments.end()){
          if (slb->second->j2>=sit->second->j2 && slb->second->V<=sit->second->V){
            Lsegments.erase(slb++);
            continue;
          }
          ++slb;
        }
        if (bestseg==0x00 || sit->second->V > bestseg->V)
          bestseg=sit->second;
        Lsegments.insert(pair<int,ediag*>(sit->second->j2,sit->second));
      }
    }
  }

//  seqident_seg_right(s1,s2,bestseg,matches,mismatches,gaps);
  if (bestseg==0x00) return;

//  cout << "s1: " << s1.seqlen << " s2: " << s2.seqlen << endl;

  seqident_rev_seg(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata);
//  cout << "rightid: " << bestseg->j << "," << s1.seqlen << " l: " << s1.seqlen-bestseg->j << " d: " << bestseg->i2-bestseg->i << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
  while (bestseg!=0x00 && bestseg->bestseg != 0x00){
    adata.matches+=bestseg->j-bestseg->i;
//    cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
    seqident_rev_seg(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata);
//    cout << "segid: " << bestseg->i << "," << bestseg->j << " l: " << bestseg->j-bestseg->i << " d: " << bestseg->i2-bestseg->i << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
    bestseg=bestseg->bestseg;
  }
  adata.matches+=bestseg->j-bestseg->i;
//  cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
//  seqident_rev_seg_left(s1,bestseg->i,s2,bestseg->i2,matches,mismatches,gaps);
//  cout << "leftid: "<< bestseg->i << " l: " << (bestseg->i<bestseg->i2?bestseg->i:bestseg->i2) << " d: " << bestseg->i2-bestseg->i << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;

  return;
}

/*
void seqident(const eseq& s1,int offset1,eintarray& kmerpos1,const eseq& s2,int offset2,eintarray& kmerpos2,ealigndata& adata,unsigned char *tncounts)
{
  int p1,p2;
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s1.seq._str);
  unsigned long *pstr2=reinterpret_cast<unsigned long*>(s2.seq._str);
  unsigned long v1,v2;
  int ndelta=-s1.seqlen-s2.seqlen;
  int lastkmerlen,lastkmerpos,lastndelta;

//  for (int i=0; i<11; ++i)
//    match[i]=0;
  adata.matches=0; adata.mismatches=0; adata.gaps=0;
  adata.s1=-1; adata.e1=-1;
  adata.s2=-1; adata.e2=-1;

  char tmp[5];
  char itmp;
  int k;

  ebasicarray<ediag> diags;

  if (s1.seqlen<s2.seqlen){
    lastkmerpos=-2*int(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1.seqlen-s2.seqlen;
    for (p1=0; p1<s1.seqlen-32; p1+=k){
      v1=pstr1[p1/32u]>>(2u*(p1%32u));
      v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
      for (k=0; k<32u-KMERSIZE; ++k,v1>>=2u){
        int kpos2=kmerpos2[v1&KMERMASK]-offset2;
        if (kpos2>-1){
          int d=p1+k-lastkmerpos;
          if (d<=KMERSIZE && ((lastndelta==kpos2-p1-k) || (v1&KMERMASK)==seqkmer(s2,p1+k+lastndelta))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p1: " << lastkmerpos-lastkmerlen+KMERSIZE<< "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=kpos2-p1-k;
          }
          lastkmerpos=p1+k;
        }
      }
    }
    if (lastkmerlen){
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
//      cout << "end p1: " << lastkmerpos-lastkmerlen+KMERSIZE << "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
    }
  }else{
    lastkmerpos=-2*int(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1.seqlen-s2.seqlen;
    for (p2=0; p2<s2.seqlen-32 && ndelta==-s1.seqlen-s2.seqlen; p2+=k){
      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-KMERSIZE; ++k,v2>>=2u){
        int kpos1=kmerpos1[v2&KMERMASK]-offset1;
        if (kpos1>-1){
          int d=p2+k-lastkmerpos;
          if (d<=KMERSIZE && ((lastndelta==kpos1-p2-k) || (v2&KMERMASK)==seqkmer(s1,p2+k-lastndelta))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=p2+k-kpos1;
          }
          lastkmerpos=p2+k;
        }
      }
    }
    if (lastkmerlen){
//      cout << "end p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
    }
  }

  // sparse dynamic programming to find longest chain of segments
  multimap<int,ediag*> segments;
  for (int j=0; j<diags.size(); ++j){
    ediag& diag(diags.at(j));
    diag.V=diag.j-diag.i;
    segments.insert(pair<int,ediag*>(diag.i,&diag));
    segments.insert(pair<int,ediag*>(diag.j,&diag));
  }

  multimap<int,ediag*> Lsegments;
  multimap<int,ediag*>::iterator sit;
  multimap<int,ediag*>::iterator slb;
  ediag *bestseg=0x00;
  for (sit=segments.begin(); sit!=segments.end(); ++sit){
//      cout << "x: " << sit->first << " start: " << (sit->first==sit->second->i?"1":"0") << " i: " << sit->second->i << " i2: " << sit->second->i2 << " Lsegs: " << Lsegments.size() << endl;
    if (sit->first==sit->second->i) { // start of rectangle
      slb=Lsegments.upper_bound(sit->second->i2-1);
      if (Lsegments.size()==0 || slb==Lsegments.begin()) continue;
      --slb;
      sit->second->V=sit->second->j-sit->second->i+slb->second->V;
      sit->second->bestseg=slb->second;
//        cout << "linking to: " << slb->second->i << endl;
    } else {  // end of rectangle
      slb=Lsegments.lower_bound(sit->second->j2);
      if (slb==Lsegments.end() || sit->second->V > slb->second->V){
        while (slb!=Lsegments.end()){
          if (slb->second->j2>=sit->second->j2 && slb->second->V<=sit->second->V){
            Lsegments.erase(slb++);
            continue;
          }
          ++slb;
        }
        if (bestseg==0x00 || sit->second->V > bestseg->V)
          bestseg=sit->second;
        Lsegments.insert(pair<int,ediag*>(sit->second->j2,sit->second));
      }
    }
  }

//  seqident_seg_right(s1,s2,bestseg,matches,mismatches,gaps);
  if (bestseg==0x00) return;

//  cout << "s1: " << s1.seqlen << " s2: " << s2.seqlen << endl;

  int l=(s2.seqlen-bestseg->j2<s1.seqlen-bestseg->j?s2.seqlen-bestseg->j2:s1.seqlen-bestseg->j);
  adata.e1=bestseg->j+l;
  adata.e2=bestseg->j2+l;

  estr ps1,ps2;

  seqident_seg(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,tncounts);
//  cout << "rightid: " << bestseg->j << "," << s1.seqlen << " l: " << s1.seqlen-bestseg->j << " d: " << bestseg->i2-bestseg->i << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
  while (bestseg!=0x00 && bestseg->bestseg != 0x00){
    adata.matches+=bestseg->j-bestseg->i;
    seqncount(s1,bestseg->i,bestseg->j,s2,bestseg->i2,bestseg->j2,tncounts);
//    cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
    seqident_seg(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata,tncounts);
//    cout << "segid: " << bestseg->i << "," << bestseg->j << " l: " << bestseg->j-bestseg->i << " d: " << bestseg->i2-bestseg->i << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
    bestseg=bestseg->bestseg;
  }
  adata.matches+=bestseg->j-bestseg->i;
  seqncount(s1,bestseg->i,bestseg->j,s2,bestseg->i2,bestseg->j2,tncounts);
//  cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
  seqident_seg_left(s1,bestseg->i,s2,bestseg->i2,adata);
//  cout << "leftid: "<< bestseg->i << " l: " << (bestseg->i<bestseg->i2?bestseg->i:bestseg->i2) << " d: " << bestseg->i2-bestseg->i << " m: " << matches << " miss: " << mismatches << " gaps: " << gaps << endl;
  l=(bestseg->i2<bestseg->i?bestseg->i2:bestseg->i);
  adata.s1=bestseg->i-l;
  adata.s2=bestseg->i2-l;

  return;
}
*/

void seqident_local_fast(const eseq& s1,eintarray& kmerpos1,const eseq& s2,ealigndata& adata,esearchws& sws,const ealignscore& as)
{
  long p1,p2;
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s1.seq._str);
  unsigned long *pstr2=reinterpret_cast<unsigned long*>(s2.seq._str);
  unsigned long v1,v2;
  long lastkmerlen,lastkmerpos,lastndelta;

//  for (int i=0; i<11; ++i)
//    match[i]=0;
  adata.matches=0; adata.mismatches=0; adata.gaps=0; adata._score=0.0;
  adata.s1=-1; adata.e1=-1;
  adata.s2=-1; adata.e2=-1;

  char tmp[5];
  char itmp;
  int k;

  ebasicarray<ediag> diags;
  ebasicarray<ediag> dheap;

  if (s1.seqlen<s2.seqlen){
    lastkmerpos=-2*long(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1.seqlen-s2.seqlen;
    for (p1=0; p1<s1.seqlen-32; p1+=k){
      v1=pstr1[p1/32u]>>(2u*(p1%32u));
      v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
      for (k=0; k<32u-KMERSIZE; ++k,v1>>=2u){
        long kpos2=sws.kmerpos2[v1&KMERMASK]-sws.offset2;
        if (kpos2>-1){
          long d=p1+k-lastkmerpos;
          
          if (d<=KMERSIZE && ((lastndelta==kpos2-p1-k) || (p1+k+lastndelta+KMERSIZE<=s2.seqlen && (v1&KMERMASK)==seqkmer(s2,p1+k+lastndelta)))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p1: " << lastkmerpos-lastkmerlen+KMERSIZE<< "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=kpos2-p1-k;
          }
          lastkmerpos=p1+k;
        }
      }
    }
    if (lastkmerlen){
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
//      cout << "end p1: " << lastkmerpos-lastkmerlen+KMERSIZE << "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
    }
  }else{
    lastkmerpos=-2*long(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1.seqlen-s2.seqlen;
    for (p2=0; p2<s2.seqlen-32; p2+=k){
      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-KMERSIZE; ++k,v2>>=2u){
        long kpos1=kmerpos1[v2&KMERMASK]-sws.offset;
        if (kpos1>-1){
          long d=p2+k-lastkmerpos;
          if (d<=KMERSIZE && ((lastndelta==kpos1-p2-k) || (p2+k-lastndelta+KMERSIZE<=s1.seqlen && (v2&KMERMASK)==seqkmer(s1,p2+k-lastndelta)))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=p2+k-kpos1;
          }
          lastkmerpos=p2+k;
        }
      }
    }
    if (lastkmerlen){
//      cout << "end p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
    }
  }

  // sparse dynamic programming to find longest chain of segments
  multimap<long,ediag*> segments;
  for (long j=0; j<diags.size(); ++j){
    ediag& diag(diags.at(j));
    diag.V=(diag.j-diag.i)*as.match;
    segments.insert(pair<long,ediag*>(diag.i,&diag));
    segments.insert(pair<long,ediag*>(diag.j,&diag));
  }

  multimap<long,ediag*> Lsegments;
  multimap<long,ediag*>::iterator sit;
  multimap<long,ediag*>::iterator slb;
  ediag *bestseg=0x00;
  for (sit=segments.begin(); sit!=segments.end(); ++sit){
//      cout << "x: " << sit->first << " start: " << (sit->first==sit->second->i?"1":"0") << " i: " << sit->second->i << " i2: " << sit->second->i2 << " Lsegs: " << Lsegments.size() << endl;
    if (sit->first==sit->second->i) { // start of rectangle
      slb=Lsegments.upper_bound(sit->second->i2-1);
      if (Lsegments.size()==0 || slb==Lsegments.begin()) continue;
      --slb;
//      sit->second->V=sit->second->j-sit->second->i+slb->second->V;
      int gapdiff=abs((sit->second->i-sit->second->i2)-(slb->second->i-slb->second->i2));
      sit->second->V=(sit->second->j-sit->second->i)*as.match+gapdiff*as.gapext+(gapdiff>0?as.gapopen:0)+ MIN(sit->second->i-slb->second->j,sit->second->i2-slb->second->j2)*(0.75*as.match+0.25*as.mismatch) +slb->second->V;
      sit->second->bestseg=slb->second;
//      cout << "linking to: " << sit->second->i << " " << slb->second->j << endl;
    } else {  // end of rectangle
      slb=Lsegments.lower_bound(sit->second->j2);
      if (slb==Lsegments.end() || sit->second->V > slb->second->V){
        while (slb!=Lsegments.end()){
          if (slb->second->j2>=sit->second->j2 && slb->second->V<=sit->second->V){
            Lsegments.erase(slb++);
            continue;
          }
          ++slb;
        }
        if (bestseg==0x00 || sit->second->V > bestseg->V)
          bestseg=sit->second;
        Lsegments.insert(pair<long,ediag*>(sit->second->j2,sit->second));
      }
    }
  }

  if (bestseg==0x00) return;

//  estr fas1,fas2;
//  cout << "full alignment" << endl;
//  cout << "ali: " << 0 <<","<<s2.seqlen<<" : " << 0 <<","<<s1.seqlen << endl; 
//  adata._score=seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,fas1,fas2);
//  print_seqali(fas1,fas2);

//  estr pas1,pas2;
//  seqident_seg_local(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,tncounts);
//  seqident_seg_sw_right(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,tncounts);
//  estr as1,as2;
//  seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,as1,as2);
//  print_seqali(as1,as2);

//  seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,adata,alignws);
//  return;

  ldieif(s2.seqlen<bestseg->j2,"segment start larger than sequence length: "+estr(s2.seqlen)+" < "+estr(bestseg->j2));
  long minright=MIN(s1.seqlen-bestseg->j,s2.seqlen-bestseg->j2);
//  seqcalign_global_norightedgegap(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,alignws);

//  seqcalign_global_norightedgegap(s1,bestseg->j,MIN(s1.seqlen,bestseg->j+2*minright),s2,bestseg->j2,MIN(s2.seqlen,bestseg->j2+2*minright),adata,sws.alignws,as);
  seqident_seg(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,as);
  LDEBUG(D_SEQIDENT,cout << "rightid: " << bestseg->j << "," << s1.seqlen << " l: " << s1.seqlen-bestseg->j << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
  while (bestseg!=0x00 && bestseg->bestseg != 0x00){
    seqncount(s1,bestseg->i,bestseg->j,s2,bestseg->i2,bestseg->j2,adata,as);
    LDEBUG(D_SEQIDENT,cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
    LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
    long gapdiff=abs((bestseg->bestseg->i-bestseg->bestseg->i2)-(bestseg->i-bestseg->i2));
    if (gapdiff==0) /* fast alignment no gaps */
      seqident_seg_nogaps(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata,as);
    else /* full sw alignment */
      seqident_seg(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata,as);
//      seqcalign_global(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata,sws.alignws,as);
    LDEBUG(D_SEQIDENT,cout << "segid: " << bestseg->bestseg->j << "," << bestseg->i << " l: " << bestseg->i-bestseg->bestseg->j << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
    LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
//    cout << "segid: " << bestseg->i << "," << bestseg->j << " l: " << bestseg->j-bestseg->i << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl;
    bestseg=bestseg->bestseg;
  }
  seqncount(s1,bestseg->i,bestseg->j,s2,bestseg->i2,bestseg->j2,adata,as);
  LDEBUG(D_SEQIDENT,cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
//  seqident_seg_left_local(s1,bestseg->i,s2,bestseg->i2,adata,tncounts);
  long minleft=MIN(bestseg->i,bestseg->i2);
  seqident_seg_left(s1,bestseg->i,s2,bestseg->i2,adata,as);
//  seqcalign_global_noleftedgegap(s1,MAX(0,bestseg->i-2*minleft),bestseg->i,s2,MAX(0,bestseg->i2-2*minleft),bestseg->i2,adata,sws.alignws,as);

//  estr as1,as2;
//  seqcalign_global_noleftedgegap(s1,0,bestseg->i,s2,0,bestseg->i2,as1,as2);
//  cout << endl;
//  cout << as1 << endl << as2 << endl;
  LDEBUG(D_SEQIDENT,cout << "leftid: "<< bestseg->i << " l: " << (bestseg->i<bestseg->i2?bestseg->i:bestseg->i2) << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
//  cout << "partial alignment" << endl;
//  print_seqali(pas1,pas2);

  return;
}

void seqident_global(const eseq& s1,eintarray& kmerpos1,const eseq& s2,ealigndata& adata,esearchws& sws,const ealignscore& as)
{
  long p1,p2;
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s1.seq._str);
  unsigned long *pstr2=reinterpret_cast<unsigned long*>(s2.seq._str);
  unsigned long v1,v2;
  long lastkmerlen,lastkmerpos,lastndelta;

//  for (int i=0; i<11; ++i)
//    match[i]=0;
  adata.matches=0; adata.mismatches=0; adata.gaps=0; adata._score=0.0;
  adata.s1=-1; adata.e1=-1;
  adata.s2=-1; adata.e2=-1;

  char tmp[5];
  char itmp;
  int k;

  ebasicarray<ediag> diags;
  ebasicarray<ediag> dheap;

  if (s1.seqlen<s2.seqlen){
    lastkmerpos=-2*long(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1.seqlen-s2.seqlen;
    for (p1=0; p1<s1.seqlen-32; p1+=k){
      v1=pstr1[p1/32u]>>(2u*(p1%32u));
      v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
      for (k=0; k<32u-KMERSIZE; ++k,v1>>=2u){
        long kpos2=sws.kmerpos2[v1&KMERMASK]-sws.offset2;
        if (kpos2>-1){
          long d=p1+k-lastkmerpos;
          
          if (d<=KMERSIZE && ((lastndelta==kpos2-p1-k) || (p1+k+lastndelta+KMERSIZE<=s2.seqlen && (v1&KMERMASK)==seqkmer(s2,p1+k+lastndelta)))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p1: " << lastkmerpos-lastkmerlen+KMERSIZE<< "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=kpos2-p1-k;
          }
          lastkmerpos=p1+k;
        }
      }
    }
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    for (k=0; p1+k+KMERSIZE<s1.seqlen; ++k,v1>>=2u){
      long kpos2=sws.kmerpos2[v1&KMERMASK]-sws.offset2;
      if (kpos2>-1){
        long d=p1+k-lastkmerpos;
        if (d<=KMERSIZE && ((lastndelta==kpos2-p1-k) || (p1+k+lastndelta+KMERSIZE<=s2.seqlen && (v1&KMERMASK)==seqkmer(s2,p1+k+lastndelta)))){
          lastkmerlen+=d;
        }else{
          if (lastkmerlen){
            diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
          }
          lastkmerlen=KMERSIZE;
          lastndelta=kpos2-p1-k;
        }
        lastkmerpos=p1+k;
      }
    }
    if (lastkmerlen){
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
//      cout << "end p1: " << lastkmerpos-lastkmerlen+KMERSIZE << "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
    }
  }else{
    lastkmerpos=-2*long(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1.seqlen-s2.seqlen;
    for (p2=0; p2<s2.seqlen-32; p2+=k){
      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-KMERSIZE; ++k,v2>>=2u){
        long kpos1=kmerpos1[v2&KMERMASK]-sws.offset;
        if (kpos1>-1){
          long d=p2+k-lastkmerpos;
          if (d<=KMERSIZE && ((lastndelta==kpos1-p2-k) || (p2+k-lastndelta+KMERSIZE<=s1.seqlen && (v2&KMERMASK)==seqkmer(s1,p2+k-lastndelta)))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=p2+k-kpos1;
          }
          lastkmerpos=p2+k;
        }
      }
    }
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
    for (k=0; p2+k+KMERSIZE<s2.seqlen; ++k,v2>>=2u){
      long kpos1=kmerpos1[v2&KMERMASK]-sws.offset;
      if (kpos1>-1){
        long d=p2+k-lastkmerpos;
        if (d<=KMERSIZE && ((lastndelta==kpos1-p2-k) || (p2+k-lastndelta+KMERSIZE<=s1.seqlen && (v2&KMERMASK)==seqkmer(s1,p2+k-lastndelta)))){
          lastkmerlen+=d;
        }else{
          if (lastkmerlen){
            diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
          }
          lastkmerlen=KMERSIZE;
          lastndelta=p2+k-kpos1;
        }
        lastkmerpos=p2+k;
      }
    }
    if (lastkmerlen){
//      cout << "end p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
    }
  }

  if (diags.size()==0){
//    lderror("no shared segments found");
    return;
  }

  // sparse dynamic programming to find longest chain of segments
  ediag *bestseg=&diags[0];
  multimap<long,ediag*> segments;
  for (long j=0; j<diags.size(); ++j){
    ediag& diag(diags.at(j));
    diag.V=(diag.j-diag.i)*as.match;
    if (diag.V > bestseg->V) bestseg=&diag;
    segments.insert(pair<long,ediag*>(diag.i,&diag));
    segments.insert(pair<long,ediag*>(diag.j,&diag));
  }

  multimap<long,ediag*> Lsegments;
  multimap<long,ediag*>::iterator sit;
  multimap<long,ediag*>::iterator slb;
  for (sit=segments.begin(); sit!=segments.end(); ++sit){
//      cout << "x: " << sit->first << " start: " << (sit->first==sit->second->i?"1":"0") << " i: " << sit->second->i << " i2: " << sit->second->i2 << " Lsegs: " << Lsegments.size() << endl;
    if (sit->first==sit->second->i) { // start of rectangle
//      if (Lsegments.size()==0) continue;
//      slb=Lsegments.upper_bound(sit->second->i2-1);
//      if (slb==Lsegments.end()) continue;
      ediag *tmpbest=0x00;
      double tmpbestscore=(sit->second->j-sit->second->i)*as.match;
      for (slb=Lsegments.begin(); slb!=Lsegments.end(); ++slb){
        if (slb->first >= sit->second->i2) break;
//      sit->second->V=sit->second->j-sit->second->i+slb->second->V;
        long gapdiff=abs((sit->second->i-sit->second->i2)-(slb->second->i-slb->second->i2));
        double tmpscore=(sit->second->j-sit->second->i)*as.match - gapdiff*as.gapext-(gapdiff>0?as.gapopen:0) + slb->second->V;
        if (tmpscore>tmpbestscore) { tmpbestscore=tmpscore; tmpbest=slb->second; }
      }
      sit->second->bestseg=tmpbest;
      sit->second->V=tmpbestscore;
      
//      sit->second->V=(sit->second->j-sit->second->i)*as.match-gapdiff*as.gapext-(gapdiff>0?as.gapopen:0)+ MIN(sit->second->j-slb->second->i,sit->second->j2-slb->second->i2)*(9.0/10.0*as.match-1.0/10.0*as.mismatch) +slb->second->V;
/*
      sit->second->bestseg=0x00;
      double newscore=-gapdiff*as.gapext-(gapdiff>0?as.gapopen:0)+ MIN(sit->second->j-slb->second->i,sit->second->j2-slb->second->i2)*(7.0/8.0*as.match-1.0/8.0*as.mismatch) +slb->second->V;
      if (sit->second->V+MIN(sit->second->i,sit->second->j)*(7.0/8.0*as.match-1.0/8.0*as.mismatch) < newscore){
        sit->second->V=newscore;
        sit->second->bestseg=slb->second;
      }else{
        sit->second->V+=MIN(sit->second->i,sit->second->j)*(7.0/8.0*as.match-1.0/8.0*as.mismatch);
      }
*/
//        double newscore=(sit->second->j-sit->second->i)*as.match - gapdiff*as.gapext-(gapdiff>0?as.gapopen:0) + slb->second->V;
//      if (sit->second->V < newscore){
//        sit->second->V=newscore;
//        sit->second->bestseg=slb->second;
//      }
//      cout << "linking to: " << sit->second->i << " " << slb->second->j << endl;
    } else {  // end of rectangle
      slb=Lsegments.lower_bound(sit->second->j2);
      if (slb==Lsegments.end() || sit->second->V > slb->second->V){
        while (slb!=Lsegments.end()){
          if (slb->second->j2>=sit->second->j2 && slb->second->V<=sit->second->V){
            Lsegments.erase(slb++);
            continue;
          }
          ++slb;
        }
        if (bestseg==0x00 || sit->second->V > bestseg->V)
          bestseg=sit->second;
        Lsegments.insert(pair<long,ediag*>(sit->second->j2,sit->second));
      }
    }
  }

//  estr fas1,fas2;
//  cout << "full alignment" << endl;
//  cout << "ali: " << 0 <<","<<s2.seqlen<<" : " << 0 <<","<<s1.seqlen << endl; 
//  adata._score=seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,fas1,fas2);
//  print_seqali(fas1,fas2);

//  estr pas1,pas2;
//  seqident_seg_local(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,tncounts);
//  seqident_seg_sw_right(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,tncounts);
//  estr as1,as2;
//  seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,as1,as2);
//  print_seqali(as1,as2);

/*  
  if (bestseg->bestseg==0x00){ // single segment hit, do full alignment
    seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,adata,sws.alignws,as);
    return;
  }
*/

  ldieif(s2.seqlen<bestseg->j2,"segment start larger than sequence length: "+estr(s2.seqlen)+" < "+estr(bestseg->j2));
  long minright=MIN(s1.seqlen-bestseg->j,s2.seqlen-bestseg->j2);
//  seqcalign_global_norightedgegap(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,alignws);
  seqcalign_global_norightedgegap(s1,bestseg->j,MIN(s1.seqlen,bestseg->j+2*minright),s2,bestseg->j2,MIN(s2.seqlen,bestseg->j2+2*minright),adata,sws.alignws,as);
  LDEBUG(D_SEQIDENT,cout << "rightid: " << bestseg->j << "," << s1.seqlen << " l: " << s1.seqlen-bestseg->j << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
  while (bestseg!=0x00 && bestseg->bestseg != 0x00){
    seqncount(s1,bestseg->i,bestseg->j,s2,bestseg->i2,bestseg->j2,adata,as);
    LDEBUG(D_SEQIDENT,cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
    LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
    long gapdiff=abs((bestseg->bestseg->i-bestseg->bestseg->i2)-(bestseg->i-bestseg->i2));
    if (gapdiff==0) /* fast alignment no gaps */
      seqident_seg_nogaps(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata,as);
    else /* full sw alignment */
      seqcalign_global(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata,sws.alignws,as);
    LDEBUG(D_SEQIDENT,cout << "segid: " << bestseg->bestseg->j << "," << bestseg->i << " l: " << bestseg->i-bestseg->bestseg->j << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << " gapdiff: " << gapdiff << endl);
    LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
//    cout << "segid: " << bestseg->i << "," << bestseg->j << " l: " << bestseg->j-bestseg->i << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl;
    bestseg=bestseg->bestseg;
  }
  seqncount(s1,bestseg->i,bestseg->j,s2,bestseg->i2,bestseg->j2,adata,as);
  LDEBUG(D_SEQIDENT,cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
//  seqident_seg_left_local(s1,bestseg->i,s2,bestseg->i2,adata,tncounts);
  long minleft=MIN(bestseg->i,bestseg->i2);
  seqcalign_global_noleftedgegap(s1,MAX(0,bestseg->i-2*minleft),bestseg->i,s2,MAX(0,bestseg->i2-2*minleft),bestseg->i2,adata,sws.alignws,as);

//  estr as1,as2;
//  seqcalign_global_noleftedgegap(s1,0,bestseg->i,s2,0,bestseg->i2,as1,as2);
//  cout << endl;
//  cout << as1 << endl << as2 << endl;
  LDEBUG(D_SEQIDENT,cout << "leftid: "<< bestseg->i << " l: " << (bestseg->i<bestseg->i2?bestseg->i:bestseg->i2) << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,0,s1.seqlen));
//  cout << "partial alignment" << endl;
//  print_seqali(pas1,pas2);

  return;
}

void seqident_local(const eseq& s1,eintarray& kmerpos1,const eseq& s2,ealigndata& adata,esearchws& sws,const ealignscore& as,long s1_start=0,long s1_end=0)
{
  long p1,p2;
  unsigned long *pstr1=reinterpret_cast<unsigned long*>(s1.seq._str);
  unsigned long *pstr2=reinterpret_cast<unsigned long*>(s2.seq._str);
  unsigned long v1,v2;
  long lastkmerlen,lastkmerpos,lastndelta;
  if (s1_end==0) s1_end=s1.seqlen;
  long s1_len=s1_end-s1_start;

//  for (int i=0; i<11; ++i)
//    match[i]=0;
  adata.matches=0; adata.mismatches=0; adata.gaps=0; adata._score=0.0;
  adata.s1=-1; adata.e1=-1;
  adata.s2=-1; adata.e2=-1;

  char tmp[5];
  char itmp;
  int k;

  ebasicarray<ediag> diags;
  ebasicarray<ediag> dheap;

  t2.reset();
  if (s1_len<s2.seqlen){
    lastkmerpos=s1_start-2*long(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1_end-s2.seqlen;
    for (p1=s1_start; p1<s1_end-32; p1+=k){
      v1=pstr1[p1/32u]>>(2u*(p1%32u));
      v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
      for (k=0; k<32u-KMERSIZE; ++k,v1>>=2u){
        long kpos2=sws.kmerpos2[v1&KMERMASK]-sws.offset2;
        if (kpos2>-1){
          long d=p1+k-lastkmerpos;
          
          if (d<=KMERSIZE && ((lastndelta==kpos2-p1-k) || (p1+k+lastndelta+KMERSIZE<=s2.seqlen && (v1&KMERMASK)==seqkmer(s2,p1+k+lastndelta)))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p1: " << lastkmerpos-lastkmerlen+KMERSIZE<< "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=kpos2-p1-k;
          }
          lastkmerpos=p1+k;
        }
      }
    }
    v1=pstr1[p1/32u]>>(2u*(p1%32u));
    v1|=(pstr1[p1/32u+1u]<<(64u-2u*(p1%32u)))&safe_shift[p1%32u];
    for (k=0; p1+k+KMERSIZE<s1_end; ++k,v1>>=2u){
      long kpos2=sws.kmerpos2[v1&KMERMASK]-sws.offset2;
      if (kpos2>-1){
        long d=p1+k-lastkmerpos;
        if (d<=KMERSIZE && ((lastndelta==kpos2-p1-k) || (p1+k+lastndelta+KMERSIZE<=s2.seqlen && (v1&KMERMASK)==seqkmer(s2,p1+k+lastndelta)))){
          lastkmerlen+=d;
        }else{
          if (lastkmerlen){
            diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
          }
          lastkmerlen=KMERSIZE;
          lastndelta=kpos2-p1-k;
        }
        lastkmerpos=p1+k;
      }
    }
    if (lastkmerlen){
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE,lastndelta,lastkmerlen));
//      cout << "end p1: " << lastkmerpos-lastkmerlen+KMERSIZE << "," << lastkmerpos+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
    }
  }else{
    lastkmerpos=-2*long(KMERSIZE);
    lastkmerlen=0;
    lastndelta=-s1_end-s2.seqlen;
    for (p2=0; p2<s2.seqlen-32; p2+=k){
      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-KMERSIZE; ++k,v2>>=2u){
        long kpos1=kmerpos1[v2&KMERMASK]-sws.offset+s1_start;
        if (kpos1>-1){
          long d=p2+k-lastkmerpos;
          if (d<=KMERSIZE && ((lastndelta==kpos1-p2-k) || (p2+k-lastndelta+KMERSIZE<=s1_end && (v2&KMERMASK)==seqkmer(s1,p2+k-lastndelta)))){
            lastkmerlen+=d;
          }else{
            if (lastkmerlen){
//              cout << "p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
              diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
            }
            lastkmerlen=KMERSIZE;
            lastndelta=p2+k-kpos1;
          }
          lastkmerpos=p2+k;
        }
      }
    }
    v2=pstr2[p2/32u]>>(2u*(p2%32u));
    v2|=(pstr2[p2/32u+1u]<<(64u-2u*(p2%32u)))&safe_shift[p2%32u];
    for (k=0; p2+k+KMERSIZE<s2.seqlen; ++k,v2>>=2u){
      long kpos1=kmerpos1[v2&KMERMASK]-sws.offset+s1_start;
      if (kpos1>-1){
        long d=p2+k-lastkmerpos;
        if (d<=KMERSIZE && ((lastndelta==kpos1-p2-k) || (p2+k-lastndelta+KMERSIZE<=s1_end && (v2&KMERMASK)==seqkmer(s1,p2+k-lastndelta)))){
          lastkmerlen+=d;
        }else{
          if (lastkmerlen){
            diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
          }
          lastkmerlen=KMERSIZE;
          lastndelta=p2+k-kpos1;
        }
        lastkmerpos=p2+k;
      }
    }
    if (lastkmerlen){
//      cout << "end p2: " << lastkmerpos-lastkmerlen+KMERSIZE << " len: " << lastkmerlen << " ndelta: " << lastndelta << endl;
      diags.add(ediag(lastkmerpos-lastkmerlen+KMERSIZE-lastndelta,lastndelta,lastkmerlen));
    }
  }
  tdp1=tdp1*0.99+t2.lap()*0.01;

  if (diags.size()==0){
//    lderror("no shared segments found");
    return;
  }

  // sparse dynamic programming to find longest chain of segments
  ediag *bestseg=&diags[0];
  multimap<long,ediag*> segments;
  for (long j=0; j<diags.size(); ++j){
    ediag& diag(diags.at(j));
    diag.V=(diag.j-diag.i)*as.match;
    if (diag.V > bestseg->V) bestseg=&diag;
    segments.insert(pair<long,ediag*>(diag.i,&diag));
    segments.insert(pair<long,ediag*>(diag.j,&diag));
  }

  multimap<long,ediag*> Lsegments;
  multimap<long,ediag*>::iterator sit;
  multimap<long,ediag*>::iterator slb;
  for (sit=segments.begin(); sit!=segments.end(); ++sit){
//      cout << "x: " << sit->first << " start: " << (sit->first==sit->second->i?"1":"0") << " i: " << sit->second->i << " i2: " << sit->second->i2 << " Lsegs: " << Lsegments.size() << endl;
    if (sit->first==sit->second->i) { // start of rectangle
//      if (Lsegments.size()==0) continue;
//      slb=Lsegments.upper_bound(sit->second->i2-1);
//      if (slb==Lsegments.end()) continue;
      ediag *tmpbest=0x00;
      double tmpbestscore=(sit->second->j-sit->second->i)*as.match;
      for (slb=Lsegments.begin(); slb!=Lsegments.end(); ++slb){
        if (slb->first >= sit->second->i2) break;
//      sit->second->V=sit->second->j-sit->second->i+slb->second->V;
        long gapdiff=abs((sit->second->i-sit->second->i2)-(slb->second->i-slb->second->i2));
        double tmpscore=(sit->second->j-sit->second->i)*as.match - gapdiff*as.gapext-(gapdiff>0?as.gapopen:0) + slb->second->V;
        if (tmpscore>tmpbestscore) { tmpbestscore=tmpscore; tmpbest=slb->second; }
      }
      sit->second->bestseg=tmpbest;
      sit->second->V=tmpbestscore;
      
//      sit->second->V=(sit->second->j-sit->second->i)*as.match-gapdiff*as.gapext-(gapdiff>0?as.gapopen:0)+ MIN(sit->second->j-slb->second->i,sit->second->j2-slb->second->i2)*(9.0/10.0*as.match-1.0/10.0*as.mismatch) +slb->second->V;
/*
      sit->second->bestseg=0x00;
      double newscore=-gapdiff*as.gapext-(gapdiff>0?as.gapopen:0)+ MIN(sit->second->j-slb->second->i,sit->second->j2-slb->second->i2)*(7.0/8.0*as.match-1.0/8.0*as.mismatch) +slb->second->V;
      if (sit->second->V+MIN(sit->second->i,sit->second->j)*(7.0/8.0*as.match-1.0/8.0*as.mismatch) < newscore){
        sit->second->V=newscore;
        sit->second->bestseg=slb->second;
      }else{
        sit->second->V+=MIN(sit->second->i,sit->second->j)*(7.0/8.0*as.match-1.0/8.0*as.mismatch);
      }
*/
//        double newscore=(sit->second->j-sit->second->i)*as.match - gapdiff*as.gapext-(gapdiff>0?as.gapopen:0) + slb->second->V;
//      if (sit->second->V < newscore){
//        sit->second->V=newscore;
//        sit->second->bestseg=slb->second;
//      }
//      cout << "linking to: " << sit->second->i << " " << slb->second->j << endl;
    } else {  // end of rectangle
      slb=Lsegments.lower_bound(sit->second->j2);
      if (slb==Lsegments.end() || sit->second->V > slb->second->V){
        while (slb!=Lsegments.end()){
          if (slb->second->j2>=sit->second->j2 && slb->second->V<=sit->second->V){
            Lsegments.erase(slb++);
            continue;
          }
          ++slb;
        }
        if (bestseg==0x00 || sit->second->V > bestseg->V)
          bestseg=sit->second;
        Lsegments.insert(pair<long,ediag*>(sit->second->j2,sit->second));
      }
    }
  }
  tdp2=tdp2*0.99+t2.lap()*0.01;

//  estr fas1,fas2;
//  cout << "full alignment" << endl;
//  cout << "ali: " << 0 <<","<<s2.seqlen<<" : " << 0 <<","<<s1.seqlen << endl; 
//  adata._score=seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,fas1,fas2);
//  print_seqali(fas1,fas2);

//  estr pas1,pas2;
//  seqident_seg_local(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,tncounts);
//  seqident_seg_sw_right(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,tncounts);
//  estr as1,as2;
//  seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,as1,as2);
//  print_seqali(as1,as2);

/*  
  if (bestseg->bestseg==0x00){ // single segment hit, do full alignment
    seqcalign_global_noedgegap(s1,0,s1.seqlen,s2,0,s2.seqlen,adata,sws.alignws,as);
    return;
  }
*/


  float tmpflanks=0.0;
  ldieif(s2.seqlen<bestseg->j2,"segment start larger than sequence length: "+estr(s2.seqlen)+" < "+estr(bestseg->j2));
  long minright=MIN(s1.seqlen-bestseg->j,s2.seqlen-bestseg->j2);
//  seqcalign_global_norightedgegap(s1,bestseg->j,s1.seqlen,s2,bestseg->j2,s2.seqlen,adata,alignws);
  seqcalign_local_rightext(s1,bestseg->j,MIN(s1.seqlen,bestseg->j+2*minright),s2,bestseg->j2,MIN(s2.seqlen,bestseg->j2+2*minright),adata,sws.alignws,as);
  tmpflanks+=t2.lap();
  adata.profile.add(AT_RIGHT);
  LDEBUG(D_SEQIDENT,cout << "rightid: " << bestseg->j << "," << s1_end << " l: " << s1_end-bestseg->j << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,s1_start,s1_end));
  while (bestseg!=0x00 && bestseg->bestseg != 0x00){
    seqncount(s1,bestseg->i,bestseg->j,s2,bestseg->i2,bestseg->j2,adata,as);
    adata.profile.add(AT_ID);
    LDEBUG(D_SEQIDENT,cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
    LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,s1_start,s1_end));
    long gapdiff=abs((bestseg->bestseg->i-bestseg->bestseg->i2)-(bestseg->i-bestseg->i2));
//    if (gapdiff==0){ /* fast alignment no gaps: has a problem when there is a deletion and insertion in the same stretch */    
//      seqident_seg_nogaps(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata,as);
//      adata.profile.add(AT_ALIGNF);
//    }else{ /* full sw alignment */
      seqcalign_global(s1,bestseg->bestseg->j,bestseg->i,s2,bestseg->bestseg->j2,bestseg->i2,adata,sws.alignws,as);
      adata.profile.add(AT_ALIGN);
//    }
    LDEBUG(D_SEQIDENT,cout << "segid: " << bestseg->bestseg->j << "," << bestseg->i << " l: " << bestseg->i-bestseg->bestseg->j << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << " gapdiff: " << gapdiff << endl);
    LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,s1_start,s1_end));
//    cout << "segid: " << bestseg->i << "," << bestseg->j << " l: " << bestseg->j-bestseg->i << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl;
    bestseg=bestseg->bestseg;
  }
  seqncount(s1,bestseg->i,bestseg->j,s2,bestseg->i2,bestseg->j2,adata,as);
  adata.profile.add(AT_ID);
  LDEBUG(D_SEQIDENT,cout << "seg: "<< bestseg->i << "," << bestseg->j << " l: " << (bestseg->j-bestseg->i) << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,s1_start,s1_end));


  tdpmd=tdpmd*0.99+t2.lap()*0.01;
//  seqident_seg_left_local(s1,bestseg->i,s2,bestseg->i2,adata,tncounts);
  long minleft=MIN(bestseg->i,bestseg->i2);
  seqcalign_local_leftext(s1,MAX(0,bestseg->i-2*minleft),bestseg->i,s2,MAX(0,bestseg->i2-2*minleft),bestseg->i2,adata,sws.alignws,as);
  adata.profile.add(AT_LEFT);
  tmpflanks+=t2.lap();
  tdpfl=tdpfl*0.99+tmpflanks*0.01;

//  estr as1,as2;
//  seqcalign_global_noleftedgegap(s1,0,bestseg->i,s2,0,bestseg->i2,as1,as2);
//  cout << endl;
//  cout << as1 << endl << as2 << endl;
  LDEBUG(D_SEQIDENT,cout << "leftid: "<< bestseg->i << " l: " << (bestseg->i<bestseg->i2?bestseg->i:bestseg->i2) << " d: " << bestseg->i2-bestseg->i << " m: " << adata.matches << " miss: " << adata.mismatches << " gaps: " << adata.gaps << endl);
  LDEBUG(D_SEQALIGNMENT,print_tncount(adata.aln,s1_start,s1_end));
//  cout << "partial alignment" << endl;
//  print_seqali(pas1,pas2);

  return;
}



void kmercount_mask(estrarrayof<eseq>& seqs,eintarray& seqids,uint64_t *kmermask,int maskid,ebasicarray<uint32_t>& idcount,unsigned int *akmers,unsigned long akmask)
{
  int k;
  long p;
  long kmerpos;
  for (int i=0; i<seqids.size(); ++i){
    if (seqids[i]<seqs.size()){
      eseq &s(seqs.values(seqids[i]));
      unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
      unsigned long v;
      kmerpos=-long(KMERSIZE2);
      for (p=0; p+32<s.seqlen; p+=32-KMERSIZE2){
        v=pstr[p/32u]>>(2u*(p%32u));
        v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
        for (k=0; k+KMERSIZE2<32; ++k,v>>=2u){
          if ((kmermask[(v&KMERMASK2)/64u]>>((v&KMERMASK2)%64u))&0x01ul){
            unsigned long d=p+k-kmerpos;
            if (d>KMERSIZE2) d=KMERSIZE2;
            idcount[i]+=d;
//            idcount[i]+=d1_lt[(p+k-kmerpos)&BMASK16];
            kmerpos=p+k;
          }
        }
      }
      v=pstr[p/32u]>>(2u*(p%32u));
      v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
      for (k=0; p+k+KMERSIZE2<s.seqlen; ++k,v>>=2u){
        if ((kmermask[(v&KMERMASK2)/64u]>>((v&KMERMASK2)%64u))&0x01ul){
          unsigned long d=p+k-kmerpos;
          if (d>KMERSIZE2) d=KMERSIZE2;
          idcount[i]+=d;
//          idcount[i]+=d1_lt[(p+k-kmerpos)&BMASK16];
          kmerpos=p+k;
        }
      }
    }else{
      eseq &s(seqs.values(seqids[i]-seqs.size()));
      unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
      unsigned long v;
      kmerpos=-long(KMERSIZE2);
      for (p=0; p+32<s.seqlen; p+=32-KMERSIZE2){
        v=pstr[p/32u]>>(2u*(p%32u));
        v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
        for (k=0; k+KMERSIZE2<32; ++k,v>>=2u){
          uint32_t tmpv=kmer_rev_lt2[v&KMERMASK2];
          if ((kmermask[tmpv/64u]>>(tmpv%64u))&0x01ul){
            unsigned long d=p+k-kmerpos;
            if (d>KMERSIZE2) d=KMERSIZE2;
            idcount[i]+=d;
//            idcount[i]+=d1_lt[(p+k-kmerpos)&BMASK16];
            kmerpos=p+k;
          }
        }
      }
      v=pstr[p/32u]>>(2u*(p%32u));
      v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
      for (k=0; p+k+KMERSIZE2<s.seqlen; ++k,v>>=2u){
        uint32_t tmpv=kmer_rev_lt[v&KMERMASK2];
        if ((kmermask[tmpv/64u]>>(tmpv%64u))&0x01ul){
          unsigned long d=p+k-kmerpos;
          if (d>KMERSIZE2) d=KMERSIZE2;
          idcount[i]+=d;
//          idcount[i]+=d1_lt[(p+k-kmerpos)&BMASK16];
          kmerpos=p+k;
        }
      }
    }
  }
}

void kmercount_mask(estrarrayof<eseq>& seqs,eintarray& seqids,eintarray& kmermask,int maskid,ebasicarray<uint32_t>& idcount)
{
  int k;
  long p;
  long kmerpos;
  for (int i=0; i<seqids.size(); ++i){
    if (seqids[i]<seqs.size()){
      eseq &s(seqs.values(seqids[i]));
      unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
      unsigned long v;
      kmerpos=-long(KMERSIZE2);
      for (p=0; p+32<s.seqlen; p+=32-KMERSIZE2){
        v=pstr[p/32u]>>(2u*(p%32u));
        v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
        for (k=0; k+KMERSIZE2<32; ++k,v>>=2u){
          if (kmermask[v&KMERMASK2]==maskid){
//          if ((kmermask[(v&KMERMASK)/64u]>>((v&KMERMASK)%64u))&0x01ul){
            unsigned long d=p+k-kmerpos;
            if (d>KMERSIZE2) d=KMERSIZE2;
            idcount[i]+=d;
            kmerpos=p+k;
          }
        }
      }
      v=pstr[p/32u]>>(2u*(p%32u));
      v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
      for (k=0; p+k+KMERSIZE2<s.seqlen; ++k,v>>=2u){
        if (kmermask[v&KMERMASK2]==maskid){
//        if ((kmermask[(v&KMERMASK)/64u]>>((v&KMERMASK)%64u))&0x01ul){
          unsigned long d=p+k-kmerpos;
          if (d>KMERSIZE2) d=KMERSIZE2;
          idcount[i]+=d;
          kmerpos=p+k;
        }
      }
    }else{
      eseq &s(seqs.values(seqids[i]-seqs.size()));
      unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
      unsigned long v;
      kmerpos=-long(KMERSIZE2);
      for (p=0; p+32<s.seqlen; p+=32-KMERSIZE2){
        v=pstr[p/32u]>>(2u*(p%32u));
        v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
        for (k=0; k+KMERSIZE2<32; ++k,v>>=2u){
//          uint32_t tmpv=kmer_rev_lt[v&KMERMASK];
//          if ((kmermask[tmpv/64u]>>(tmpv%64u))&0x01ul){
          if (kmermask[kmer_rev_lt2[v&KMERMASK2]]==maskid){
            unsigned long d=p+k-kmerpos;
            if (d>KMERSIZE2) d=KMERSIZE2;
            idcount[i]+=d;
            kmerpos=p+k;
          }
        }
      }
      v=pstr[p/32u]>>(2u*(p%32u));
      v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
      for (k=0; p+k+KMERSIZE2<s.seqlen; ++k,v>>=2u){
//        uint32_t tmpv=kmer_rev_lt[v&KMERMASK];
//        if ((kmermask[tmpv/64u]>>(tmpv%64u))&0x01ul){
        if (kmermask[kmer_rev_lt2[v&KMERMASK2]]==maskid){
          unsigned long d=p+k-kmerpos;
          if (d>KMERSIZE2) d=KMERSIZE2;
          idcount[i]+=d;
          kmerpos=p+k;
        }
      }
    }
  }
}

void setkmermask(uint64_t *kmermask,eseq& s,unsigned int *akmers,unsigned long akmask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
//  unsigned int akmask=0xFul;
  int k;
  long p;
  for (p=0; p<s.seqlen-32; p+=32-KMERSIZE2){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE2; ++k,v>>=2u){
      if (akmers[v&akmask]==0x0u) continue;
      uint32_t tmpv=v&KMERMASK2;
      kmermask[tmpv/64u]|=0x1u<<(tmpv%64u);
    }
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k<s.seqlen-long(KMERSIZE2); ++k,v>>=2u){
    if (akmers[v&akmask]==0x0u) continue;
    uint32_t tmpv=v&KMERMASK2;
    kmermask[tmpv/64u]|=0x1u<<(tmpv%64u);
  }
}

void setkmermask(eintarray& kmermask,eseq& s,int maskid)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
//  unsigned int akmask=0xFul;
  int k;
  long p;
  for (p=0; p<s.seqlen-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u)
      kmermask[v&KMERMASK]=maskid;
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k<s.seqlen-long(KMERSIZE); ++k,v>>=2u)
    kmermask[v&KMERMASK]=maskid;
}

void setkmermask(eintarray& kmermask,eseq& s,int maskid,unsigned int *akmers,unsigned long akmask,unsigned long start,unsigned long end)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
//  unsigned int akmask=0xFul;
  int k;
  long p;
  for (p=start; p+32<end; p+=32-KMERSIZE2){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k+KMERSIZE2<32; ++k,v>>=2u){
      if (akmers[v&akmask]==0x0u) continue;
      kmermask[v&KMERMASK2]=maskid;
    }
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE2<end; ++k,v>>=2u){
    if (akmers[v&akmask]==0x0u) continue;
    kmermask[v&KMERMASK2]=maskid;
  }
}

void setkmermask(eintarray& kmermask,eseq& s,int maskid,unsigned int *akmers,unsigned long akmask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
//  unsigned int akmask=0xFul;
  int k;
  long p;
  for (p=0; p+32<s.seqlen; p+=32-KMERSIZE2){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE2; ++k,v>>=2u){
      if (akmers[v&akmask]==0x0u) continue;
      kmermask[v&KMERMASK2]=maskid;
    }
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k<s.seqlen-int(KMERSIZE2); ++k,v>>=2u){
    if (akmers[v&akmask]==0x0u) continue;
    kmermask[v&KMERMASK2]=maskid;
  }
}

void setkmerpos(eintarray& kmerpos,eseq& s,int offset,unsigned long start,unsigned long end)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p;
  for (p=start; p+32<end; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k+KMERSIZE<32; ++k,v>>=2u)
      kmerpos[v&KMERMASK]=offset+(p+k-start);
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE<end; ++k,v>>=2u)
    kmerpos[v&KMERMASK]=offset+(p+k-start);
}

void setkmerpos(eintarray& kmerpos,eseq& s,int offset)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p;
  for (p=0; p+32<s.seqlen; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k+KMERSIZE<32; ++k,v>>=2u)
      kmerpos[v&KMERMASK]=offset+p+k;
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE<s.seqlen; ++k,v>>=2u)
    kmerpos[v&KMERMASK]=offset+p+k;
}



void kmercount_single3(ebasicarray<deque<int>*>& otukmers,eseq& s,uint32_t sp,uint32_t rp,ebasicarray<uint32_t>& idcount,eintarray& kmerpos)
//void kmercount_single2(ebasicarray<deque<int>*>& otukmers,eseq& s,uint32_t rp,eintarray& idcount,eintarray& kmerpos)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  int p=0;
  sp<<=10u;
  for (; p<long(s.seqlen)-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if (otukmers[v&KMERMASK]){
        deque<int> &otuind(*otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l){
          unsigned int tl=otuind[l];
          unsigned int d=rp+p+k-kmerpos[tl];
          unsigned int &idc(idcount[tl]);
//          d=d1_lt[d&BMASK16]|d2_lt[(d>>16u)&BMASK16];
          if (d>KMERSIZE) d=KMERSIZE;
          uint32_t tmp=(id_lt[((idc^sp)>>10u)&BMASK11]*id_lt[((idc^sp)>>21u)&BMASK11]*idc&BMASK10);
          tmp=(((tmp+d)&BMASK10)<tmp?1024u:(tmp+d)&BMASK10);
//          if (d==1 || d==KMERSIZE)
          idcount[tl]=sp|tmp;
//          idcount[tl]+=d;
          kmerpos[tl]=rp+p+k;
        }
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k<long(s.seqlen)-long(KMERSIZE); ++k,v>>=2u){
    if (otukmers[v&KMERMASK]){
      deque<int> &otuind(*otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        unsigned int tl=otuind[l];
        unsigned int d=rp+p+k-kmerpos[tl];
        unsigned int &idc(idcount[tl]);
        if (d>KMERSIZE) d=KMERSIZE;
//        d=d1_lt[d&BMASK16]|d2_lt[(d>>16u)&BMASK16];
//        idcount[tl]+=d;
        idcount[tl]=sp|((id_lt[((idc^sp)>>10u)&BMASK11]*id_lt[((idc^sp)>>21u)&BMASK11]*idc)&BMASK10)+d;
        kmerpos[tl]=rp+p+k;
      }
    }
  }
}

void kmercount_single2(ebasicarray<deque<int>*>& otukmers,eseq& s,uint32_t rp,ebasicarray<uint32_t>& idcount,eintarray& kmerpos)
//void kmercount_single2(ebasicarray<deque<int>*>& otukmers,eseq& s,uint32_t rp,eintarray& idcount,eintarray& kmerpos)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  int p=0;
  for (; p<long(s.seqlen)-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if (otukmers[v&KMERMASK]){
        deque<int> &otuind(*otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l){
          unsigned int tl=otuind[l];
          unsigned int d=rp+p+k-kmerpos[tl];
//          d=d1_lt[d&BMASK16]|d2_lt[(d>>16u)&BMASK16];
          if (d>KMERSIZE) d=KMERSIZE;
//          if (d==1 || d==KMERSIZE)
//          idcount[tl]=sp|((id_lt[((idc^sp)>>10u)&BMASK11]*id_lt[((idc^sp)>>21u)&BMASK11]*idc)&BMASK10)+d;
          idcount[tl]+=d;
          kmerpos[tl]=rp+p+k;
        }
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k<long(s.seqlen)-long(KMERSIZE); ++k,v>>=2u){
    if (otukmers[v&KMERMASK]){
      deque<int> &otuind(*otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        unsigned int tl=otuind[l];
        unsigned int d=rp+p+k-kmerpos[tl];
        unsigned int &idc(idcount[tl]);
        if (d>KMERSIZE) d=KMERSIZE;
//        d=d1_lt[d&BMASK16]|d2_lt[(d>>16u)&BMASK16];
        idcount[tl]+=d;
        kmerpos[tl]=rp+p+k;
      }
    }
  }
}

void kmercount_single(ebasicarray<deque<uint32_t> >& otukmers,eseq& s,ebasicarray<uint32_t>& idcount,eintarray& kmerpos)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p=0;
  for (; p<long(s.seqlen)-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
//      if (otukmers[v&KMERMASK]){
        deque<uint32_t> &otuind(otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l){
          unsigned int tl=otuind[l];
          unsigned long d=p+k-kmerpos[tl];
          if (d>KMERSIZE) d=KMERSIZE;
//          if (d==1 || d==KMERSIZE)
            idcount[tl]+=d;
          kmerpos[tl]=p+k;
        }
//      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k<long(s.seqlen)-long(KMERSIZE); ++k,v>>=2u){
//    if (otukmers[v&KMERMASK]){
      deque<uint32_t> &otuind(otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        unsigned int tl=otuind[l];
        unsigned long d=p+k-kmerpos[tl];
        if (d>KMERSIZE) d=KMERSIZE;
//        if (d==1 || d==KMERSIZE)
          idcount[tl]+=d;
        kmerpos[tl]=p+k;
      }
//    }
  }
}

uint32_t bmask_lt[2]={0x0u,0xffffffffu};
int16_t idc_lt[2]={1,-1};

void kmercount_both_nopos2_skip_plist(int scount,ebasicarray<deque<uint32_t> >& otukmers,eseq& s,unsigned long start,unsigned long end,ebasicarray<uint32_t>& idcount,uint32_t* &plist,unsigned int *akmers,unsigned long akmask,int skipn)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p=start;
  for (; p+32<end; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if ((p+k)%skipn) continue;
//      cout << kmer2str(v&KMERMASK) << " " << kmer2str(kmer_rev_lt[v&KMERMASK]) << endl;
//      if (otukmers[v&KMERMASK]){
      if (akmers[v&akmask]){
        deque<uint32_t> &otuind(otukmers[v&KMERMASK]);
        unsigned int oisize=otuind.size();
        for (unsigned int l=0; l<oisize; ++l){
          *plist=otuind[l];
          uint32_t &idc(idcount[*plist]);
          plist+=(~idc)&0x1u;
          idc|=0x1u;
          idc+=2u;
        }
      }
      if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
//      if (otukmers[kmer_rev_lt[v&KMERMASK]]){
        deque<uint32_t> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
        unsigned int oisize=otuind.size();
        for (unsigned int l=0; l<oisize; ++l){
          *plist=otuind[l]+scount;
          uint32_t &idc(idcount[*plist]);
          plist+=(~idc)&0x1u;
          idc|=0x1u;
          idc+=2u;
        }
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE<end; ++k,v>>=2u){
    if ((p+k)%skipn) continue;
//    if (otukmers[v&KMERMASK]){
    if (akmers[v&akmask]){
      deque<uint32_t> &otuind(otukmers[v&KMERMASK]);
      unsigned int oisize=otuind.size();
      for (unsigned int l=0; l<oisize; ++l){
        *plist=otuind[l];
        uint32_t &idc(idcount[*plist]);
        plist+=(~idc)&0x1u;
        idc|=0x1u;
        idc+=2u;
      }
    }
//    if (otukmers[kmer_rev_lt[v&KMERMASK]]){
    if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
      deque<uint32_t> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
      unsigned int oisize=otuind.size();
      for (unsigned int l=0; l<oisize; ++l){
        *plist=otuind[l]+scount;
        uint32_t &idc(idcount[*plist]);
        plist+=(~idc)&0x1u;
        idc|=0x1u;
        idc+=2u;
      }
    }
  }
}

void kmercount_both_nopos2_skip(int scount,ebasicarray<deque<uint32_t> >& otukmers,eseq& s,unsigned long start,unsigned long end,ebasicarray<uint32_t>& idcount,unsigned int *akmers,unsigned long akmask,int skipn)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p=start;
  for (; p+32<end; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if ((p+k)%skipn) continue;
//      cout << kmer2str(v&KMERMASK) << " " << kmer2str(kmer_rev_lt[v&KMERMASK]) << endl;
//      if (otukmers[v&KMERMASK]){
      if (akmers[v&akmask]){
        deque<uint32_t> &otuind(otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l)
          ++idcount[otuind[l]];
      }
      if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
//      if (otukmers[kmer_rev_lt[v&KMERMASK]]){
        deque<uint32_t> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
        for (unsigned int l=0; l<otuind.size(); ++l)
          ++idcount[scount+otuind[l]];
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE<end; ++k,v>>=2u){
    if ((p+k)%skipn) continue;
//    if (otukmers[v&KMERMASK]){
    if (akmers[v&akmask]){
      deque<uint32_t> &otuind(otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l)
        ++idcount[otuind[l]];
    }
//    if (otukmers[kmer_rev_lt[v&KMERMASK]]){
    if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
      deque<uint32_t> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
      for (unsigned int l=0; l<otuind.size(); ++l)
        ++idcount[scount+otuind[l]];
    }
  }
}

void kmercount_both_nopos2(int scount,ebasicarray<deque<uint32_t> >& otukmers,eseq& s,eintarray& idcount,eintarray& kmerpos,unsigned int *akmers,unsigned long akmask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p=0;
  for (; p+32<s.seqlen; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
//      cout << kmer2str(v&KMERMASK) << " " << kmer2str(kmer_rev_lt[v&KMERMASK]) << endl;
//      if (otukmers[v&KMERMASK]){
      if (akmers[v&akmask]){
        deque<uint32_t> &otuind(otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l)
          idcount[(otuind[l]&BMASK31)]+=idc_lt[otuind[l]>>31u];
      }
      if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
//      if (otukmers[kmer_rev_lt[v&KMERMASK]]){
        deque<uint32_t> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
        for (unsigned int l=0; l<otuind.size(); ++l)
          idcount[scount+(otuind[l]&BMASK31)]+=idc_lt[otuind[l]>>31u];
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE<s.seqlen; ++k,v>>=2u){
//    if (otukmers[v&KMERMASK]){
    if (akmers[v&akmask]){
      deque<uint32_t> &otuind(otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l)
        idcount[(otuind[l]&BMASK31)]+=idc_lt[otuind[l]>>31u];
    }
//    if (otukmers[kmer_rev_lt[v&KMERMASK]]){
    if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
      deque<uint32_t> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
      for (unsigned int l=0; l<otuind.size(); ++l)
        idcount[scount+(otuind[l]&BMASK31)]+=idc_lt[otuind[l]>>31u];
    }
  }
}

void kmercount_both_nopos(int scount,ebasicarray<deque<int> >& otukmers,eseq& s,eintarray& idcount,eintarray& kmerpos,unsigned int *akmers,unsigned long akmask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p=0;
  for (; p+32<s.seqlen; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
//      cout << kmer2str(v&KMERMASK) << " " << kmer2str(kmer_rev_lt[v&KMERMASK]) << endl;
//      if (otukmers[v&KMERMASK]){
      if (akmers[v&akmask]){
        deque<int> &otuind(otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l)
          ++idcount[otuind[l]];
      }
      if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
//      if (otukmers[kmer_rev_lt[v&KMERMASK]]){
        deque<int> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
        for (unsigned int l=0; l<otuind.size(); ++l)
          ++idcount[scount+otuind[l]];
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE<s.seqlen; ++k,v>>=2u){
//    if (otukmers[v&KMERMASK]){
    if (akmers[v&akmask]){
      deque<int> &otuind(otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l)
        ++idcount[otuind[l]];
    }
//    if (otukmers[kmer_rev_lt[v&KMERMASK]]){
    if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
      deque<int> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
      for (unsigned int l=0; l<otuind.size(); ++l)
        ++idcount[scount+otuind[l]];
    }
  }
}

void kmercount_both2(int scount,ebasicarray<deque<int> >& otukmers,eseq& s,eintarray& idcount,eintarray& kmerpos,unsigned int *akmers,unsigned long akmask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p=0;
  for (; p+32<s.seqlen; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
//      cout << kmer2str(v&KMERMASK) << " " << kmer2str(kmer_rev_lt[v&KMERMASK]) << endl;
//      if (otukmers[v&KMERMASK]){
      if (akmers[v&akmask]){
        deque<int> &otuind(otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l){
          unsigned int tl=otuind[l];
          unsigned int d=p+k+KMERSIZE-kmerpos[tl];
          if (d>KMERSIZE) d=KMERSIZE;
//          idcount[tl]+=d1_lt[(p+k+KMERSIZE-kmerpos[tl])&BMASK16];
          idcount[tl]+=d;
          kmerpos[tl]=p+k+KMERSIZE;
        }
      }
      if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
//      if (otukmers[kmer_rev_lt[v&KMERMASK]]){
        deque<int> &otuind2(otukmers[kmer_rev_lt[v&KMERMASK]]);
        for (unsigned int l=0; l<otuind2.size(); ++l){
          uint32_t tl=scount+otuind2[l];
          unsigned int d=p+k+KMERSIZE-kmerpos[tl];
          if (d>KMERSIZE) d=KMERSIZE;
          idcount[tl]+=d;
          kmerpos[tl]=p+k+KMERSIZE;
        }
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE<s.seqlen; ++k,v>>=2u){
//    if (otukmers[v&KMERMASK]){
    if (akmers[v&akmask]){
      deque<int> &otuind(otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        uint32_t tl=otuind[l];
        unsigned int d=p+k+KMERSIZE-kmerpos[tl];
        if (d>KMERSIZE) d=KMERSIZE;
        idcount[tl]+=d;
        kmerpos[tl]=p+k+KMERSIZE;
      }
    }
//    if (otukmers[kmer_rev_lt[v&KMERMASK]]){
    if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
      deque<int> &otuind2(otukmers[kmer_rev_lt[v&KMERMASK]]);
      for (unsigned int l=0; l<otuind2.size(); ++l){
        uint32_t tl=scount+otuind2[l];
        unsigned int d=p+k+KMERSIZE-kmerpos[tl];
        if (d>KMERSIZE) d=KMERSIZE;
        idcount[tl]+=d;
        kmerpos[tl]=p+k+KMERSIZE;
      }
    }
  }
}

void kmercount_both(int scount,ebasicarray<deque<int> >& otukmers,eseq& s,ebasicarray<uint32_t>& idcount,uint64_t *bitmask,unsigned int *akmers,unsigned long akmask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p=0;
  for (; p+32<s.seqlen; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
//      cout << kmer2str(v&KMERMASK) << " " << kmer2str(kmer_rev_lt[v&KMERMASK]) << endl;
      if (akmers[v&akmask]){
        deque<int> &otuind(otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l){
          unsigned int tl=otuind[l];
          uint16_t bt=(bitmask[tl/64u]>>(tl%64u))&0x1ul;
          bitmask[tl/64u]|=(0x1ul<<(tl%64u));
//          unsigned int d=p+k+KMERSIZE-(kmerpos[tl]&bmask_lt[bt]);
//          uint32_t idc=idcount[tl]&bmask_lt[bt];
          uint32_t idc=idcount[tl]&bmask_lt[bt];
          uint32_t d=d1_lt[p+k+KMERSIZE-(idc>>16u)];
//          if (d>KMERSIZE) d=KMERSIZE;
          idcount[tl]=((idc+d)&BMASK16)|((p+k+KMERSIZE)<<16u);
//          if (d>KMERSIZE) d=KMERSIZE;
//          idcount[tl]=((idc+d)&BMASK16)|((p+k+KMERSIZE)<<16u);
//          kmerpos[tl]=p+k+KMERSIZE;
//          idcount[tl]+=d;
//          kmerpos[tl]=p+k;
        }
      }
      if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
        deque<int> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
        for (unsigned int l=0; l<otuind.size(); ++l){
          uint32_t tl=scount+otuind[l];
          uint16_t bt=(bitmask[tl/64u]>>(tl%64u))&0x1ul;
          bitmask[tl/64u]|=(0x1ul<<(tl%64u));
          uint32_t idc=idcount[tl]&bmask_lt[bt];
          uint32_t d=d1_lt[p+k+KMERSIZE-(idc>>16u)];
//          uint32_t d=p+k+KMERSIZE-(idc>>16u);
//          if (d>KMERSIZE) d=KMERSIZE;
          idcount[tl]=((idc+d)&BMASK16)|((p+k+KMERSIZE)<<16u);
//          uint32_t idc=idcount[tl]&bmask_lt[bt];
//          uint32_t d=p+k+KMERSIZE-(idc>>16u);
//          if (d>KMERSIZE) d=KMERSIZE;
//          idcount[tl]=((idc+d)&BMASK16)|((p+k+KMERSIZE)<<16u);
        }
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k+KMERSIZE<s.seqlen; ++k,v>>=2u){
    if (akmers[v&akmask]){
      deque<int> &otuind(otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        uint32_t tl=otuind[l];
        uint16_t bt=(bitmask[tl/64u]>>(tl%64u))&0x1ul;
        bitmask[tl/64u]|=(0x1ul<<(tl%64u));
        uint32_t idc=idcount[tl]&bmask_lt[bt];
        uint32_t d=d1_lt[p+k+KMERSIZE-(idc>>16u)];
//        uint32_t d=p+k+KMERSIZE-(idc>>16u);
//        if (d>KMERSIZE) d=KMERSIZE;
        idcount[tl]=((idc+d)&BMASK16)|((p+k+KMERSIZE)<<16u);
//        uint32_t idc=idcount[tl]&bmask_lt[bt];
//        uint32_t d=p+k+KMERSIZE-(idc>>16u);
//        idcount[tl]=((idc+d)&BMASK16)|((p+k+KMERSIZE)<<16u);
      }
    }
    if (akmers[kmer_rev_lt[v&KMERMASK]&akmask]){
      deque<int> &otuind(otukmers[kmer_rev_lt[v&KMERMASK]]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        uint32_t tl=scount+otuind[l];
        uint16_t bt=(bitmask[tl/64u]>>(tl%64u))&0x1ul;
        bitmask[tl/64u]|=(0x1ul<<(tl%64u));
        uint32_t idc=idcount[tl]&bmask_lt[bt];
        uint32_t d=d1_lt[p+k+KMERSIZE-(idc>>16u)];
//        uint32_t d=p+k+KMERSIZE-(idc>>16u);
//        if (d>KMERSIZE) d=KMERSIZE;
        idcount[tl]=((idc+d)&BMASK16)|((p+k+KMERSIZE)<<16u);
//        uint32_t d=p+k+KMERSIZE-(idc>>16u);
      }
    }
  }
}

void kmercount(earray<eintarray>& otus,ebasicarray<eintarray*>& otukmers,int seqcount,eseq& s,eintarray& idcount,eintarray& kmerpos,eintarray& tmpkmers,eintarray& seqkmers,eintarray& revtmpkmers,eintarray& revseqkmers)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  int p=0;
  for (; p<int(s.seqlen)-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if (otukmers[v&KMERMASK]){
        eintarray &otuind(*otukmers[v&KMERMASK]);
        for (unsigned int l=0; l<otuind.size(); ++l){
          unsigned int tl=otuind[l];
          unsigned int d=p+k-kmerpos[tl];
          if (d>KMERSIZE) d=KMERSIZE;
          idcount[tl]+=d;
          kmerpos[tl]=p+k;
        }
      }
      if (tmpkmers[v&KMERMASK]!=seqcount){
        tmpkmers[v&KMERMASK]=seqcount;
        seqkmers[v&KMERMASK]=p+k;
      }
      if (otukmers[kmer_rev_lt[v&KMERMASK]]){
        eintarray &otuind(*otukmers[kmer_rev_lt[v&KMERMASK]]);
        for (unsigned int l=0; l<otuind.size(); ++l){
          unsigned int tl=otus.size()+otuind[l];
          unsigned int d=p+k-kmerpos[tl];
          if (d>KMERSIZE) d=KMERSIZE;
          idcount[tl]+=d;
          kmerpos[tl]=p+k;
        }
      }
      if (revtmpkmers[kmer_rev_lt[v&KMERMASK]]!=seqcount){
        revtmpkmers[kmer_rev_lt[v&KMERMASK]]=seqcount;
        revseqkmers[kmer_rev_lt[v&KMERMASK]]=s.seqlen-p-k-1-KMERSIZE;
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k<int(s.seqlen)-int(KMERSIZE); ++k,v>>=2u){
    if (otukmers[v&KMERMASK]){
      eintarray &otuind(*otukmers[v&KMERMASK]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        unsigned int tl=otuind[l];
        unsigned int d=p+k-kmerpos[tl];
        if (d>KMERSIZE) d=KMERSIZE;
        idcount[tl]+=d;
        kmerpos[tl]=p+k;
      }
    }
    if (tmpkmers[v&KMERMASK]!=seqcount){
      tmpkmers[v&KMERMASK]=seqcount;
      seqkmers[v&KMERMASK]=p+k;
    }
    if (otukmers[kmer_rev_lt[v&KMERMASK]]){
      eintarray &otuind(*otukmers[kmer_rev_lt[v&KMERMASK]]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        unsigned int tl=otus.size()+otuind[l];
        unsigned int d=p+k-kmerpos[tl];
        if (d>KMERSIZE) d=KMERSIZE;
        idcount[tl]+=d;
        kmerpos[tl]=p+k;
      }
    }
    if (revtmpkmers[kmer_rev_lt[v&KMERMASK]]!=seqcount){
      revtmpkmers[kmer_rev_lt[v&KMERMASK]]=seqcount;
      revseqkmers[kmer_rev_lt[v&KMERMASK]]=s.seqlen-p-k-1-KMERSIZE;
    }
  }
}


void kmercount_single_nopos(ebasicarray<deque<uint32_t> >& otukmers,eseq& s,uint64_t *bitmask,ebasicarray<uint32_t>& idcount,eintarray& tmpkmers,int ti,int& kmercount,uint32_t& bid,short& bcount,unsigned int *akmers,unsigned long akmask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  int k;
  long p=0;
  kmercount=0;
  bcount=0; bid=0;
  for (; p<long(s.seqlen)-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (k=0; k<32-KMERSIZE; ++k,v>>=2u){
      unsigned long kmer=v&KMERMASK;
      if (akmers[v&akmask]==0x0u || tmpkmers[kmer]==ti) continue;
      ++kmercount;
      tmpkmers[kmer]=ti;
      deque<uint32_t> &otuind(otukmers[kmer]);
      for (unsigned int l=0; l<otuind.size(); ++l){
        uint32_t oi=otuind[l]&BMASK31;
        uint16_t bt=(bitmask[oi/64u]>>(oi%64u))&0x1ul;
        bitmask[oi/64u]|=(0x1ul<<(oi%64u));
        idcount[oi]=(idcount[oi]&bmask_lt[bt])+idc_lt[otuind[l]>>31u];
        if (bcount<idcount[oi]){
          bcount=idcount[oi];
          bid=oi;
        }
      }
    }
  }

//  cout << "p: " << p << " seqlen: " << s.seqlen << endl;

  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (k=0; p+k<int(s.seqlen)-int(KMERSIZE); ++k,v>>=2u){
    unsigned long kmer=v&KMERMASK;
    if (akmers[v&akmask]==0x0u || tmpkmers[kmer]==ti) continue;
    ++kmercount;
    tmpkmers[kmer]=ti;
    deque<uint32_t> &otuind(otukmers[kmer]);
    for (unsigned int l=0; l<otuind.size(); ++l){
      uint32_t oi=otuind[l]&BMASK31;
      uint16_t bt=(bitmask[oi/64u]>>(oi%64u))&0x1ul;
      bitmask[oi/64u]|=(0x1ul<<(oi%64u));
      idcount[oi]=(idcount[oi]&bmask_lt[bt])+idc_lt[otuind[l]>>31u];
      if (bcount<idcount[oi]){
        bcount=idcount[oi];
        bid=oi;
      }
    }
  }
}



void otuaddkmerdiff(ebasicarray<deque<uint32_t> >& otukmers,eseq& s1,eseq& s2,eintarray& tmpkmers,int ti,int pi,unsigned int *akmers,unsigned long akmask)
{
  unsigned long v;
  long p;
  int k;
  unsigned long *pstr;
  
  pstr=reinterpret_cast<unsigned long*>(s1.seq._str);
  for (p=0; p<s1.seqlen-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (int k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if (akmers[v&akmask]==0x0u || tmpkmers[v&KMERMASK]==-ti) continue; // only add allowed kmers 
      if (tmpkmers[v&KMERMASK]!=ti) // kmer not in new sequence
        otukmers[v&KMERMASK].push_back(pi|(1u<<31u));
      tmpkmers[v&KMERMASK]=-ti;
    }
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (int k=0; p+k<long(s1.seqlen)-long(KMERSIZE); ++k,v>>=2u){
    if (akmers[v&akmask]==0x0u || tmpkmers[v&KMERMASK]==-ti) continue; // only add allowed kmers 
    if (tmpkmers[v&KMERMASK]!=ti) // kmer not in new sequence
      otukmers[v&KMERMASK].push_back(pi|(1u<<31u));
    tmpkmers[v&KMERMASK]=-ti;
  }

  // add kmer increments in new leaf
  pstr=reinterpret_cast<unsigned long*>(s2.seq._str);
  for (p=0; p<s2.seqlen-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (int k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if (akmers[v&akmask]==0x0u || tmpkmers[v&KMERMASK]==-ti) continue; // only add allowed kmers 
      tmpkmers[v&KMERMASK]=-ti;
      otukmers[v&KMERMASK].push_back(pi);
    }
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (int k=0; p+k<long(s2.seqlen)-long(KMERSIZE); ++k,v>>=2u){
    if (akmers[v&akmask]==0x0u || tmpkmers[v&KMERMASK]==-ti) continue; // only add allowed kmers 
    tmpkmers[v&KMERMASK]=-ti;
    otukmers[v&KMERMASK].push_back(pi);
  }
}




void otukmeradd(ebasicarray<deque<uint32_t> >& otukmers,int i,eseq& s,eintarray& tmpkmers,int ti,unsigned int *akmers,unsigned long akmask,ebasicarray<char>& otukmermask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  long p;
  int k;
  for (p=0; p<s.seqlen-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (int k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if (akmers[v&akmask]==0x0u || otukmermask[v&KMERMASK]==1) continue; // only add allowed kmers 
      if (tmpkmers[v&KMERMASK]==ti) continue;
      tmpkmers[v&KMERMASK]=ti;
//      if (otukmers[v&KMERMASK]==0x00){
//        otukmers[v&KMERMASK]=new deque<int>();
//        otukmers[v&KMERMASK]->reserve(1000);
//      }
      otukmers[v&KMERMASK].push_back(i);
    }
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (int k=0; p+k<long(s.seqlen)-long(KMERSIZE); ++k,v>>=2u){
    if (akmers[v&akmask]==0x0u || otukmermask[v&KMERMASK]==1) continue; // only add allowed kmers 
    if (tmpkmers[v&KMERMASK]==ti) continue;
    tmpkmers[v&KMERMASK]=ti;
//    if (otukmers[v&KMERMASK]==0x00){
//      otukmers[v&KMERMASK]=new deque<int>();
//      otukmers[v&KMERMASK]->reserve(1000);
//    }
    otukmers[v&KMERMASK].push_back(i);
  }
}

void otukmeradd(ebasicarray<deque<int> >& otukmers,int i,eseq& s,eintarray& tmpkmers,int ti,unsigned int *akmers,unsigned long akmask)
{
  unsigned long *pstr=reinterpret_cast<unsigned long*>(s.seq._str);
  unsigned long v;
  long p;
  int k;
  for (p=0; p<s.seqlen-32; p+=32-KMERSIZE){
    v=pstr[p/32u]>>(2u*(p%32u));
    v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
    for (int k=0; k<32-KMERSIZE; ++k,v>>=2u){
      if (akmers[v&akmask]==0x0u) continue; // only add allowed kmers 
      if (tmpkmers[v&KMERMASK]==ti) continue;
      tmpkmers[v&KMERMASK]=ti;
//      if (otukmers[v&KMERMASK]==0x00){
//        otukmers[v&KMERMASK]=new deque<int>();
//        otukmers[v&KMERMASK]->reserve(1000);
//      }
      otukmers[v&KMERMASK].push_back(i);
    }
  }
  v=pstr[p/32u]>>(2u*(p%32u));
  v|=(pstr[p/32u+1u]<<(64u-2u*(p%32u)))&safe_shift[p%32u];
  for (int k=0; p+k<int(s.seqlen)-int(KMERSIZE); ++k,v>>=2u){
    if (akmers[v&akmask]==0x0u) continue; // only add allowed kmers 
    if (tmpkmers[v&KMERMASK]==ti) continue;
    tmpkmers[v&KMERMASK]=ti;
//    if (otukmers[v&KMERMASK]==0x00){
//      otukmers[v&KMERMASK]=new deque<int>();
//      otukmers[v&KMERMASK]->reserve(1000);
//    }
    otukmers[v&KMERMASK].push_back(i);
  }
}

void randomize(eintarray& seqids)
{
  for (int i=seqids.size()-1; i>0; --i){
    int r=int(rnd.uniform()*(i+1));
    if (r!=i) seqids.swap(r,i);
  }
}

/*
void seqsearch_benchmark(const estr& str2id,estrarrayof<eseq>& seqs,earray<eintarray>& otus,ebasicarray<eintarray*>& otukmers,eseq& s,eintarray& otukmerpos,eintarray& idcount,int& maskid,eintarray& kmermask,int& offset,eintarray& kmerpos,eintarray& kmerposrev,int& offset2,eintarray& kmerpos2,ebasicarray<ealigndata>& matchcounts,eintarray& truetax)
{
  eseq srev(s);
  srev.revcompl();

  matchcounts.clear();
  otukmerpos.init(otus.size()*2,-int(KMERSIZE));
  idcount.init(otus.size()*2,0);

  kmercount_both(otus.size(),otukmers,s,idcount,otukmerpos);

//  eintarray best;
//  ebasicarray<ealigndata> matchcounts;

  eintarray best;
//  best.clear();
  best.add(0);
  for (int l=1; l<idcount.size(); ++l){
//    if (idcount[l]<idcount[best[best.size()-1]] && (idcount[l]<0.8*idcount[ibest] && best.size()>=20 || best.size()==100)) continue; // worst id than the bottom of the list
    if (idcount[l]<idcount[best[best.size()-1]] && best.size()==10) continue; // worst id than the bottom of the list

    int j;
    for (j=0; j<best.size() && idcount[l]<idcount[best[j]]; ++j);

    if (best.size()<10)
      best.add(l);
    for (int tj=best.size()-1; tj>j; --tj)
      best[tj]=best[tj-1];
    best[j]=l;
  }

  eintarray seqids;
  for (int l=0; l<best.size(); ++l){
    int ibest=best[l];
    if (ibest<otus.size()){
      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest][0]) << "\t" << idcount[ibest] << endl;
      for (int l2=0; l2<otus[ibest].size(); ++l2)
        seqids.add(otus[ibest][l2]);
    }else{
      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest-otus.size()][0]) << "\t" << idcount[ibest] << " reversed" << endl;
      ibest-=otus.size();
      for (int l2=0; l2<otus[ibest].size(); ++l2)
        seqids.add(seqs.size()+otus[ibest][l2]);
    }
  }
//  randomize(seqids);
  idcount.init(seqids.size(),0);

  if (maskid+1<maskid){
    maskid=-1;
    kmermask.init(KMERSIZE,-1);
  }
  ++maskid;
  setkmermask(kmermask,s,maskid,akmers,0xFul);
  kmercount_mask(seqs,seqids,kmermask,maskid,idcount);

  best.clear();
  for (int l=0; l<idcount.size(); ++l){
    eseq& tmps(seqs.values(seqids[l]<seqs.size()?seqids[l]:seqids[l]-seqs.size()));
    ldieif(tmps.tax!=0x00 && tmps.tax->size()==0,"error in tax");
    if ((best.size()>0 && idcount[l]<idcount[best[best.size()-1]] && best.size()==20) || (tmps.tax!=0x00 && (*tmps.tax)[tmps.tax->size()-1]==truetax[truetax.size()-1])) continue; // worst id than the bottom of the list, ignore same species

    int j;
    for (j=0; j<best.size() && idcount[l]<idcount[best[j]]; ++j);

    if (best.size()<20)
      best.add(l);
    for (int tj=best.size()-1; tj>j; --tj)
      best[tj]=best[tj-1];
    best[j]=l;
  }
  if (idcount[best[0]]<=0) return;  // no hits in db found, skip query

  if (offset+s.seqlen<offset){
    offset=0;
    kmerpos.init(MAXSIZE,-1);
    kmerposrev.init(MAXSIZE,-1);
  }
  setkmerpos(kmerpos,s,offset);
  setkmerpos(kmerposrev,srev,offset);
  
  for (int l=0; l<best.size(); ++l){
    int sbest=(seqids[best[l]]<seqs.size()?seqids[best[l]]:seqids[best[l]]-seqs.size());
    if (offset2+seqs.values(sbest).seqlen<offset2){
      offset2=0;
      kmerpos2.init(MAXSIZE,-1);
    }
    setkmerpos(kmerpos2,seqs.values(sbest),offset2);

    ealigndata adata;
    adata.seqid=sbest;
    adata.revcompl=(seqids[best[l]]>=seqs.size());
    adata.kmercount=idcount[best[l]];
    if (seqids[best[l]]<seqs.size())
      seqident_local(s,offset,kmerpos,seqs.values(sbest),offset2,kmerpos2,adata); //matches,mismatches,gaps);
    else
      seqident_local(srev,offset,kmerposrev,seqs.values(sbest),offset2,kmerpos2,adata); //matches,mismatches,gaps);
    matchcounts.add(adata);
    offset2+=seqs.values(sbest).seqlen;
  }
  offset+=s.seqlen;

  heapsort(matchcounts);
}
*/

class epredinfotax {
 public:
  eintarray tax;
  efloatarray cf;
  efloatarray cutoff;
};

class epredinfo {
 public:
  int seqid;
  ealigndata tophit;
  ebasicarray<ealigndata> matchcounts;
//  earrayof<float,int> ptax;
  ebasicarray<epredinfotax> ptax;
  efloatarray cutoff;
  epredinfo(const ealignscore& as): tophit(as) {}
  void clear() { matchcounts.clear(); tophit.seqid=-1; }
};


/*
void calc_score_bs(unsigned char* tncounts,int p,int e,edoublearray &scores,unsigned char *bscoefs,int bscount)
{
//  uint64_t *pcount=reinterpret_cast<uint64_t*>(tncounts);
  uint64_t vid;
  scores.init(scores.size(),0.0l);
  double tmp;
  char tmpstr[32];

  for (int i=p; i<e; ++i){
    vid=tncounts[i];
    switch (vid&0x3u){ 
      case 0x0u: tmp=gapcost; break;
      case 0x1u: tmp=matchcost; break;
      case 0x2u: tmp=misscost; break;
    }
    for (int j=0; j<scores.size(); ++j) scores[j]+=double(bscoefs[j+i*bscount])*tmp;
  }
}
*/

const int SEQSEGSIZE=1500;

void seqsearch_reponly(const estr& str2id,eseqdb& db,eseq& s,earray<epredinfo>& pinfoarr,esearchws& sws)
{
  t1.reset();
  pinfoarr.clear();
//  memset(sws.bitmask,0,((db.otus.size()*2)/64+1)*sizeof(uint64_t));
//  kmercount_both(db.otus.size(),db.otukmers,s,sws.idcount2,sws.bitmask,akmers,0x0Fu);
//  kmercount_both2(db.otus.size(),db.otukmers,s,sws.idcount,sws.otukmerpos,akmers,0x0Fu);
//  kmercount_both_nopos2(db.otus.size(),db.otukmers,s,sws.idcount,sws.otukmerpos,akmers,0x0Fu);
  for (long segi=0; segi<s.seqlen/SEQSEGSIZE+1; ++segi){
    epredinfo pinfo(as);
    long s_start=segi*SEQSEGSIZE;
    long s_end=(s_start+SEQSEGSIZE<s.seqlen?s_start+SEQSEGSIZE:s.seqlen);
    long s_len=s_end-s_start;
    eseq srev;
    srev.setrevcompl(s,s_start,s_end);
//    cout << "# " << str2id << " start: " << s_start << " end: " << s_end << " len: " << s_len << " seqlen: " << s.seqlen << " srev: " << srev.seqlen << endl;
  
    sws.otukmerpos.init(db.otus.size()*2,0);
    sws.idcount.init(db.otus.size()*2,0);

//    cout << "# counting" << endl;
    kmercount_both_nopos2_skip(db.otus.size(),db.otukmers,s,s_start,s_end,sws.idcount,akmers,0x0Fu,s_len/100+1);
    ti=ti*0.99+t1.lap()*0.01;
  
    eintarray best;
  //  ebasicarray<ealigndata> matchcounts;
  
    // choosing sequences for kmercounting step
  //  eintarray& best(sws.best);
  //  eintarray bestcount;
  //  best.clear();
    int l;
    for (l=0; l<sws.idcount.size(); ++l){
  //    uint16_t bt=(sws.bitmask[l/64u]>>(l%64u))&0x1ul;
      if (sws.idcount[l]>0) { best.add(l); break; }
    }
    for (; l<sws.idcount.size(); ++l){
  //    uint16_t bt=(sws.bitmask[l/64u]>>(l%64u))&0x1ul;
      if (sws.idcount[l]==0 || (sws.idcount[l]<sws.idcount[best[best.size()-1]] && best.size()==3*topotus)) continue; // worst id than the bottom of the list
  
      int j;
      for (j=0; j<best.size() && sws.idcount[l]<sws.idcount[best[j]]; ++j);
  
      if (best.size()<3*topotus)
        best.add(l);
      for (int tj=best.size()-1; tj>j; --tj)
        best[tj]=best[tj-1];
      best[j]=l;
    }
    if (best.size()==0) { // no overlap found to any seq group representatives
      continue;
//      pinfo.seqid=-2;
//      return;
    }
    ts=ts*0.99+t1.lap()*0.01;
  
    eintarray seqids;
    for (int l=0; l<best.size() && l<topotus; ++l){
      int ibest=best[l];
      if (ibest<db.otus.size()){
  //      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest][0]) << "\t" << idcount[ibest] << endl;
        seqids.add(db.otus[ibest][0]);
      }else{
  //      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest-otus.size()][0]) << "\t" << idcount[ibest] << " reversed" << endl;
        ibest-=db.otus.size();
        seqids.add(db.seqs.size()+db.otus[ibest][0]);
      }
    }
  
//    cout << "# 2nd counting" << endl;

    randomize(seqids);
    sws.idcount.init(seqids.size(),0);
  
    if (sws.maskid+1<sws.maskid){
      sws.maskid=-1;
      sws.kmermask.init(KMERSIZE,-1);
    }
    ++sws.maskid;
  
    ti2=ti2*0.99+t1.lap()*0.01;
  //  memset(sws.kmerbitmask,0,(MAXSIZE2/64+1)*sizeof(uint64_t));
  //  setkmermask(sws.kmerbitmask,s,akmers,0xFul);
  //  kmercount_mask(db.seqs,seqids,sws.kmerbitmask,sws.maskid,sws.idcount,akmers,0xFul);
//    cout << "# 2nd counting -- setkmermask" << endl;
    setkmermask(sws.kmermask,s,sws.maskid,akmers,0xFul,s_start,s_end);
//    cout << "# 2nd counting -- kmercount_mask" << endl;
    kmercount_mask(db.seqs,seqids,sws.kmermask,sws.maskid,sws.idcount);
    ts2=ts2*0.99+t1.lap()*0.01;
  
  /*
    int ibest=0;
    for (int l=1; l<sws.idcount.size(); ++l)
      if (sws.idcount[l]>sws.idcount[ibest]) ibest=l;
  
    if (sws.idcount[ibest]<1) {
      pinfo.seqid=-2;
      return;  // no hits in db found, skip query
    }
  */
  
    best.clear();
    for (l=0; l<sws.idcount.size(); ++l){
      if (sws.idcount[l]>0){ best.add(l); break; }
    }
    for (; l<sws.idcount.size(); ++l){
  //    if (idcount[l]<idcount[best[best.size()-1]] && (idcount[l]<0.8*idcount[ibest] && best.size()>=20 || best.size()==100)) continue; // worse id than the bottom of the list
      if (sws.idcount[l]==0 || (sws.idcount[l]<sws.idcount[best[best.size()-1]] && best.size()>=tophits)) continue; // zero counts or worse id than the bottom of the list
  
      int j;
      for (j=0; j<best.size() && sws.idcount[l]<sws.idcount[best[j]]; ++j);
  
  //    if (best.size()<20 || idcount[l]<0.8*idcount[ibest] && best.size()<100)
      if (best.size()<tophits)
        best.add(l);
      for (int tj=best.size()-1; tj>j; --tj)
        best[tj]=best[tj-1];
      best[j]=l;
    }

    if (best.size()==0){
//      pinfo.seqid=-2;
//      return;  // no hits in db found, skip query
      continue;
    }
  
    if (sws.offset+int(s_len)<sws.offset){ // need an int here otherwise the comparison is made in long and the offset is not correctly reset
      sws.offset=0;
      sws.kmerpos.init(MAXSIZE,-1);
      sws.kmerposrev.init(MAXSIZE,-1);
    }

//    cout << "# 2nd counting -- setkmerpos" << endl;
    setkmerpos(sws.kmerpos,s,sws.offset,s_start,s_end);
//    cout << "# 2nd counting -- setkmerposrev" << endl;
    setkmerpos(sws.kmerposrev,srev,sws.offset);
   
  /*
    // add last added seqs (worst kmercounts) to improve confidence estimation by keeping lower scoring seqs for full alignment
    // TODO: might not need full alignment but just estimated alignment score from shared kmers
    eintarray chosenflag;
    chosenflag.init(seqids.size(),0);
  
    for (int l=0; l<best.size(); ++l) chosenflag[best[l]]=1;
    for (int l=0,c=seqids.size()-1; c>=0 && l<10; --c){
      if (chosenflag[c]==1 || sws.idcount[c]==0) continue; // no shared kmer or already in list
      chosenflag[c]=1;
      best.add(c);
      ++l;
    }
  */
  
  //  eintarray seqboth;
  //  seqboth.init(db.seqs.size(),-1);
//    cout << "# aligning" << endl;
  
    for (int l=0; l<best.size(); ++l){
      int sbest=seqids[best[l]]%db.seqs.size();
      eseq &sdb(db.seqs.values(sbest));
      if (sws.offset2+int(db.seqs.values(sbest).seqlen)<sws.offset2){  // need an int here otherwise the comparison is made in long and the offset is not correctly reset
        sws.offset2=0;
        sws.kmerpos2.init(MAXSIZE,-1);
      }
      setkmerpos(sws.kmerpos2,db.seqs.values(sbest),sws.offset2);
  
      ealigndata adata(as);
      adata.seqid=sbest;
      adata.revcompl=(seqids[best[l]]>=db.seqs.size());
      adata.kmercount=sws.idcount[best[l]];
    
      if (seqids[best[l]]<db.seqs.size()){
//        cout << "# forward align" << endl;
        seqident_local(s,sws.kmerpos,sdb,adata,sws,as,s_start,s_end);
  //      seqident_global(s,sws.kmerpos,sdb,adata,sws,as);
  //      seqcalign_global_noedgegap(s,0,s.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
  /*
        seqcalign_global_noedgegap(s,0,s.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
        seqident_local(s,sws.kmerpos,sdb,adata2,sws,as);
  
        cout << str2id << endl;
        cout << "# " << " M:"<<adata.matches << " X:" << adata.mismatches << " G:"<<adata.gaps << " S:"<< adata._score << " "<< adata.profile << endl;
        if (!(adata==adata2))
          cout << "#!" << " M:"<<adata2.matches << " X:" << adata2.mismatches << " G:"<<adata2.gaps << " S:"<<adata2._score << " K:" << adata2.kmercount << " " << adata2.profile << endl;
  */
      }else{
//        cout << "# reverse align" << endl;
        seqident_local(srev,sws.kmerposrev,sdb,adata,sws,as);
        // flip 
        int tmp=srev.seqlen-adata.s1+s_start; adata.s1=srev.seqlen-adata.e1+s_start; adata.e1=tmp; 
        adata.revcompl=true;
  //      seqident_global(srev,sws.kmerposrev,sdb,adata,sws,as);
  //      seqcalign_global_noedgegap(srev,0,srev.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
  /*
        seqcalign_global_noedgegap(srev,0,srev.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
        seqident_local(srev,sws.kmerposrev,sdb,adata2,sws,as);
  
        cout << str2id << endl;
        cout << "# " << " M:"<<adata.matches << " X:" << adata.mismatches << " G:"<<adata.gaps << " S:" <<adata._score << " " << adata.profile << endl;
        if (!(adata==adata2))
          cout << "#!" << " M:"<<adata2.matches << " X:" << adata2.mismatches << " G:"<<adata2.gaps << " S:" << adata2._score << " K:" << adata2.kmercount << " " << adata2.profile << endl;
  */
  
      }
  //    LDEBUG(D_SEQALIGNMENT,print_tncount(&tncounts[l*NCOUNT_MAXLEN],0,s.seqlen));
      if (adata.matches+adata.mismatches>0 && adata.score()>=30)
        pinfo.matchcounts.add(adata);
  
      sws.offset2+=db.seqs.values(sbest).seqlen;
    }
    sws.offset+=s_len;
 
    if (pinfo.matchcounts.size()>0){
      heapsort(pinfo.matchcounts);
    //  for (int i=0; i<pinfo.matchcounts.size(); ++i)
    //    cout << "#best: " << i << " " << seqs.keys(pinfo.matchcounts[i].seqid) << " " << pinfo.matchcounts[i].score() << endl;
      pinfo.tophit=pinfo.matchcounts[pinfo.matchcounts.size()-1];
      pinfo.seqid=pinfo.tophit.seqid;
      pinfoarr.add(pinfo);
    }
  }
  ta=ta*0.99+t1.lap()*0.01;
}

class eidavg
{
 public:
  int count;
  float sum;
  float sums;
  void add(float v,float s) { ++count; sum+=v; sums+=s; }
  float avg_id() { return(count>0?sum/count:0.0); }
  float avg_score() { return(count>0?sums/count:0.0); }
  eidavg(): count(0),sum(0.0),sums(0.0) {}
};

void seqsearch2(const estr& str2id,eseqdb& db,eseq& s,earray<epredinfo>& pinfoarr,esearchws& sws)
{
  t1.reset();
  pinfoarr.clear();
//  memset(sws.bitmask,0,((db.otus.size()*2)/64+1)*sizeof(uint64_t));
//  kmercount_both(db.otus.size(),db.otukmers,s,sws.idcount2,sws.bitmask,akmers,0x0Fu);
//  kmercount_both2(db.otus.size(),db.otukmers,s,sws.idcount,sws.otukmerpos,akmers,0x0Fu);
//  kmercount_both_nopos2(db.otus.size(),db.otukmers,s,sws.idcount,sws.otukmerpos,akmers,0x0Fu);
  for (long segi=0; segi<s.seqlen/SEQSEGSIZE+1; ++segi){
    epredinfo pinfo(as);
    long s_start=segi*SEQSEGSIZE;
    long s_end=(s_start+SEQSEGSIZE<s.seqlen?s_start+SEQSEGSIZE:s.seqlen);
    long s_len=s_end-s_start;
    eseq srev;
    srev.setrevcompl(s,s_start,s_end);
//    cout << "# " << str2id << " start: " << s_start << " end: " << s_end << " len: " << s_len << " seqlen: " << s.seqlen << " srev: " << srev.seqlen << endl;
    t1.reset();
  
    sws.otukmerpos.init(db.otus.size()*2,0);
    sws.idcount.init(db.otus.size()*2,0);

//    cout << "# counting" << endl;
    kmercount_both_nopos2_skip(db.otus.size(),db.otukmers,s,s_start,s_end,sws.idcount,akmers,0x0Fu,s_len/100+1);
    ti=ti*0.99+t1.lap()*0.01;
  
    eintarray best;
  //  ebasicarray<ealigndata> matchcounts;
  
    // choosing sequences for kmercounting step
  //  eintarray& best(sws.best);
  //  eintarray bestcount;
  //  best.clear();
    int l;
    for (l=0; l<sws.idcount.size(); ++l){
  //    uint16_t bt=(sws.bitmask[l/64u]>>(l%64u))&0x1ul;
      if (sws.idcount[l]>0) { best.add(l); break; }
    }
    for (; l<sws.idcount.size(); ++l){
  //    uint16_t bt=(sws.bitmask[l/64u]>>(l%64u))&0x1ul;
      if (sws.idcount[l]==0 || (sws.idcount[l]<sws.idcount[best[best.size()-1]] && best.size()==3*topotus)) continue; // worst id than the bottom of the list
  
      int j;
      for (j=0; j<best.size() && sws.idcount[l]<sws.idcount[best[j]]; ++j);
  
      if (best.size()<3*topotus)
        best.add(l);
      for (int tj=best.size()-1; tj>j; --tj)
        best[tj]=best[tj-1];
      best[j]=l;
    }
    if (best.size()==0) { // no overlap found to any seq group representatives
      continue;
//      pinfo.seqid=-2;
//      return;
    }
    ts=ts*0.99+t1.lap()*0.01;
  
    eintarray seqids;
    for (int l=0; l<best.size() && l<topotus; ++l){
      int ibest=best[l];
      if (ibest<db.otus.size()){
  //      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest][0]) << "\t" << idcount[ibest] << endl;
        for (int l2=0; l2<db.otus[ibest].size() && l2<20; ++l2)
          seqids.add(db.otus[ibest][l2]);
      }else{
  //      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest-otus.size()][0]) << "\t" << idcount[ibest] << " reversed" << endl;
        ibest-=db.otus.size();
        for (int l2=0; l2<db.otus[ibest].size() && l2<20; ++l2)
          seqids.add(db.seqs.size()+db.otus[ibest][l2]);
      }
    }
  
    // add representatives from all non-chosen OTUS (may improve confidence and novel otu estimation)
  /*
    for (int l=topotus; l<best.size(); ++l){
      int ibest=best[l];
      if (ibest<db.otus.size()){
        if (db.otus[ibest].size()==0 || sws.idcount[ibest]==0) continue;
        seqids.add(db.otus[ibest][0]);
      } else {
  //      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest-otus.size()][0]) << "\t" << idcount[ibest] << " reversed" << endl;
        ibest-=db.otus.size();
        if (db.otus[ibest].size()==0 || sws.idcount[ibest]==0) continue;
        seqids.add(db.seqs.size()+db.otus[ibest][0]);
      }
    }
  */
  
//    cout << "# 2nd counting" << endl;

    randomize(seqids);
    sws.idcount.init(seqids.size(),0);
  
    if (sws.maskid+1<sws.maskid){
      sws.maskid=-1;
      sws.kmermask.init(KMERSIZE,-1);
    }
    ++sws.maskid;
  
    ti2=ti2*0.99+t1.lap()*0.01;
  //  memset(sws.kmerbitmask,0,(MAXSIZE2/64+1)*sizeof(uint64_t));
  //  setkmermask(sws.kmerbitmask,s,akmers,0xFul);
  //  kmercount_mask(db.seqs,seqids,sws.kmerbitmask,sws.maskid,sws.idcount,akmers,0xFul);
//    cout << "# 2nd counting -- setkmermask" << endl;
    setkmermask(sws.kmermask,s,sws.maskid,akmers,0xFul,s_start,s_end);
//    cout << "# 2nd counting -- kmercount_mask" << endl;
    kmercount_mask(db.seqs,seqids,sws.kmermask,sws.maskid,sws.idcount);
    ts2=ts2*0.99+t1.lap()*0.01;
  
  /*
    int ibest=0;
    for (int l=1; l<sws.idcount.size(); ++l)
      if (sws.idcount[l]>sws.idcount[ibest]) ibest=l;
  
    if (sws.idcount[ibest]<1) {
      pinfo.seqid=-2;
      return;  // no hits in db found, skip query
    }
  */
  
    best.clear();
    for (l=0; l<sws.idcount.size(); ++l){
      if (sws.idcount[l]>0){ best.add(l); break; }
    }
    for (; l<sws.idcount.size(); ++l){
  //    if (idcount[l]<idcount[best[best.size()-1]] && (idcount[l]<0.8*idcount[ibest] && best.size()>=20 || best.size()==100)) continue; // worse id than the bottom of the list
      if (sws.idcount[l]==0 || (sws.idcount[l]<sws.idcount[best[best.size()-1]] && best.size()>=tophits)) continue; // zero counts or worse id than the bottom of the list
  
      int j;
      for (j=0; j<best.size() && sws.idcount[l]<sws.idcount[best[j]]; ++j);
  
  //    if (best.size()<20 || idcount[l]<0.8*idcount[ibest] && best.size()<100)
      if (best.size()<tophits)
        best.add(l);
      for (int tj=best.size()-1; tj>j; --tj)
        best[tj]=best[tj-1];
      best[j]=l;
    }
    if (best.size()==0){
//      pinfo.seqid=-2;
//      return;  // no hits in db found, skip query
      continue;
    }
  
    if (sws.offset+int(s_len)<sws.offset){ // need an int here otherwise the comparison is made in long and the offset is not correctly reset
      sws.offset=0;
      sws.kmerpos.init(MAXSIZE,-1);
      sws.kmerposrev.init(MAXSIZE,-1);
    }
//    cout << "# 2nd counting -- setkmerpos" << endl;
    setkmerpos(sws.kmerpos,s,sws.offset,s_start,s_end);
//    cout << "# 2nd counting -- setkmerposrev" << endl;
    setkmerpos(sws.kmerposrev,srev,sws.offset);
   
  /*
    // add last added seqs (worst kmercounts) to improve confidence estimation by keeping lower scoring seqs for full alignment
    // TODO: might not need full alignment but just estimated alignment score from shared kmers
    eintarray chosenflag;
    chosenflag.init(seqids.size(),0);
  
    for (int l=0; l<best.size(); ++l) chosenflag[best[l]]=1;
    for (int l=0,c=seqids.size()-1; c>=0 && l<10; --c){
      if (chosenflag[c]==1 || sws.idcount[c]==0) continue; // no shared kmer or already in list
      chosenflag[c]=1;
      best.add(c);
      ++l;
    }
  */
  
  //  eintarray seqboth;
  //  seqboth.init(db.seqs.size(),-1);
//    cout << "# aligning" << endl;
  
    for (int l=0; l<best.size(); ++l){
      int sbest=seqids[best[l]]%db.seqs.size();
      eseq &sdb(db.seqs.values(sbest));
      if (sws.offset2+int(db.seqs.values(sbest).seqlen)<sws.offset2){  // need an int here otherwise the comparison is made in long and the offset is not correctly reset
        sws.offset2=0;
        sws.kmerpos2.init(MAXSIZE,-1);
      }
      setkmerpos(sws.kmerpos2,db.seqs.values(sbest),sws.offset2);
  
      ealigndata adata(as);
      adata.seqid=sbest;
      adata.revcompl=(seqids[best[l]]>=db.seqs.size());
      adata.kmercount=sws.idcount[best[l]];
    
      if (seqids[best[l]]<db.seqs.size()){
//        cout << "# forward align" << endl;
        seqident_local(s,sws.kmerpos,sdb,adata,sws,as,s_start,s_end);
  //      seqident_global(s,sws.kmerpos,sdb,adata,sws,as);
  //      seqcalign_global_noedgegap(s,0,s.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
  /*
        seqcalign_global_noedgegap(s,0,s.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
        seqident_local(s,sws.kmerpos,sdb,adata2,sws,as);
  
        cout << str2id << endl;
        cout << "# " << " M:"<<adata.matches << " X:" << adata.mismatches << " G:"<<adata.gaps << " S:"<< adata._score << " "<< adata.profile << endl;
        if (!(adata==adata2))
          cout << "#!" << " M:"<<adata2.matches << " X:" << adata2.mismatches << " G:"<<adata2.gaps << " S:"<<adata2._score << " K:" << adata2.kmercount << " " << adata2.profile << endl;
  */
      }else{
//        cout << "# reverse align" << endl;
        seqident_local(srev,sws.kmerposrev,sdb,adata,sws,as);
        // flip 
        int tmp=srev.seqlen-adata.s1+s_start; adata.s1=srev.seqlen-adata.e1+s_start; adata.e1=tmp; 
        adata.revcompl=true;
  //      seqident_global(srev,sws.kmerposrev,sdb,adata,sws,as);
  //      seqcalign_global_noedgegap(srev,0,srev.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
  /*
        seqcalign_global_noedgegap(srev,0,srev.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
        seqident_local(srev,sws.kmerposrev,sdb,adata2,sws,as);
  
        cout << str2id << endl;
        cout << "# " << " M:"<<adata.matches << " X:" << adata.mismatches << " G:"<<adata.gaps << " S:" <<adata._score << " " << adata.profile << endl;
        if (!(adata==adata2))
          cout << "#!" << " M:"<<adata2.matches << " X:" << adata2.mismatches << " G:"<<adata2.gaps << " S:" << adata2._score << " K:" << adata2.kmercount << " " << adata2.profile << endl;
  */
  
      }
  //    LDEBUG(D_SEQALIGNMENT,print_tncount(&tncounts[l*NCOUNT_MAXLEN],0,s.seqlen));
      if (adata.matches+adata.mismatches>0 && adata.score()>=30)
        pinfo.matchcounts.add(adata);
  
      sws.offset2+=db.seqs.values(sbest).seqlen;
    }
    sws.offset+=s_len;
 
    if (pinfo.matchcounts.size()>0){
      heapsort(pinfo.matchcounts);
    //  for (int i=0; i<pinfo.matchcounts.size(); ++i)
    //    cout << "#best: " << i << " " << seqs.keys(pinfo.matchcounts[i].seqid) << " " << pinfo.matchcounts[i].score() << endl;
      pinfo.tophit=pinfo.matchcounts[pinfo.matchcounts.size()-1];
      pinfo.seqid=pinfo.tophit.seqid;
      pinfoarr.add(pinfo);
    }
    ta=ta*0.99+t1.lap()*0.01;
  }
}

void seqsearch(const estr& str2id,eseqdb& db,eseq& s,earray<epredinfo>& pinfoarr,esearchws& sws)
{
  t1.reset();
  pinfoarr.clear();
//  memset(sws.bitmask,0,((db.otus.size()*2)/64+1)*sizeof(uint64_t));
//  kmercount_both(db.otus.size(),db.otukmers,s,sws.idcount2,sws.bitmask,akmers,0x0Fu);
//  kmercount_both2(db.otus.size(),db.otukmers,s,sws.idcount,sws.otukmerpos,akmers,0x0Fu);
//  kmercount_both_nopos2(db.otus.size(),db.otukmers,s,sws.idcount,sws.otukmerpos,akmers,0x0Fu);
  for (long segi=0; segi<s.seqlen/SEQSEGSIZE+1; ++segi){
    epredinfo pinfo(as);
    long s_start=segi*SEQSEGSIZE;
    long s_end=(s_start+SEQSEGSIZE<s.seqlen?s_start+SEQSEGSIZE:s.seqlen);
    long s_len=s_end-s_start;
    eseq srev;
    srev.setrevcompl(s,s_start,s_end);
//    cout << "# " << str2id << " start: " << s_start << " end: " << s_end << " len: " << s_len << " seqlen: " << s.seqlen << " srev: " << srev.seqlen << endl;
  
//    sws.otukmerpos.init(db.otus.size()*2,0);
//    sws.idcount.init(db.otus.size()*2,0);
    t1.reset();

    uint32_t *plist=&sws.idlist[0];
//    cout << "# counting" << endl;
    kmercount_both_nopos2_skip_plist(db.otus.size(),db.otukmers,s,s_start,s_end,sws.idcount,plist,akmers,0x0Fu,s_len/100+1);
    ti=ti*0.99+t1.lap()*0.01;

//    cout << plist-&sws.idlist[0] << endl;
  
    eintarray best;
    eintarray bestcount;
  //  ebasicarray<ealigndata> matchcounts;
  
    // choosing sequences for kmercounting step
  //  eintarray& best(sws.best);
  //  eintarray bestcount;
  //  best.clear();
    int l;
    uint32_t *pblist;
    for (pblist=&sws.idlist[0]; pblist!=plist; ++pblist){
      if (best.size()!=0 && sws.idcount[*pblist]<bestcount[best.size()-1] && best.size()==3*topotus) {
        sws.idcount[*pblist]=0u;
        continue; // worst id than the bottom of the list
      }
  
      int j;
      for (j=0; j<best.size() && sws.idcount[*pblist]<bestcount[j]; ++j);
  
      if (best.size()<3*topotus){
        best.add(*pblist);
        bestcount.add(sws.idcount[*pblist]);
      }
      for (int tj=best.size()-1; tj>j; --tj){
        best[tj]=best[tj-1];
        bestcount[tj]=bestcount[tj-1];
      }
      best[j]=*pblist;
      bestcount[j]=sws.idcount[*pblist];
      sws.idcount[*pblist]=0u;
    }
    if (best.size()==0) { // no overlap found to any seq group representatives
      continue;
//      pinfo.seqid=-2;
//      return;
    }
    ts=ts*0.99+t1.lap()*0.01;
  
    eintarray seqids;
    for (int l=0; l<best.size() && l<topotus; ++l){
      int ibest=best[l];
      if (ibest<db.otus.size()){
  //      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest][0]) << "\t" << idcount[ibest] << endl;
        for (int l2=0; l2<db.otus[ibest].size() && l2<otulim; ++l2)
          seqids.add(db.otus[ibest][l2]);
      }else{
  //      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest-otus.size()][0]) << "\t" << idcount[ibest] << " reversed" << endl;
        ibest-=db.otus.size();
        for (int l2=0; l2<db.otus[ibest].size() && l2<otulim; ++l2)
          seqids.add(db.seqs.size()+db.otus[ibest][l2]);
      }
    }
  
    // add representatives from all non-chosen OTUS (may improve confidence and novel otu estimation)
  /*
    for (int l=topotus; l<best.size(); ++l){
      int ibest=best[l];
      if (ibest<db.otus.size()){
        if (db.otus[ibest].size()==0 || sws.idcount[ibest]==0) continue;
        seqids.add(db.otus[ibest][0]);
      } else {
  //      cout << "#\t" << str2id << "\t" << seqs.keys(otus[ibest-otus.size()][0]) << "\t" << idcount[ibest] << " reversed" << endl;
        ibest-=db.otus.size();
        if (db.otus[ibest].size()==0 || sws.idcount[ibest]==0) continue;
        seqids.add(db.seqs.size()+db.otus[ibest][0]);
      }
    }
  */
  
//    cout << "# 2nd counting" << endl;

    randomize(seqids);
  
    if (sws.maskid+1<sws.maskid){
      sws.maskid=-1;
      sws.kmermask.init(KMERSIZE,-1);
    }
    ++sws.maskid;
  
    ti2=ti2*0.99+t1.lap()*0.01;
  //  memset(sws.kmerbitmask,0,(MAXSIZE2/64+1)*sizeof(uint64_t));
  //  setkmermask(sws.kmerbitmask,s,akmers,0xFul);
  //  kmercount_mask(db.seqs,seqids,sws.kmerbitmask,sws.maskid,sws.idcount,akmers,0xFul);
//    cout << "# 2nd counting -- setkmermask" << endl;
    setkmermask(sws.kmermask,s,sws.maskid,akmers,0xFul,s_start,s_end);
//    cout << "# 2nd counting -- kmercount_mask" << endl;
    kmercount_mask(db.seqs,seqids,sws.kmermask,sws.maskid,sws.idcount);
    ts2=ts2*0.99+t1.lap()*0.01;
 
    best.clear();
    for (l=0; l<seqids.size(); ++l){
      if (sws.idcount[l]>0){ best.add(l); break; }
    }
    for (; l<seqids.size(); ++l){
  //    if (idcount[l]<idcount[best[best.size()-1]] && (idcount[l]<0.8*idcount[ibest] && best.size()>=20 || best.size()==100)) continue; // worse id than the bottom of the list
      if (sws.idcount[l]==0 || (sws.idcount[l]<sws.idcount[best[best.size()-1]] && best.size()>=tophits)) continue; // zero counts or worse id than the bottom of the list
  
      int j;
      for (j=0; j<best.size() && sws.idcount[l]<sws.idcount[best[j]]; ++j);
  
  //    if (best.size()<20 || idcount[l]<0.8*idcount[ibest] && best.size()<100)
      if (best.size()<tophits)
        best.add(l);
      for (int tj=best.size()-1; tj>j; --tj)
        best[tj]=best[tj-1];
      best[j]=l;
    }
    for (l=0; l<seqids.size(); ++l)
      sws.idcount[l]=0u;

    if (best.size()==0){
//      pinfo.seqid=-2;
//      return;  // no hits in db found, skip query
      continue;
    }
  
    if (sws.offset+int(s_len)<sws.offset){ // need an int here otherwise the comparison is made in long and the offset is not correctly reset
      sws.offset=0;
      sws.kmerpos.init(MAXSIZE,-1);
      sws.kmerposrev.init(MAXSIZE,-1);
    }
//    cout << "# 2nd counting -- setkmerpos" << endl;
    setkmerpos(sws.kmerpos,s,sws.offset,s_start,s_end);
//    cout << "# 2nd counting -- setkmerposrev" << endl;
    setkmerpos(sws.kmerposrev,srev,sws.offset);
   
 
    for (int l=0; l<best.size(); ++l){
      int sbest=seqids[best[l]]%db.seqs.size();
      eseq &sdb(db.seqs.values(sbest));
      if (sws.offset2+int(db.seqs.values(sbest).seqlen)<sws.offset2){  // need an int here otherwise the comparison is made in long and the offset is not correctly reset
        sws.offset2=0;
        sws.kmerpos2.init(MAXSIZE,-1);
      }
      setkmerpos(sws.kmerpos2,db.seqs.values(sbest),sws.offset2);
  
      ealigndata adata(as);
      adata.seqid=sbest;
      adata.revcompl=(seqids[best[l]]>=db.seqs.size());
      adata.kmercount=sws.idcount[best[l]];
    
      if (seqids[best[l]]<db.seqs.size()){
//        cout << "# forward align" << endl;
        seqident_local(s,sws.kmerpos,sdb,adata,sws,as,s_start,s_end);
  //      seqident_global(s,sws.kmerpos,sdb,adata,sws,as);
  //      seqcalign_global_noedgegap(s,0,s.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
  /*
        seqcalign_global_noedgegap(s,0,s.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
        seqident_local(s,sws.kmerpos,sdb,adata2,sws,as);
  
        cout << str2id << endl;
        cout << "# " << " M:"<<adata.matches << " X:" << adata.mismatches << " G:"<<adata.gaps << " S:"<< adata._score << " "<< adata.profile << endl;
        if (!(adata==adata2))
          cout << "#!" << " M:"<<adata2.matches << " X:" << adata2.mismatches << " G:"<<adata2.gaps << " S:"<<adata2._score << " K:" << adata2.kmercount << " " << adata2.profile << endl;
  */
      }else{
//        cout << "# reverse align" << endl;
        seqident_local(srev,sws.kmerposrev,sdb,adata,sws,as);
        // flip 
        int tmp=srev.seqlen-adata.s1+s_start; adata.s1=srev.seqlen-adata.e1+s_start; adata.e1=tmp; 
        adata.revcompl=true;
  //      seqident_global(srev,sws.kmerposrev,sdb,adata,sws,as);
  //      seqcalign_global_noedgegap(srev,0,srev.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
  /*
        seqcalign_global_noedgegap(srev,0,srev.seqlen,sdb,0,sdb.seqlen,adata,sws.alignws,as);
        seqident_local(srev,sws.kmerposrev,sdb,adata2,sws,as);
  
        cout << str2id << endl;
        cout << "# " << " M:"<<adata.matches << " X:" << adata.mismatches << " G:"<<adata.gaps << " S:" <<adata._score << " " << adata.profile << endl;
        if (!(adata==adata2))
          cout << "#!" << " M:"<<adata2.matches << " X:" << adata2.mismatches << " G:"<<adata2.gaps << " S:" << adata2._score << " K:" << adata2.kmercount << " " << adata2.profile << endl;
  */
  
      }
  //    LDEBUG(D_SEQALIGNMENT,print_tncount(&tncounts[l*NCOUNT_MAXLEN],0,s.seqlen));
      if (adata.matches+adata.mismatches>0 && adata.score()>=30)
        pinfo.matchcounts.add(adata);
  
      sws.offset2+=db.seqs.values(sbest).seqlen;
    }
    sws.offset+=s_len;
 
    if (pinfo.matchcounts.size()>0){
      heapsort(pinfo.matchcounts);
    //  for (int i=0; i<pinfo.matchcounts.size(); ++i)
    //    cout << "#best: " << i << " " << seqs.keys(pinfo.matchcounts[i].seqid) << " " << pinfo.matchcounts[i].score() << endl;
      pinfo.tophit=pinfo.matchcounts[pinfo.matchcounts.size()-1];
      int topi=-1;
      for (int t=pinfo.matchcounts.size()-2; t>=0; --t){
        if (pinfo.tophit.cmpoverlap(pinfo.matchcounts[t]))
          topi=t;
      }
      if (topi>0)
        pinfo.tophit=pinfo.matchcounts[topi];
      pinfo.seqid=pinfo.tophit.seqid;
      pinfoarr.add(pinfo);
    }
    ta=ta*0.99+t1.lap()*0.01;
  }
}

void findtaxcutoff(efloatarray& taxcutoff,epredinfo& pinfo2,earrayof<double,int>& predtax,etax& tax)
{
  ealigndata refhit=pinfo2.matchcounts[pinfo2.matchcounts.size()-1];
//      cout << "# ref: " << str2id << "\t" << seqs.keys(refhit.seqid) << "\t" << refhit.score() << endl;
  for (int l=pinfo2.matchcounts.size()-1; l>=0; --l){
    int sbest=pinfo2.matchcounts[l].seqid;
    if (tax.seqs[sbest]==0x00) continue;
    eseqtax &taxhit(*tax.seqs[sbest]);
    for (int k=0; k<taxhit.tl.size(); ++k){
      if (taxhit.tl[k].tid!=predtax.keys(k) && taxcutoff[k]==0.0) taxcutoff[k]=pinfo2.matchcounts[l].identity();
    }
  }
  // if computed thresholds are lower than what was predicted for the silver hit, then increase them to silver hit values
//        for (int k=0; k<tmptaxhit.tl.size(); ++k){
//          if (taxcutoff[k]>tmptaxhit.tl[k].ncf && tmptaxhit.tl[k].ncf>0.0) taxcutoff[k]=tmptaxhit.tl[k].ncf;
//        }
}

void taxscore(earrayof<double,int>& ptax,epredinfo& pinfo,etax& tax,ebasicarray<eintarray>& taxcounts)
{
  edoublearray taxscores;
  taxscores.init(tax.names.size(),0.0);
  for (int i=0; i<tax.names.size(); ++i)
    ptax.add(-1,0.0);

  int tophitl=-1;
  for (int l=pinfo.matchcounts.size()-1; l>=0; --l){
    int sbest=pinfo.matchcounts[l].seqid;
    if (tax.seqs[sbest]!=0x00) { tophitl=l; break; }
  }
  if (tophitl==-1) return;

  ealigndata &tophit(pinfo.matchcounts[tophitl]);
  pinfo.tophit=tophit; 
  ldieif(tax.seqs[tophit.seqid]==0x00,"top hit does not have taxonomy: "+estr(tophit.seqid));
  eseqtax &toptax(*tax.seqs[tophit.seqid]);

  double topscore=tophit.score();

  int taxid=0;
  taxcounts.init(tax.names.size());
  for (int i=0; i<taxcounts.size(); ++i)
    taxcounts[i].init(tax.names[i].size(),-1);
  for (int l=pinfo.matchcounts.size()-1; l>=0; --l){
    int sbest=pinfo.matchcounts[l].seqid;
    if (tax.seqs[sbest]==0x00) continue;
//    if (pinfo.matchcounts[l].score()<=0.0) break; // do not use alignments with less than or zero score
    eseqtax &taxhit(*tax.seqs[sbest]);
    for (int k=0; k<taxhit.tl.size(); ++k){
      if (taxcounts[k][taxhit.tl[k].tid]==taxid) continue;
      taxcounts[k][taxhit.tl[k].tid]=taxid;
      taxscores[k]+=exp((1.0l-topscore/pinfo.matchcounts[l].score())*20.0l);
//      taxscores[k]+=exp(log(pinfo.matchcounts[l].identity()*100.0)*30.0);
    }
/*
      if (seqids[best[imc[l]]]<seqs.size())
        cout << str2id << "\t" << seqs.keys(sbest) << "\t" << matchcounts[imc[l]] << "\t" << s.seqlen << endl;
      else
        cout << str2id << "\t" << seqs.keys(sbest) << "\t" << matchcounts[imc[l]] << "\t" << s.seqlen << "\treversed" << endl;
*/
  }

  
  for (int l=0; l<toptax.tl.size(); ++l){
    ldieif(toptax.tl[l].tid>=tax.names[l].size(),estr("key out of tax: ")+tophit.seqid+" "+estr(l)+" "+toptax.tl[l].tid+" "+tax.names[l].size());
    ptax.keys(l)=toptax.tl[l].tid;
    ptax.values(l)=1.0l/taxscores[l];
//    cout << "# i: " << tophit.identity() << " exp: " << exp(log(tophit.identity()*100.0l)*30.0l) << " sum: " << taxscores[l] << " = " << pinfo.predtax2.values(l) << endl;
//    pinfo.predtax2.add(toptax[l],exp(log(tophit.score())*30.0l)/taxscores[l]);
//    cout << "\t" << tax[l].values(toptax[l]) << "\t" << exp(log(tophit.score())*30.0l)/taxscores[l] << "\t" << (noveltaxa?estr(taxcutoff[l]/hitlen):estr("-"));
  }
}

void taxscore2(earrayof<double,int>& ptax,epredinfo& pinfo,etax& tax,ebasicarray<eintarray>& taxcounts)
{
  edoublearray taxscores;
  taxscores.init(tax.names.size(),0.0);
  for (int i=0; i<tax.names.size(); ++i)
    ptax.add(-1,0.0);

  int tophitl=-1;
  for (int l=pinfo.matchcounts.size()-1; l>=0; --l){
    int sbest=pinfo.matchcounts[l].seqid;
    if (tax.seqs[sbest]!=0x00) { tophitl=l; break; }
  }
  if (tophitl==-1) return;

  ealigndata &tophit(pinfo.matchcounts[tophitl]);
  pinfo.tophit=tophit; 
  ldieif(tax.seqs[tophit.seqid]==0x00,"top hit does not have taxonomy: "+estr(tophit.seqid));
  eseqtax &toptax(*tax.seqs[tophit.seqid]);

  int taxid=0;
  taxcounts.init(tax.names.size());
  for (int i=0; i<taxcounts.size(); ++i)
    taxcounts[i].init(tax.names[i].size(),-1);
  for (int l=pinfo.matchcounts.size()-1; l>=0; --l){
    int sbest=pinfo.matchcounts[l].seqid;
    if (tax.seqs[sbest]==0x00) continue;
//    if (pinfo.matchcounts[l].score()<=0.0) break; // do not use alignments with less than or zero score
    eseqtax &taxhit(*tax.seqs[sbest]);
    for (int k=0; k<taxhit.tl.size(); ++k){
      if (taxcounts[k][taxhit.tl[k].tid]==taxid) continue;
      taxcounts[k][taxhit.tl[k].tid]=taxid;
//      taxscores[k]+=exp(log(pinfo.matchcounts[l].score())*30.0);
      taxscores[k]+=exp(log(pinfo.matchcounts[l].identity()*100.0)*30.0);
    }
/*
      if (seqids[best[imc[l]]]<seqs.size())
        cout << str2id << "\t" << seqs.keys(sbest) << "\t" << matchcounts[imc[l]] << "\t" << s.seqlen << endl;
      else
        cout << str2id << "\t" << seqs.keys(sbest) << "\t" << matchcounts[imc[l]] << "\t" << s.seqlen << "\treversed" << endl;
*/
  }

  
  for (int l=0; l<toptax.tl.size(); ++l){
    ldieif(toptax.tl[l].tid>=tax.names[l].size(),estr("key out of tax: ")+tophit.seqid+" "+estr(l)+" "+toptax.tl[l].tid+" "+tax.names[l].size());
    ptax.keys(l)=toptax.tl[l].tid;
    ptax.values(l)=exp(log(tophit.identity()*100.0l)*30.0l)/taxscores[l];
//    cout << "# i: " << tophit.identity() << " exp: " << exp(log(tophit.identity()*100.0l)*30.0l) << " sum: " << taxscores[l] << " = " << pinfo.predtax2.values(l) << endl;
//    pinfo.predtax2.add(toptax[l],exp(log(tophit.score())*30.0l)/taxscores[l]);
//    cout << "\t" << tax[l].values(toptax[l]) << "\t" << exp(log(tophit.score())*30.0l)/taxscores[l] << "\t" << (noveltaxa?estr(taxcutoff[l]/hitlen):estr("-"));
  }
}

void load_taxa(const estr& taxfile,eseqdb& db)
{
  efile f;
  estr line;
  estrarray parts,parts2;
  int notfound=0;
  f.open(taxfile,"r");
  etax& tax(db.taxa.add(etax()));
  tax.seqs.init(db.seqs.size(),0x00);
  int taxind=0;
  while (!f.eof()){
    f.readln(line);
    if (line.len()==0) continue;
    if (line[0]=='#'){
      parts=line.explode(" ");
      if (parts[0]=="#cutoff:"){
        for (int i=1; i<parts.size(); ++i){
          parts2=parts[i].explode(":");
          ldieif(parts2.size()<2,"not enough parts on #cutoff line, i.e.: 0.97:0.02");
          tax.cutoff.add(parts2[0].f());
          tax.cutoffcoef.add(parts2[1].f());
        }
      }else if (parts[0]=="#name:" && parts.size()>1){
        tax.name=parts[1];
      }else if (parts[0]=="#levels:" && parts.size()>1){
        for (int i=1; i<parts.size(); ++i)
          tax.levels.add(parts[i]);
      }
      continue; 
    }
    parts=line.explode("\t");
    ldieif(parts.size()<2,"loading taxonomy, not enough fields in line: "+line);
    if (db.seqind.exists(parts[0])){
      if (parts.size()==2){ // simple taxonomy file
        parts2=parts[1].explode(";");
        eseqtax *newtax=new eseqtax();
        newtax->tl.reserve(parts2.size());
        for (int i=0; i<parts2.size(); ++i){
          if (i>=tax.ind.size() || i>=tax.names.size()){
            tax.names.add(estrarray());
            tax.ind.add(estrhashof<int>());
          }
          if (!tax.ind[i].exists(parts2[i])){
            tax.ind[i].add(parts2[i],tax.names[i].size());
            tax.names[i].add(parts2[i]);
          }
          newtax->tl.add(eseqtaxlevel(tax.ind[i][parts2[i]]));
        }
        tax.seqs[db.seqind[parts[0]]]=newtax;
      }else if (parts.size()>11) { // indirect taxonomy file with confidences
        
        if (taxind==0){
          for (int i=1; i<parts.size(); ++i)
            if (parts[i].len()==0)
              { taxind=i+1; break; }
        }
        int dbhit=(db.seqind.exists(parts[1])?db.seqind[parts[1]]:-1);
//        lerrorif(dbhit==-1,"gold hit not found: "+parts[1]);
        eseqtax *newtax=new eseqtax(parts[3].f());
        int taxfields=(parts.size()-taxind)/3;
        newtax->tl.reserve(taxfields);
//        cout << parts[0] << " " << parts[3] << " " << taxfields;
        for (int i=0; i<taxfields; ++i){
          int field=i*3+taxind;
          if (i>=tax.ind.size() || i>=tax.names.size()){
            tax.names.add(estrarray());
            tax.ind.add(estrhashof<int>());
          }
          if (!tax.ind[i].exists(parts[field])){
            tax.ind[i].add(parts[field],tax.names[i].size());
            tax.names[i].add(parts[field]);
          }
          newtax->tl.add(eseqtaxlevel(tax.ind[i][parts[field]],parts[field+2].f()));
//          cout << parts[field] << "(" << taxind[i][parts[field]] << ") " << parts[field+1] << " " << parts[field+2] << endl;
        }
//        cout << endl;
        tax.seqs[db.seqind[parts[0]]]=newtax;
      }
    }else{
      ++notfound;
    }
  }
  if (tax.name.len()==0)
    tax.name=basename(taxfile);
  lwarnif(notfound>0,"loading taxonomy, "+estr(notfound)+" sequences not found in sequence database");
  f.close();
}

void make_taxcluster(etax& tax,const estr& rank,eseqdb& db)
{
  int rid=tax.levels.find(rank);
  ldieif(rid==-1,"rank not found: "+rank);
  db.otus2.init(tax.names[rid].size());
  db.seqotu2.init(db.seqs.size(),-1);
//  eintarray tmpkmers;
//  tmpkmers.init(MAXSIZE,-1);
  int tmpi=0;
  for (long i=0; i<db.seqs.size(); ++i){
    ldieif(tax.seqs[i]==0x00,"sequence missing taxa: "+db.seqs.keys(i));
    int otuid=tax.seqs[i]->tl[rid].tid;
    db.seqotu2[i]=otuid;
    db.otus2[otuid].add(i);
/*
    if (db.otus2[otuid].size()==1){
      otukmeradd(db.otukmers,otuid,db.seqs.values(i),tmpkmers,tmpi,akmers,0xFul);
      ++tmpi;
    }
*/
  }
  ldieif(db.otus2.size()==0,"no clusters in database");
}

template <class T>
class ebasicmtqueue
{
 public:
  queue<T> q;
  emutex m;
  bool finished;
  econdsig qSignal;
  bool get(T& v);
  void put(T v);
  void finish();

  ebasicmtqueue();
};

template <class T>
void ebasicmtqueue<T>::finish()
{
  m.lock();
  finished=true;
  qSignal.signal();
  m.unlock();
}

template <class T>
ebasicmtqueue<T>::ebasicmtqueue(): finished(false) {}

template <class T>
bool ebasicmtqueue<T>::get(T& v)
{
  m.lock();
  while (q.size()==0 && !finished) qSignal.wait(m);
  if (q.size()==0 && finished) { m.unlock(); return(false);}
  v=q.front();
  q.pop();
  m.unlock();
  return(true);
}

template <class T>
void ebasicmtqueue<T>::put(T b)
{
  m.lock();
  q.push_back(b);
  qSignal.signal();
  m.unlock();
}


template <class T>
class emtqueue
{
 public:
  queue<T> full;
  stack<T> empty;
  emutex m;
  bool finished;
  econdsig fullSignal;
  econdsig emptySignal;
  T getEmpty();
  T getFull();
  emtqueue();

  void returnFull(T b);
  void returnEmpty(T b);

  void finish();
};

template <class T>
emtqueue<T>::emtqueue(): finished(false) {}

template <class T>
void emtqueue<T>::finish()
{
  m.lock();
  finished=true;
  fullSignal.signal();
  m.unlock();
}

template <class T>
T emtqueue<T>::getEmpty()
{
  m.lock();
  while (empty.size()==0 && !finished) emptySignal.wait(m);
  if (finished) { m.unlock(); return(0x00);}

  T pbuf=empty.top();
  empty.pop();
  m.unlock();
  return(pbuf);
}

template <class T>
T emtqueue<T>::getFull()
{
  m.lock();
  while (full.size()==0 && !finished) fullSignal.wait(m);
  if (full.size()==0 && finished) { m.unlock(); return(0x00);}

  T pbuf=full.front();
  full.pop();
  m.unlock();
  return(pbuf);
}

template <class T>
void emtqueue<T>::returnFull(T b)
{
  m.lock();
  full.push(b);
  fullSignal.signal();
  m.unlock();
}

template <class T>
void emtqueue<T>::returnEmpty(T b)
{
  m.lock();
  empty.push(b);
  emptySignal.signal();
  m.unlock();
}


struct emtdata {
  emtqueue<estrarrayof<eseq>*> squeue; 
  emutex m;
  eseqdb *seqdb;
//  ebasicarray<etax> *taxa;
  bool noveltaxa;
  bool print_align;
  bool print_hits;
} mtdata;

class epredtax
{
 public:
  earrayof<double,int> tax;
  efloatarray cutoff;
  
};

estr tax2str(emtdata& mtdata,long seqid)
{
  estr tmpstr;
  for (int t=0; t<mtdata.seqdb->taxa.size(); ++t){
    etax& tax(mtdata.seqdb->taxa.at(t));
    if (tax.seqs[seqid]==0x00)
      tmpstr+="-";
    else
      for (int l=0; l<tax.seqs[seqid]->tl.size(); ++l)
        tmpstr+=estr(";")+(tax.seqs[seqid]->tl[l].tid==-1?estr("NA"):tax.names[l].values(tax.seqs[seqid]->tl[l].tid));
    tmpstr+="\t";
  }
  tmpstr.del(0,1);
  return(tmpstr);
}

typedef estr (*outfmt_fd)(const etax&,const earrayof<double,int>&,const efloatarray&);

eoption<outfmt_fd> outfmt;

estr outfmt_simple(const etax& tax,const earrayof<double,int>& ptax,const efloatarray& mcfarr)
{
  estr res;
  for (int l=0; l<ptax.size() && l<tax.names.size(); ++l){
    if (mcfarr[l]<0.5) break;
    res+=(ptax.keys(l)==-1?estr("NA"):tax.names[l].values(ptax.keys(l)))+";";
  }
  res.del(-1);
  return(res);
}

estr outfmt_confidences(const etax& tax,const earrayof<double,int>& ptax,const efloatarray& mcfarr)
{
  estr res;
  for (int l=0; l<ptax.size() && l<tax.names.size(); ++l)
    res+=(ptax.keys(l)==-1?estr("NA"):tax.names[l].values(ptax.keys(l)))+"\t"+mcfarr[l]+"\t"+ptax.values(l)+"\t";
  res.del(-1);
  return(res);
}

void make_kmermask(eseqdb& db)
{
  int fcount=0;
  int okmercount=0;
  for (int i=0; i<db.otukmers.size(); ++i){
    ++okmercount;
    if (db.otukmers[i].size()>db.otus.size()*0.50){
      db.otukmermask[i]=1;
      db.otukmers[i].clear();
      ++fcount;
    }
  }
  cerr << "# filtered: " << fcount << endl;
}

template <class T>
void parallel_sortTask(T& arr,ebasicmtqueue<int>& q,int nthreads)
{
  int v;
  while (q.get(v)){
    std::sort(&arr[v*arr.size()/nthreads],&arr[(v+1)*arr.size()/nthreads]);
  }
}

template <class T>
void parallel_sortTaskMerge(T& arr,ebasicmtqueue<int>& q,int nthreads)
{
  int v;
  while (q.get(v)){
    std::inplace_merge(&arr[v*arr.size()/nthreads],&arr[(v+1)*arr.size()/nthreads],&arr[(v+2)*arr.size()/nthreads]);
  }
}


template <class T>
void parallel_sort(T& arr,int nthreads)
{
  nthreads-=nthreads%2; // make even
  if (nthreads<=1 || arr.size()<=19){
    std::sort(&arr[0],&arr[arr.size()]);
    return;
  }
  if (arr.size()/nthreads <= 10)
    nthreads=arr.size()/10;
  nthreads-=nthreads%2; // make even

  ebasicmtqueue<int> q;
  for (int i=0; i<nthreads; ++i)
    q.put(i);
  q.finish();

  ethreads t;
  t.run(parallel_sortTask<T>,evararray(arr,q),nthreads);
  t.wait();

  for (int tasks=nthreads; tasks>2; tasks>>=1){
    for (int i=0; i<tasks; i+=2)
      q.put(i);
    q.finish();
    t.run(parallel_sortTaskMerge<T>,evararray(arr,q),tasks);
    t.wait();
  }
  std::inplace_merge(&arr[0],&arr[arr.size()/2],&arr[arr.size()]);
}

void mtseqcluster(eseqdb& db,eseqdb& tmpdb,esearchws& searchws,estr& sid,eseq& s,estr& outstr)
{
  if (db.seqs.size()==0) ldie("db is empty");

  earray<epredinfo> pinfoarr;

  //TODO: improve count reset by using a bitmask and list for ids encountered to have a kmer from the query.
  //  everytime an element matching a kmer is found, set the bit in the bitmask and increment the pointer in the buffer
  // in the end iterate through the list and reset all bits and counts to zero. this list can be used to find the best matches
  seqsearch(sid,db,s,pinfoarr,searchws);
  if (pinfoarr.size()==0) {
    outstr+="!"+sid+"\t-\t-\n";
    return;
  }
  if (pinfoarr.size()>1){
    outstr+="!!"+sid+"\t-\t-\n";
    lwarn("clustering not made for long sequences: "+sid+" length: "+s.seqlen);
    return;
  }
  epredinfo &pinfo(pinfoarr[0]);
  einthashof<float> otuscore;
  float tweight=0.0;
  float topscore=pinfo.tophit.score();
  for (int l=pinfo.matchcounts.size()-1; l>=0; --l){
    int taxid=db.seqotu2[pinfo.matchcounts[l].seqid];
    if (!otuscore.exists(taxid)){
      tweight+=exp((1.0l-topscore/pinfo.matchcounts[l].score())*20.0l);
      otuscore.add(taxid,1.0);
    }
  }
  float cf=1.0/tweight;
  float idcf=(pinfo.tophit.identity()-0.96+0.02)/0.06+0.4;
  cf=(cf<idcf?cf:idcf);
  cf=(cf<0.0?0.0:(cf>1.0?1.0:cf));

  if (cf>=0.4){
    int bestotu=db.seqotu2[pinfo.tophit.seqid];
    int ibest=db.seqotu[pinfo.tophit.seqid];

//    db.seqs.add(sid,s);
//    db.seqotu.add(ibest);
//    db.otus[ibest].add(db.seqs.size()-1);
    tmpdb.seqotu2.add(bestotu);
//    db.otus2[bestotu].add(db.seqs.size()-1);
//        outstr+=pbuf->keys(i)+"\t"+db.seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.identity()+"\n";
    outstr+=sid+"\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+pinfo.tophit.identity()+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+(s.seqstart+pinfo.tophit.s1)+"\t"+(s.seqstart+pinfo.tophit.e1)+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t"+float(pinfo.tophit.e2-pinfo.tophit.s2)/mtdata.seqdb->seqs.values(pinfo.tophit.seqid).seqlen+"\n";
//        outstr+=sid+"\t"+db.seqs.keys(db.otus2[bestotu][0])+"\t"+pinfo.tophit.identity()+"\t"+(1.0/tweight)+"\t"+idcf+"\t\n";
  }else if (idcf<0.4) {
//    db.seqs.add(sid,s);
    db.seqotu2.add(db.otus2.size());
//    db.otus2.add(eintarray(db.seqs.size()-1));

    int ibest=db.seqotu[pinfo.tophit.seqid];
    if (pinfo.tophit.identity()*pinfo.tophit.alen1()/s.seqlen >= ccutoff){
      db.otus[ibest].add(db.seqs.size()-1);
      db.seqotu.add(ibest);
    } else {
      db.otus.add(eintarray(db.seqs.size()-1));
      ibest=db.otus.size()-1;
      db.seqotu.add(ibest);
      otukmeradd(db.otukmers,ibest,s,searchws.tmpkmers,ibest+1,akmers,0xFul,db.otukmermask);
    }
    outstr+=sid+"\t-\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+pinfo.tophit.identity()+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+(s.seqstart+pinfo.tophit.s1)+"\t"+(s.seqstart+pinfo.tophit.e1)+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t"+float(pinfo.tophit.e2-pinfo.tophit.s2)/mtdata.seqdb->seqs.values(pinfo.tophit.seqid).seqlen+"\n";
    if (mtdata.print_hits){
      outstr+="#\t";
      for (int l=pinfo.matchcounts.size()-1; l>0; --l){
        ealigndata& adata(pinfo.matchcounts[l]);
   //          outstr+="# "+mtdata.seqdb->seqs.keys(adata.seqid)+"\t"+adata.score()+"\t"+adata.identity()+"\t"+adata.matches+"\t"+adata.mismatches+"\t"+adata.gaps+"\t"+adata.s2+"\t"+adata.e2+"\t"+tax2str(mtdata,adata.seqid)+"\n";
        outstr+=mtdata.seqdb->seqs.keys(adata.seqid)+"\t"+adata.score()+"\t"+adata.identity()+"\t";
      }
      outstr+="\n";
    }
/*
    pinfo.tophit.profile.inv();
    outstr+=pinfo.tophit.profile.str() + "\n";
    outstr+=pinfo.tophit.align_str(s,mtdata.seqdb->seqs.values(pinfo.tophit.seqid));
    outstr+="\n";
*/
  }else {
    outstr+="#"+sid+"\t-\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+pinfo.tophit.identity()+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+(s.seqstart+pinfo.tophit.s1)+"\t"+(s.seqstart+pinfo.tophit.e1)+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t"+float(pinfo.tophit.e2-pinfo.tophit.s2)/mtdata.seqdb->seqs.values(pinfo.tophit.seqid).seqlen+"\n";
//        outstr+="#"+sid+"\t-\t"+pinfo.tophit.identity()+"\t"+(1.0/tweight)+"\t"+idcf+"\n";
  }
}

void taskSeqcluster2()
{
  estr outstr;
  esearchws searchws;
  searchws.init(MAXSIZE);
  eseqdb& db(*mtdata.seqdb);
  eseqdb tmpdb;

  // do normal clustering for a block of sequences
  while(1){
    estrarrayof<eseq> *pbuf=mtdata.squeue.getFull();
    if (pbuf==0x00) return;
//    sort(&(*pbuf).values(0),&(*pbuf).values(0)+pbuf->size());
    eintarray len_si(iheapsort(*pbuf));
    outstr.clear();

    for (int i=0; i<pbuf->size(); ++i){
      eseq& s(pbuf->values(len_si[i]));
      estr& sid(pbuf->keys(len_si[i]));
      mtseqcluster(db,tmpdb,searchws,sid,s,outstr);
      cout << outstr;
      outstr.clear();
//      eseq& s(pbuf->values(i));
//      estr& sid(pbuf->keys(i));
//      cout << s.seqlen << endl;
    }
    mtdata.squeue.returnEmpty(pbuf);

    mtdata.m.lock();
    cout << outstr;
    mtdata.m.unlock();
  }

  // for all new sequences do all agaisnt all

  // at the end add all remaining sequences to db and restart



//  delete searchws.bitmask;
}



void seqcluster(eseqdb& db,esearchws& searchws,estr& sid,eseq& s,estr& outstr)
{
  if (db.seqs.size()==0){
    db.seqs.add(sid,s);
    db.seqotu.add(db.otus.size());
    db.seqotu2.add(db.otus2.size());
    db.otus.add(eintarray(db.seqs.size()-1));
    db.otus2.add(eintarray(db.seqs.size()-1));
    otukmeradd(db.otukmers,db.otus.size()-1,s,searchws.tmpkmers,db.otus.size(),akmers,0xFul,db.otukmermask);
    outstr+=sid+"\t-\t-\n";
    return;
  }
  earray<epredinfo> pinfoarr;

  seqsearch(sid,db,s,pinfoarr,searchws);
  if (pinfoarr.size()==0) {
    outstr+="!"+sid+"\t-\t-\n";
    return;
  }
  if (pinfoarr.size()>1){
    outstr+="!!"+sid+"\t-\t-\n";
    lwarn("clustering not made for long sequences: "+sid+" length: "+s.seqlen);
    return;
  }
  epredinfo &pinfo(pinfoarr[0]);
  einthashof<float> otuscore;
  float tweight=0.0;
  float topscore=pinfo.tophit.score();
  for (int l=pinfo.matchcounts.size()-1; l>=0; --l){
    int taxid=db.seqotu2[pinfo.matchcounts[l].seqid];
    if (!otuscore.exists(taxid)){
      tweight+=exp((1.0l-topscore/pinfo.matchcounts[l].score())*20.0l);
      otuscore.add(taxid,1.0);
    }
  }
  float cf=1.0/tweight;
  float idcf=(pinfo.tophit.identity()-0.96+0.02)/0.06+0.4;
  cf=(cf<idcf?cf:idcf);
  cf=(cf<0.0?0.0:(cf>1.0?1.0:cf));

  if (cf>=0.4){
    int bestotu=db.seqotu2[pinfo.tophit.seqid];
    int ibest=db.seqotu[pinfo.tophit.seqid];

    db.seqs.add(sid,s);
    db.seqotu.add(ibest);
    db.otus[ibest].add(db.seqs.size()-1);
    db.seqotu2.add(bestotu);
    db.otus2[bestotu].add(db.seqs.size()-1);
//        outstr+=pbuf->keys(i)+"\t"+db.seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.identity()+"\n";
    outstr+=sid+"\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+pinfo.tophit.identity()+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+(s.seqstart+pinfo.tophit.s1)+"\t"+(s.seqstart+pinfo.tophit.e1)+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t"+float(pinfo.tophit.e2-pinfo.tophit.s2)/mtdata.seqdb->seqs.values(pinfo.tophit.seqid).seqlen+"\n";
//        outstr+=sid+"\t"+db.seqs.keys(db.otus2[bestotu][0])+"\t"+pinfo.tophit.identity()+"\t"+(1.0/tweight)+"\t"+idcf+"\t\n";
  }else if (idcf<0.4) {
    db.seqs.add(sid,s);
    db.seqotu2.add(db.otus2.size());
    db.otus2.add(eintarray(db.seqs.size()-1));

    int ibest=db.seqotu[pinfo.tophit.seqid];
    if (pinfo.tophit.identity()*pinfo.tophit.alen1()/s.seqlen >= ccutoff){
      db.otus[ibest].add(db.seqs.size()-1);
      db.seqotu.add(ibest);
    } else {
      db.otus.add(eintarray(db.seqs.size()-1));
      ibest=db.otus.size()-1;
      db.seqotu.add(ibest);
      otukmeradd(db.otukmers,ibest,s,searchws.tmpkmers,ibest+1,akmers,0xFul,db.otukmermask);
    }
    outstr+=sid+"\t-\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+pinfo.tophit.identity()+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+(s.seqstart+pinfo.tophit.s1)+"\t"+(s.seqstart+pinfo.tophit.e1)+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t"+float(pinfo.tophit.e2-pinfo.tophit.s2)/mtdata.seqdb->seqs.values(pinfo.tophit.seqid).seqlen+"\n";
    if (mtdata.print_hits){
      outstr+="#\t";
      for (int l=pinfo.matchcounts.size()-1; l>0; --l){
        ealigndata& adata(pinfo.matchcounts[l]);
   //          outstr+="# "+mtdata.seqdb->seqs.keys(adata.seqid)+"\t"+adata.score()+"\t"+adata.identity()+"\t"+adata.matches+"\t"+adata.mismatches+"\t"+adata.gaps+"\t"+adata.s2+"\t"+adata.e2+"\t"+tax2str(mtdata,adata.seqid)+"\n";
        outstr+=mtdata.seqdb->seqs.keys(adata.seqid)+"\t"+adata.score()+"\t"+adata.identity()+"\t";
      }
      outstr+="\n";
    }
/*
    pinfo.tophit.profile.inv();
    outstr+=pinfo.tophit.profile.str() + "\n";
    outstr+=pinfo.tophit.align_str(s,mtdata.seqdb->seqs.values(pinfo.tophit.seqid));
    outstr+="\n";
*/
  }else {
    outstr+="#"+sid+"\t-\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+pinfo.tophit.identity()+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+(s.seqstart+pinfo.tophit.s1)+"\t"+(s.seqstart+pinfo.tophit.e1)+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t"+float(pinfo.tophit.e2-pinfo.tophit.s2)/mtdata.seqdb->seqs.values(pinfo.tophit.seqid).seqlen+"\n";
//        outstr+="#"+sid+"\t-\t"+pinfo.tophit.identity()+"\t"+(1.0/tweight)+"\t"+idcf+"\n";
  }
}


void taskSeqcluster()
{
  estr outstr;
  esearchws searchws;
  searchws.init(MAXSIZE);
//  sws.idcount.reserve(1000000);
  searchws.idcount.init(2000000,0);
  searchws.idlist.init(2000000,0);
//  sws.idcount.init(db.otus.size()*2,0);
  eseqdb& db(*mtdata.seqdb);

  while(1){
    estrarrayof<eseq> *pbuf=mtdata.squeue.getFull();
    if (pbuf==0x00) return;
//    sort(&(*pbuf).values(0),&(*pbuf).values(0)+pbuf->size());
    eintarray len_si(iheapsort(*pbuf));
    outstr.clear();

    for (long i=0; i<pbuf->size(); ++i){
      eseq& s(pbuf->values(len_si[i]));
      estr& sid(pbuf->keys(len_si[i]));
      seqcluster(db,searchws,sid,s,outstr);
      cout << outstr;
      outstr.clear();
      if (i%100==0 && isatty(2))
//        fprintf(stderr,"\r# processing input... %li otus: %li otus2: %li",i,db.otus.size(),db.otus2.size());
        fprintf(stderr,"\r# processing input... %li o: %li o2: %li ti: %f ts: %f ti2: %f ts2: %f ta: %f",i,db.otus.size(),db.otus2.size(),ti,ts,ti2,ts2,ta);
//      eseq& s(pbuf->values(i));
//      estr& sid(pbuf->keys(i));
//      cout << s.seqlen << endl;
    }
    mtdata.squeue.returnEmpty(pbuf);

    mtdata.m.lock();
    cout << outstr;
    mtdata.m.unlock();
  }
//  delete searchws.bitmask;
}

void esearchws::init(int mxsize)
{
//  searchws.seqkmers.init(MAXSIZE,-1);
//  searchws.revseqkmers.init(MAXSIZE,-1);
//  searchws.kmerbitmask=new uint64_t[MAXSIZE/64+1];
//  searchws.bitmask=new uint64_t[(mtdata.seqdb->otus.size()*2)/64+1];
  otukmerpos.init(mtdata.seqdb->otus.size()*2,0);
  idcount.init(mtdata.seqdb->otus.size()*2,0);
//  idcount2.init(mtdata.seqdb->otus.size()*2,0);
  kmerpos.init(mxsize,-1);
  kmerpos2.init(mxsize,-1);
  kmerposrev.init(mxsize,-1);
  offset=0;
  offset2=0;
  kmermask.init(mxsize,-1);
  maskid=-1;
  tmpkmers.init(mxsize,-1);
}


void taskSeqsearch()
{
  estr outstr;
  esearchws searchws;
//  searchws.seqkmers.init(MAXSIZE,-1);
//  searchws.revseqkmers.init(MAXSIZE,-1);
//  searchws.kmerbitmask=new uint64_t[MAXSIZE/64+1];
//  searchws.bitmask=new uint64_t[(mtdata.seqdb->otus.size()*2)/64+1];
  searchws.otukmerpos.init(mtdata.seqdb->otus.size()*2,0);
  searchws.idcount.init(mtdata.seqdb->otus.size()*2,0);
//  searchws.idcount2.init(mtdata.seqdb->otus.size()*2,0);
  searchws.kmerpos.init(MAXSIZE,-1);
  searchws.kmerpos2.init(MAXSIZE,-1);
  searchws.kmerposrev.init(MAXSIZE,-1);
  searchws.offset=0;
  searchws.offset2=0;
  searchws.kmermask.init(MAXSIZE,-1);
  searchws.maskid=-1;
  while(1){
    estrarrayof<eseq> *pbuf=mtdata.squeue.getFull();
    if (pbuf==0x00) break;
    outstr.clear();

    for (int i=0; i<pbuf->size(); ++i){
      eseq& s(pbuf->values(i));
      earray<epredinfo> pinfoarr;
  //    ebasicarray<ealigndata> matchcounts;
      seqsearch(pbuf->keys(i),*mtdata.seqdb,s,pinfoarr,searchws);
      if (pinfoarr.size()==0) continue;
//      if (pinfo.tophit.seqid<0) continue;
  
      for (int pi=0; pi<pinfoarr.size(); ++pi){
        epredinfo& pinfo(pinfoarr[pi]);
        float taxcutoffmin=pinfo.matchcounts[0].identity();
        float bid=pinfo.tophit.identity();
        if (mtdata.seqdb->taxa.size() && mtdata.seqdb->taxa.at(0).seqs[pinfo.tophit.seqid]!=0x00){
          eseqtax &tmptaxhit(*mtdata.seqdb->taxa.at(0).seqs[pinfo.tophit.seqid]);
          // adjust id to closest gold hit
  //        if (tmptaxhit.bid>0.0) bid=tmptaxhit.bid*bid;
          if (tmptaxhit.bid>0.0 && bid>tmptaxhit.bid) bid=tmptaxhit.bid;
        }
      
//        outstr+=pbuf->keys(i)+"\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+bid+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t"+taxcutoffmin+"\t"+pinfo.matchcounts.size()+"\t";
  
//        outstr+=pbuf->keys(i)+"\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+bid+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+(s.seqstart+pinfo.tophit.s1)+"\t"+(s.seqstart+pinfo.tophit.e1)+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t"+taxcutoffmin+"\t"+pinfo.matchcounts.size()+"\t";
        outstr+=pbuf->keys(i)+"\t"+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+"\t"+pinfo.tophit.score()+"\t"+bid+"\t"+pinfo.tophit.matches+"\t"+pinfo.tophit.mismatches+"\t"+pinfo.tophit.gaps+"\t"+(s.seqstart+pinfo.tophit.s1)+"\t"+(s.seqstart+pinfo.tophit.e1)+"\t"+pinfo.tophit.s2+"\t"+pinfo.tophit.e2+"\t";
        for (int t=0; t<mtdata.seqdb->taxa.size(); ++t){
          etax& tax(mtdata.seqdb->taxa.at(t));
      
          earrayof<double,int> ptax;
          taxscore(ptax,pinfo,tax,searchws.taxcounts);
      
          bid=pinfo.tophit.identity();
          if (mtdata.seqdb->taxa.size() && mtdata.seqdb->taxa.at(t).seqs[pinfo.tophit.seqid]!=0x00){
            eseqtax &tmptaxhit(*mtdata.seqdb->taxa.at(t).seqs[pinfo.tophit.seqid]);
            // adjust id to closest gold hit
  //          if (tmptaxhit.bid>0.0 && bid>tmptaxhit.bid) bid=tmptaxhit.bid;
  //          if (tmptaxhit.bid>0.0) bid=0.5*(bid<tmptaxhit.bid?bid:tmptaxhit.bid)+0.5*bid*tmptaxhit.bid; 
            if (tmptaxhit.bid>0.0) {
              bid=bid*tmptaxhit.bid;
              for (int l=0; l<ptax.size(); ++l)
                ptax.values(l)=ptax.values(l)*tmptaxhit.tl[l].cf;
  //            if (ptax.values(l)>tmptaxhit.tl[l].cf) ptax.values(l)=tmptaxhit.tl[l].cf;
  //          ptax.values(l)=ptax.values(l)*tmptaxhit.tl[l].cf;
            }
          }
  
  /*
          efloatarray taxcutoff;
          if (mtdata.noveltaxa){
            taxcutoff.init(tax.names.size(),0.0);
            if (pinfo2.matchcounts.size()>0)
              findtaxcutoff(taxcutoff,pinfo2,ptax,tax);
          }
  */
      
  //        int hitlen=pinfo.tophit.e2-pinfo.tophit.s2;
     
          efloatarray mcfarr;
          mcfarr.init(ptax.size());
          float lastmcf=0.0;
          //  adjust computed confidences
          for (int l=MIN(ptax.size(),tax.names.size())-1; l>=0; --l){
            ldieif(ptax.keys(l)!=-1 && ptax.keys(l)>=tax.names[l].size(),"key out of tax: "+mtdata.seqdb->seqs.keys(pinfo.tophit.seqid)+" "+estr(l)+" "+ptax.keys(l)+" "+tax.names[l].size());
            mcfarr[l]=ptax.values(l);
  
            // if both fixed thresholds and precalculated identity thresholds exist, mix both in 3:1 ratio
            if (tax.cutoff.size()>0){ // if only fixed id threshold exists
              float ncf=(bid-tax.cutoff[l]+0.02)/tax.cutoffcoef[l];
              mcfarr[l]=ncf<ptax.values(l)?ncf:ptax.values(l);
            }
  //            mcf=ptax.values(l)+0.5*(bid-adjcutoff+0.02)/0.2;
  //            mcf=0.5*ptax.values(l)+0.5*(bid-adjcutoff+0.02)/0.05;
  //            mcf=0.5*ptax.values(l)+0.5*(bid-adjcutoff+0.02)/0.1;
  //            mcf=0.5*ptax.values(l)+0.5*(bid-adjcutoff+0.02)/0.2;
  
            if (mcfarr[l]>1.0) mcfarr[l]=1.0; else if (mcfarr[l]<0.0) mcfarr[l]=0.0;
            if (mcfarr[l]<lastmcf) mcfarr[l]=lastmcf;  // do not let confidences get smaller
            lastmcf=mcfarr[l];
          }

/*      
          for (int l=0; l<ptax.size() && l<tax.names.size(); ++l){
  //          outstr+=estr("\t")+(ptax.keys(l)==-1?estr("NA"):tax.names[l].values(ptax.keys(l)))+"\t"+ptax.values(l)+"\t"+(taxcutoff.size()?taxcutoff[l]:estr("-"));
            outstr+=estr("\t")+(ptax.keys(l)==-1?estr("NA"):tax.names[l].values(ptax.keys(l)))+"\t"+mcfarr[l]+"\t"+ptax.values(l);
          }
*/
          outstr+="\t";
          outstr+=outfmt.value()(tax,ptax,mcfarr);
          outstr+="\t";
        }
        outstr+="\n";
        mtdata.m.lock();
        cout << outstr; outstr.clear();
        mtdata.m.unlock();

        if (mtdata.print_hits){
          outstr+="#\t";
          etax& tax(mtdata.seqdb->taxa.at(0));
          for (int l=pinfo.matchcounts.size()-1; l>0; --l){
            ealigndata& adata(pinfo.matchcounts[l]);
  //          outstr+="# "+mtdata.seqdb->seqs.keys(adata.seqid)+"\t"+adata.score()+"\t"+adata.identity()+"\t"+adata.matches+"\t"+adata.mismatches+"\t"+adata.gaps+"\t"+adata.s2+"\t"+adata.e2+"\t"+tax2str(mtdata,adata.seqid)+"\n";
            outstr+=mtdata.seqdb->seqs.keys(adata.seqid)+"\t"+adata.score()+"\t"+adata.identity()+"\t";
          }
          outstr+="\n";
        }
  //      outstr+=pinfo.tophit.profile.str() + "\n";
        if (mtdata.print_align){
          pinfo.tophit.profile.inv();
          outstr+=pinfo.tophit.profile.str() + "\n";
          outstr+=pinfo.tophit.align_str(s,mtdata.seqdb->seqs.values(pinfo.tophit.seqid));
          outstr+="\n";
        }
//      if (pinfo.tophit.pair)
//        outstr+=estr("#pair: ")+pinfo.tophit.pair->score()+"\t"+pinfo.tophit.pair->identity()+"\t"+pinfo.tophit.pair->matches+"\t"+pinfo.tophit.pair->mismatches+"\t"+pinfo.tophit.pair->gaps+"\t"+pinfo.tophit.pair->s2+"\t"+pinfo.tophit.pair->e2+"\n";
//      outstr+="# "+pinfo.tophit.profile.str()+"\n";

/*
      for (int j=pinfo.matchcounts.size()-1; j>=0 && j>=pinfo.matchcounts.size()-10; --j){
        ealigndata &ad(pinfo.matchcounts[j]);
        outstr+="# "+mtdata.seqdb->seqs.keys(ad.seqid)+" "+ad.identity()+" "+ad.score()+" ";
        ldieif(mtdata.seqdb->taxa.size()==0,"no taxa");
        ldieif(ad.seqid>mtdata.seqdb->taxa[0].seqs.size(),"seqid larger than taxa: "+estr(ad.seqid)+" > "+mtdata.seqdb->taxa[0].seqs.size());
        if (mtdata.seqdb->taxa[0].seqs[ad.seqid]==0x00) outstr+="-";
        else{ 
          eseqtax& st(*mtdata.seqdb->taxa[0].seqs[ad.seqid]);
          for (int k=0; k<st.tl.size(); ++k)
            outstr+=";"+mtdata.seqdb->taxa[0].names[k][st.tl[k].tid];
        }
        outstr+="\n";
      }
*/
      }
    }

    mtdata.squeue.returnEmpty(pbuf);
    mtdata.m.lock();
    cout << outstr;
    mtdata.m.unlock();
  }
//  delete searchws.bitmask;
}

void make_cluster(const estr& cfile,eseqdb& db)
{
  eintarray kmerpos;
  ebasicarray<uint32_t> idcount;

  eintarray tmpkmers;
  tmpkmers.init(MAXSIZE,-1);
  fprintf(stderr,"# no clustering file found, performing clustering\n");

  // phase 1: find all cluster seeds, but do not add non-seeds in first phase
  eintarray len_si(iheapsort(db.seqs));
  db.otus.reserve(1000000l);
  db.otus.add(eintarray(len_si[0]));
  db.seqotu.init(db.seqs.size(),-1);


  tmpkmers.init(MAXSIZE,-1);
  otukmeradd(db.otukmers,0,db.seqs.values(db.otus[0][0]),tmpkmers,0,akmers,0xFul,db.otukmermask);
  idcount.reserve(1000000l);
  kmerpos.reserve(1000000l);
  idcount.init(db.otus.size(),0u);
  kmerpos.init(db.otus.size(),0u);

  int i;
  for (i=1; i<db.seqs.size(); ++i){
    eseq& s(db.seqs.values(len_si[i]));
    if (i%100==0) 
      fprintf(stderr,"\rseq: %i otus: %li",i,db.otus.size());
    idcount.init(db.otus.size(),0);
    kmerpos.init(db.otus.size(),0u);

    kmercount_single(db.otukmers,s,idcount,kmerpos);

    int ibest=0;
    for (int l=1; l<idcount.size(); ++l){
      if (idcount[l]>idcount[ibest]) ibest=l;
    }
    if (float(idcount[ibest])/s.seqlen >= 0.80) {
    } else {
      db.otus.add(eintarray(len_si[i]));
      ibest=db.otus.size()-1;
      db.seqotu[len_si[i]]=ibest;
      otukmeradd(db.otukmers,ibest,s,tmpkmers,ibest,akmers,0xFul,db.otukmermask);
    }
  }
  fprintf(stderr,"\rseq: %i clusters: %li\n",i,db.otus.size());

  fprintf(stderr,"phase 2:\n");
  // phase 2: after finding all seeds, add all non-seeds to the most similar clusters
  for (i=0; i<db.seqotu.size(); ++i){
    eseq& s(db.seqs.values(i));
    if (i%100==0)
      fprintf(stderr,"\rseq: %i",i);
    if (db.seqotu[i]!=-1) continue;

    idcount.init(db.otus.size(),0);
    kmerpos.init(db.otus.size(),0u);

    kmercount_single(db.otukmers,s,idcount,kmerpos);

    int ibest=0;

    for (int l=1; l<idcount.size(); ++l){
      if (idcount[l]>idcount[ibest]) ibest=l;
    }
    db.otus[ibest].add(i);
    db.seqotu[i]=ibest;
  }

  fprintf(stderr,"\rseq: %li clusters: %li\n",db.seqs.size(),(long)db.otus.size());
  if (cfile.len()){
    efile f(cfile,"w");
    for (int i=0; i<db.otus.size(); ++i){
      f.write(estr(i));
      for (int j=0; j<db.otus[i].size(); ++j)
        f.write(" "+estr(db.otus[i][j]));
      f.write("\n");
    }
    f.close();
  }

}

void load_cluster(const estr& cfile,eseqdb& db)
{
  eintarray tmpkmers;
  tmpkmers.init(MAXSIZE,-1);
  efile f;
  f.open(cfile,"r");
  estr line;
  tmpkmers.init(MAXSIZE,0);
  db.seqotu.init(db.seqs.size(),-1);
  while (!f.eof() && f.readln(line)){
    if (line.len()==0) continue;
    estrarray parts(line.explode(" "));

    int tmpi,repid=-1;
    for (int i=1; i<parts.size(); ++i){
      tmpi=parts[i].i();
//        if (ignseqs.size() && !ignseqs.exists(seqs.keys(tmpi)) || seqs.values(tmpi).tax!=0x00) {repid=tmpi; break;}
//      if (ignseqs.size() && ignseqs.exists(db.seqs.keys(tmpi))) continue;
      repid=tmpi;
      break;
    }
    if (repid==-1) continue; // no sequence found with taxonomic annotation

    db.otus.add(eintarray());
    eintarray &tmpo(db.otus[db.otus.size()-1]);
    tmpo.reserve(parts.size()-1);
    tmpi=parts[1].i();
//      otukmeradd(otukmers,otus.size()-1,seqs.values(tmpi),tmpkmers,tmpi,akmers,0xFFFFul);
//      if (db.otus.size()==0)
//        otukmeradd(db.otukmers,db.otus.size()-1,db.seqs.values(tmpi),tmpkmers,tmpi,akmers,0xFul);
//      else{
//      memset(bitmask,0x0u,(db.seqs.size()/64+1)*sizeof(uint64_t));
    uint32_t ibest=0u;
    short bcount=0;
    int kmerlen=0;
    eseq &s(db.seqs.values(tmpi));
//      kmercount_single_nopos(db.otukmers,s,bitmask,idcount,tmpkmers,db.otus.size(),kmerlen,ibest,bcount,akmers,0xFul);
//      if (idcount.size()==0){
    otukmeradd(db.otukmers,db.otus.size()-1,s,tmpkmers,db.otus.size(),akmers,0xFul,db.otukmermask);
//        tmpkmers.init(MAXSIZE,0);
//        kmerhash(db.otukmers,db.kmerlast,s,tmpkmers,ti++,db.seqs.size(),akmers,0xFul);
//      }else
//        otuaddkmerdiff(db.otukmers,db.seqs.values(db.otus[ibest][0]),s,tmpkmers,i,db.otus.size()-1,akmers,0xFul);
//      idcount.add(0);

    

    for (int i=1; i<parts.size(); ++i){
      tmpi=parts[i].i();
//        if (ignseqs.size() && !ignseqs.exists(seqs.keys(tmpi)) || seqs.values(tmpi).tax==0x00) continue; // do not add sequence if there is no taxonomic annotation
//      if (ignseqs.size() && ignseqs.exists(db.seqs.keys(tmpi))) continue; // do not add sequence if there is no taxonomic annotation
      tmpo.add(tmpi);
      db.seqotu[tmpi]=db.otus.size()-1;
    }
  }
  ldieif(db.otus.size()==0,"no clusters in database");
}

emtqueue<estrarray*> strqueue;
emutex m2;
estrarrayof<eseq> query;

void mtload_seqs_task()
{
  ebasicarray<eseq> tmpseqs;

  cout << "# starting" << endl;
  while (1) {
    estrarray *pbuf=strqueue.getFull();
    if (pbuf==0x00) break;
//    cout << "# reading" << endl;

    tmpseqs.init(pbuf->size());
    for (int i=0; i<pbuf->size(); ++i)
      tmpseqs[i].setseq(pbuf->values(i));
    m2.lock();
    for (int i=0; i<pbuf->size(); ++i)
      query.add(pbuf->keys(i),tmpseqs[i]);
    m2.unlock();
    strqueue.returnEmpty(pbuf);
  }
  cout << "# quitting" << endl;
}

void mtload_seqs(egzfile& f,estrarrayof<eseq>& seqs,int nthreads)
{
  ethreads t;
//  emutex m;
//  emtqueue<estrarray*> strqueue;
  estrarray *pbuf=0x00;
  int pbufi=0;

  estr str2id,str2seq,line;
  int maxseqlen=0;

  for (int i=0; i<nthreads*2; ++i){
    estrarray *sarr=new estrarray;
    for (int j=0; j<100; ++j) sarr->add(estr(),estr());
    strqueue.empty.push(sarr);
  }
  t.run(mtload_seqs_task,evararray(),nthreads);

  int count=0;
  f.readln(line);
  while (!f.eof()){

    if (pbuf==0x00){
      pbuf=strqueue.getEmpty();
      ldieif(pbuf==0x00,"pbuf empty");
      pbufi=0;
    }
    str2id=line;
    ldieif(str2id.len()==0 || str2id[0]!='>',"Unexpected line: "+str2id);

    str2seq.clear();
    while (f.readln(line) && line.len() && line[0]!='>') str2seq+=line;
    str2id.del(0,1);

    int i=str2id.findchr(" \t");
    if (i!=-1l) str2id.del(i); // only keep id up to first white space

    ++count;
    if (count%100000==0)
      printf("%8i\r",count);
    pbuf->keys(pbufi)=str2id;
    pbuf->values(pbufi)=str2seq;
    ++pbufi;
    if (pbufi>=pbuf->size()){
      strqueue.returnFull(pbuf);
      pbuf=0x00;
    }

//    db.seqind.add(str2id,db.seqs.size());
//    seqs.add(str2id,eseq(str2seq));
//    if (str2seq.len()>maxseqlen) maxseqlen=str2seq.len();
  }
  if (pbuf){
    while (pbuf->size()>pbufi)
      pbuf->erase(pbuf->size()-1);
    strqueue.returnFull(pbuf);
  }
  printf("%8i\n",count);
  f.close();
  strqueue.finish();
  t.wait();
}

void load_seqs(egzfile& f,estrarrayof<eseq>& seqs)
{
  estr str2id,str2seq,line;
  int maxseqlen=0;

  f.readln(line);
  while (!f.eof()){
    str2id=line;
    ldieif(str2id.len()==0 || str2id[0]!='>',"Unexpected line: "+str2id);

    str2seq.clear();
    while (f.readln(line) && line.len() && line[0]!='>') str2seq+=line;
    str2id.del(0,1);

    int i=str2id.findchr(" \t");
    if (i!=-1l) str2id.del(i); // only keep id up to first white space

//    db.seqind.add(str2id,db.seqs.size());
    seqs.add(str2id,eseq(str2seq));

    if (str2seq.len()>maxseqlen) maxseqlen=str2seq.len();
  }
  f.close();
}


void load_seqs(const estr& dbfile,eseqdb& db)
{
  estr str2id,str2seq,line;
  int maxseqlen=0;
  efile f;
  f.open(dbfile,"r");
  f.readln(line);
  while (!f.eof()){
 //f.readln(str2id)){
    str2id=line;
    ldieif(str2id.len()==0 || str2id[0]!='>',"Unexpected line: "+str2id);

    str2seq.clear();
    while (f.readln(line) && line.len() && line[0]!='>') str2seq+=line;
    str2id.del(0,1);

    int i=str2id.findchr(" \t");
    if (i!=-1l) str2id.del(i); // only keep id up to first white space

    db.seqind.add(str2id,db.seqs.size());
    db.seqs.add(str2id,eseq(str2seq));

    if (str2seq.len()>maxseqlen) maxseqlen=str2seq.len();
//    cout << ">" << str2id << endl;
//    cout << str2seq << endl;
//    cout << seqs.values(seqs.size()-1) << endl;
  }
  f.close();
}

void load_seqs_newotus(const estr& dbfile,eseqdb& db)
{
  load_seqs(dbfile,db);
//  db.seqotu.reserve(db.seqs.size());
//  eintarray tmpkmers;
//  tmpkmers.init(MAXSIZE,-1);
//  int tmpi=0;
  for (long i=db.seqotu2.size(); i<db.seqs.size(); ++i){
    db.seqotu2.add(db.otus.size());
    db.otus2.add(eintarray(i));
//    otukmeradd(db.otukmers,db.otus.size()-1,db.seqs.values(i),tmpkmers,tmpi,akmers,0xFul);
//    ++tmpi;
  }
}


void help()
{
  printf("MAPclust v%s\n",MAPCLUST_PACKAGE_VERSION);
  printf("by Joao F. Matias Rodrigues\n");
  printf("Institute of Molecular Life Sciences, University of Zurich, Switzerland\n");
  printf("\n");
  printf("Usage:\n");
  printf("    %s input.fa [<db.fa> [level]]\n",efile(getParser().args[0]).basename()._str);
  printf("\n");
  printf("Cluster a fasta file\n");
  printf("Example: mapclust rawreads.fa\n"); 
  printf("\n"); 

//  printf("After clustering:\n");
//  printf("%10s    %s\n","-makeotus <alignment> <mergelog> <threshold>","generate an OTU file at a given threshold");
//  printf("%10s    %s\n","-makeotus_mothur <alignment> <mergelog> <threshold>","generate a MOTHUR compatible OTU file at a given threshold");
//  printf("%10s    %s\n","-makereps <alignment> <otu>","generate a fasta file of OTU representatives. Sequences chosen have the minimum average distance to other sequences in the OTU.");
//  printf("\n");

  printf("Report bugs to: jfmrod@gmail.com\n");
  printf("http://meringlab.org/software/mapseq/\n");

  exit(0);
}

int emain()
{
  cout << "# mapclust v"<< MAPCLUST_PACKAGE_VERSION << " (" << __DATE__ << ")" << endl;
  initdlt();
  bool print_hits=false;
  bool print_align=false;
  epregister(print_hits);
  epregister(print_align);
  epregister(ccutoff);
  outfmt.choice=0;
  outfmt.add("simple",outfmt_simple);
  outfmt.add("confidences",outfmt_confidences);

  epregisterClassInheritance(eoption<outfmt_fd>,ebaseoption);
//  epregisterClassMethod4(eoption<outfmt_fd>,operator=,int,(const estr& val),"=");

  epregister(outfmt);
  epregister2(as.dropoff,"dropoff");

  epregisterFunc(help);
  
//  epregister(galign);
  bool benchmark=false;
//  epregister(benchmark);

  bool noveltaxa=false;
//  epregister(noveltaxa);

  int nthreads=4;
  epregister(nthreads);

  epregister(tophits);
  epregister(topotus);
  epregister(otulim);
  estr ffile;
  epregister(ffile);
  int step=0;

//  epregister(step);

//  estr ignfile;
//  epregister(ignfile);

//  epregister(kmer);
  eparseArgs();
  mtdata.print_align=print_align;
  mtdata.print_hits=print_hits;

  ldieif(getParser().args.size()<2,"syntax: mapseq <query> [db] [tax] [tax2] ...");
  estr cfile;

  for (unsigned int i=0; i<MAXSIZE; ++i)
    akmers[i]=0x00u;
  akmers[0x08ul]=1u; // AG
  akmers[0x0Cul]=1u; // AC
  akmers[0x09ul]=1u; // TG
  akmers[0x0Dul]=1u; // TC
  
  for (unsigned int i=0u; i<(1u<<4u); ++i){
    if (akmers[i]==1u) continue;
    for (unsigned int j=0u; j<(1u<<4u); ++j){
      if (akmers[j]!=1u) continue;
      akmers[(i<<4u)|j]=1u;
//      cout << kmer2str((i<<4u)|j) << endl;
    }
  }

//  cout << "4mer shift" << endl;

  for (unsigned int i=0u; i<(1u<<8u); ++i){
    if (akmers[i]==1u) continue;
    for (unsigned int j=0u; j<(1u<<8u); ++j){
      if (akmers[j]!=1u) continue;
      akmers[(i<<8u)|j]=1u;
//      cout << kmer2str((i<<8u)|j) << endl;
    }
  }


  
  estrhashof<int> filterhash;
  egzfile f;
  egzfile gzf;
  estr line;

  eseqdb db;

  if (ffile.len()){
    gzf.open(ffile,"r");
    while (!gzf.eof() && gzf.readln(line))
      filterhash.add(line,0);
  }

//  estrarrayof<eseq> seqs;
//  ebasicarray<etax> taxa;
//  ebasicarray<eintarray> taxcounts;


/*
  estrhash ignseqs;
  if (ignfile.len()>0){
    f.open(ignfile,"r");
    while (!f.eof()){
      f.readln(line);
      ignseqs.add(line,estr());
    }
  }
*/

//  estrarrayof<eseq> query;
/*
  while(f.readln(line)){
    if (line.len()==0 || line[0]=='#') continue;
    args=line.explode(" ");
    ldieif(args.size()<2,"id and sequence not found in line: "+line);
    seqs.add(args[0],args[1]);
  }

  f.open(argv[2],"r");
  while(f.readln(line)){
    if (line.len()==0 || line[0]=='#') continue;
    args=line.explode(" ");
    ldieif(args.size()<2,"id and sequence not found in line: "+line);
    query.add(args[0],args[1]);
  }
*/

/*
  estr dbfile=estr(DATAPATH)+"/mapseqref.fasta";
  if (getParser().args.size()>2)
    dbfile=getParser().args[2];
  else if (!efile(dbfile).exists())
    dbfile=dirname(getSystem().getExecutablePath())+"/share/mapseq/mapseqref.fasta";
  if (!efile(dbfile).exists()) ldie("fasta db not found: "+dbfile);

  load_seqs(dbfile,db);

  cerr << "# loaded " << db.seqs.size() << " sequences" << endl;
  ldieif(db.seqs.size()==0,"empty database");
*/

  etimer t1;
  t1.reset();

//  earray<eintarray> otus;
//  ebasicarray<eintarray*> otukmers;

  db.otukmermask.init(MAXSIZE,0);
  db.otukmers.init(MAXSIZE);

//  eintarray tmpkmers;
//  eintarray kmercounts;
  int i;

  if (getParser().args.size()>2) {
    cerr << "# loading sequence db: " << getParser().args[2] << endl;
    load_seqs(getParser().args[2],db);
    if (efile(getParser().args[2]+".mscluster").exists())
      load_cluster(getParser().args[2]+".mscluster",db);
    else
      make_cluster(getParser().args[2]+".mscluster",db);

    ldieif(getParser().args.size()<=3,"need taxa level to use for mapping");
    cerr << "# loading taxonomy file: " << getParser().args[3] << endl;
    load_taxa(getParser().args[3],db);
    ldieif(getParser().args.size()<=4,"need taxa level to use for mapping");
    make_taxcluster(db.taxa[0],getParser().args[4],db);
    if (getParser().args.size()>5)
      load_seqs_newotus(getParser().args[2],db);
  }
  
  cerr << "# taxonomies: " << db.taxa.size() << endl;
  for (int i=0; i<db.taxa.size(); ++i){
    cerr << "# tax levels: " << db.taxa[i].names.size() << endl;
//    taxcounts.init(taxa[i].names.size());
    for (int j=0; j<db.taxa[i].names.size(); ++j){
      cerr << "# tax: " << i << " level: " << j << " ("<< db.taxa[i].names[j].size() << ")" << endl;
//      taxcounts[i][j].init(taxa[i].names[j].size(),-1);
    }
  }



//  eintarray seqotu;
//  db.seqotu.init(db.seqs.size(),-1);

//  tmpkmers.init(MAXSIZE,-1);



  const int MAXSEQS=1000000;
  uint64_t bitmask[MAXSEQS/64+1];

/*
  cfile=dbfile+".mscluster";
  if (efile(cfile).exists()){ // load clustering
    load_cluster(cfile,db);
  }else{ // perform clustering
    make_cluster(cfile,db);
  }

  int fcount=0;
  int okmercount=0;
  for (i=0; i<db.otukmers.size(); ++i){
//    if (db.otukmers[i]!=0x00){
      ++okmercount;
      if (db.otukmers[i].size()>db.otus.size()*0.50){
        db.otukmers[i].clear();
        ++fcount;
      }
//    }
  }
*/

//  cerr << "# fcount: " << fcount << " otukmercount: " << okmercount << endl;



//  earray<estrarray> tax;
//  earray<estrhashof<int> > taxind;

/*
  int alen;
  int match[11];
  int j,k;
*/

  t1.reset();

/*
  eintarray seqkmers;
  eintarray revseqkmers,revtmpkmers;
  eintarray tmpkmers2;
*/

  ethreads t;
  emutex m;
 
//  if (estr(getParser().args[1])=="-")
//    f.open(stdin);
//  else
    f.open(getParser().args[1],"r");

  mtdata.seqdb=&db;
  mtdata.noveltaxa=noveltaxa;


  t1.reset();
//  taskSeqsearchRead(f);
//  cerr << "# done processing " << " seqs (" << t1.lap()*0.001 << "s)" << endl;
  
  nthreads=1;
  t.run(taskSeqcluster,evararray(),nthreads);

  for (int i=0; i<nthreads*2; ++i){
    estrarrayof<eseq> *sarr=new estrarrayof<eseq>;
    for (int j=0; j<100; ++j) sarr->add(estr(),eseq());
    mtdata.squeue.empty.push(sarr);
  }

  cout << "#query\tdbhit\tbitscore\tidentity\tmatches\tmismatches\tgaps\tquery_start\tquery_end\tdbhit_start\tdbhit_end\t";
  if (outfmt.key()=="simple"){
    for (int i=0; i<db.taxa.size(); ++i)
      cout << "\t" << db.taxa[i].name;
  }else if (outfmt.key()=="confidences"){
    for (int i=0; i<db.taxa.size(); ++i){
      etax& tax(db.taxa[i]);
      if (tax.levels.size()==0){
        if (tax.names.size())
          cout << "\t"<<tax.name<<":taxlevel0\tcombined_cf\tscore_cf";
        for (int j=1; j<tax.names.size(); ++j)
          cout << "\ttaxlevel" << j << "\tcombined_cf\tscore_cf";
      }else{
        cout << "\t"<< tax.name <<":" << tax.levels[0] << "\tcombined_cf\tscore_cf";
        for (int j=1; j<tax.levels.size(); ++j)
          cout << "\t" << tax.levels[j] << "\tcombined_cf\tscore_cf";
      }
    }

  }
  cout << endl;

//  mtload_seqs(f,query,10);
  load_seqs(f,query);
  mtdata.squeue.returnFull(&query);
  mtdata.squeue.finish();
 
  fprintf(stderr,"\r# processing input... %li otus: %li\n",query.size(),db.otus2.size());
  t.wait();
  cout << "# otus: " << db.otus2.size() << endl;
  cerr << "# done processing " << query.size() << " seqs (" << t1.lap()*0.001 << "s)" << endl;
  return(0);


//  seqcluster(query,db);
  

  estr str2id,str2seq;
//  cerr << "# processing input... ";
  long seqcount=0l;
  fprintf(stderr,"# processing input... ");
  fflush(stderr);
  lerrorif(!f.readln(line),"unable to read query file");
  estrarrayof<eseq> *cbuf=0x00;
  int cbufind=0;
  const int MAXSEQLEN=100000;
  long seqstart;
  str2seq.reserve(MAXSEQLEN);
  while (!f.eof()){
    if (cbuf==0x00){
      cbuf=mtdata.squeue.getEmpty();
      cbufind=0;
    }

    seqstart=0l;
    str2id=line;
    ldieif(str2id.len()==0 || str2id[0]!='>',"Unexpected line: "+str2id);
    str2id.del(0,1);
    int i=str2id.findchr(" \t");
    if (i!=-1l) str2id.del(i); // only keep id up to first white space

    str2seq.clear();

    if (filterhash.size() && !filterhash.exists(str2id)){
      while (!f.eof() && f.readln(line) && line.len() && line[0]!='>');
      continue;
    }

    //TODO: read a limited number of nucleotides each time, so the code does not depend on line breaks. For every nucleotide "chunk", compress the nucleotides and add to sequence
    while (!f.eof() && f.readln(line) && line.len() && line[0]!='>'){
      if (str2seq.len()+line.len()>MAXSEQLEN){
        cbuf->keys(cbufind)=str2id;
        eseq& s(cbuf->values(cbufind));
        s.setseq(str2seq);
        s.seqstart=seqstart;
        ++cbufind;
        if (cbufind==cbuf->size()){
          mtdata.squeue.returnFull(cbuf);
          cbuf=mtdata.squeue.getEmpty();
          cbufind=0;
        }
        seqstart+=s.seqlen;
        str2seq.clear();
      }
      str2seq+=line;
    }

    cbuf->keys(cbufind)=str2id;
    eseq& s(cbuf->values(cbufind));
    s.setseq(str2seq);
    s.seqstart=seqstart;
    ++cbufind;
    if (cbufind==cbuf->size()){
      mtdata.squeue.returnFull(cbuf);
      cbuf=0x00;
    }
    ++seqcount;
    if (seqcount%100==0 && isatty(2))
      fprintf(stderr,"\r# processing input... %li otus: %li otus2: %li",seqcount,db.otus.size(),db.otus2.size());
//      fprintf(stderr,"\r# processing input... %li  ti: %f ts: %f ti2: %f ts2: %f ta: %f tdp1: %f tdp2: %f tdpfl: %f tdpmd: %f",seqcount,ti,ts,ti2,ts2,ta,tdp1,tdp2,tdpfl,tdpmd);
 }
  if (cbuf!=0x00){
    for (int i=cbuf->size()-1; i>=cbufind; --i) cbuf->erase(i);
//    mtdata.seqs.add(cbuf);
    mtdata.squeue.returnFull(cbuf);
  }
  mtdata.squeue.finish();
  
  fprintf(stderr,"\r# processing input... %li otus: %li\n",seqcount,db.otus2.size());
  f.close();
  t.wait();
  cout << "# otus: " << db.otus2.size() << endl;
  cerr << "# done processing " << seqcount << " seqs (" << t1.lap()*0.001 << "s)" << endl;

  return(0);
}