SHORE API
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
suffix.hpp
Go to the documentation of this file.
1 
2 /*
3  * Copyright 2008,2009,2010,2011,2012 Stephan Ossowski, Korbinian Schneeberger,
4  * Felix Ott, Joerg Hagmann, Alf Scotland, Sebastian Bender
5  *
6  * This file is part of SHORE.
7  *
8  * SHORE is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * SHORE is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with SHORE. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
25 
26 #ifndef SHORE_ALGO_SUFFIX_HPP__
27 #define SHORE_ALGO_SUFFIX_HPP__
28 
29 #include <algorithm>
30 #include <iterator>
31 #include <set>
32 #include <string>
33 #include <utility>
34 
35 #include "shore/container/intpack.hpp"
36 #include "shore/base/stringops.hpp"
37 
38 namespace shore {
39 
44 template<typename Iter,typename Oter>
45 void count_elements(Iter ibeg,Iter iend,Oter obeg,Oter oend,
46  bool initialize_counts=true,int shift=0)
47 {
48  if(initialize_counts)
49  std::fill(obeg,oend,0);
50 
51  for(Iter i=ibeg;i!=iend;++i)
52  {
53  ++(obeg[(*i)+shift]);
54  }
55 }
56 
57 template<typename Iter>
58 void counts_to_intervals(Iter beg,Iter end)
59 {
60  typedef typename std::iterator_traits<Iter>::value_type value_type;
61 
62  if(beg==end)
63  return;
64 
65  Iter last=beg;
66  Iter curr=beg;
67  ++curr;
68 
69  while(curr!=end)
70  {
71  (*curr)+=(*last);
72  last=curr;
73  ++curr;
74  }
75 }
76 
83 template<typename Iter,typename Oter>
84 void suffix_array(Iter ibeg,Iter iend,Oter obeg,bool initialize_sa=true,
85  typename std::iterator_traits<Iter>::value_type imax=
86  std::numeric_limits<typename std::iterator_traits<Iter>::value_type>::max())
87 {
88  //std::cerr<<"."<<std::endl;
89 
90  typedef typename std::iterator_traits<Iter>::value_type ivalue_type;
91  typedef typename std::iterator_traits<Oter>::value_type ovalue_type;
92 
93  typedef typename std::reverse_iterator<Oter> ROter;
94 
95  const size_t n=iend-ibeg;
96  const Oter oend=obeg+n;
97 
98  if(initialize_sa)
99  std::fill(obeg,oend,0);
100 
101  if(n<2)
102  return;
103 
104  const ROter r_obeg(oend);
105  const ROter r_oend(obeg);
106 
107  // store L/S types
108  intpack types(n);
109  types.repack(1);
110 
111  ivalue_type prev=*(iend-1);
112  ivalue_type s_type=ivalue_type(0);
113  ovalue_type nlms=0;
114 
115  if(prev<0)
116  throw std::runtime_error("suffix_array: negative value "
117  +shore::to_string(prev));
118 
119  for(size_t i=n-2;;--i)
120  {
121  const ivalue_type curr=ibeg[i];
122 
123  if(curr<0)
124  throw std::runtime_error("suffix_array: negative value "
125  +shore::to_string(curr));
126 
127  if(curr<(prev+s_type))
128  {
129  types[i]=s_type=ivalue_type(1);
130  }
131  else if(s_type!=0)
132  {
133  s_type=ivalue_type(0);
134  ++nlms;
135  }
136  prev=curr;
137 
138  if(i==0)
139  break;
140  }
141 
142  // create bucket ptrs
143  // todo?: buckets may be stored in array space for most recursions
144  intpack buckets(imax+2);
145  buckets.repack_for(n);
146  count_elements(ibeg,iend,buckets.begin(),buckets.end(),false);
147  counts_to_intervals(buckets.begin(),buckets.end());
148 
149  if(nlms>0)
150  {
151  const Oter lms_end=obeg+nlms;
152 
153  // put lmses into buckets
154  s_type=ivalue_type(intpack::int_type(types[0]));
155  for(size_t i=1;i<n-1;++i)
156  {
157  const ivalue_type tmp=ivalue_type(intpack::int_type(types[i]));
158  if(tmp>s_type)
159  {
160  obeg[--buckets[ibeg[i]]]=i;
161  }
162  s_type=tmp;
163  }
164 
165  // induced sort
166  count_elements(ibeg,iend,buckets.begin(),buckets.end(),true,1);
167  counts_to_intervals(buckets.begin(),buckets.end());
168 
169  Iter ilast=iend-1;
170  obeg[buckets[*ilast]]=n-1;
171  ++buckets[*ilast];
172 
173  for(Oter i=obeg;i!=oend;++i)
174  {
175  const ovalue_type x=*i;
176  if((x!=0)&&(types[x-1]==0))
177  {
178  const ivalue_type c=ibeg[x-1];
179  obeg[buckets[c]]=x-1;
180  ++buckets[c];
181  }
182  }
183 
184  count_elements(ibeg,iend,buckets.begin(),buckets.end(),true);
185  counts_to_intervals(buckets.begin(),buckets.end());
186 
187  for(ROter i=r_obeg;i!=r_oend;++i)
188  {
189  const ovalue_type x=*i;
190  if((x!=0)&&(types[x-1]==1))
191  {
192  const ivalue_type c=ibeg[x-1];
193  --buckets[c];
194  obeg[buckets[c]]=x-1;
195  }
196  }
197 
198  // move lms indeces to start
199  for(Oter i=obeg,j=obeg;i!=oend;++i)
200  {
201  const ovalue_type idx=*i;
202  if((idx>0)&&(types[idx]>types[idx-1]))
203  {
204  (*j)=(*i);
205  ++j;
206 
207  if(j==lms_end)
208  break;
209  }
210  }
211 
212  // check lms for equality
213  Oter nam_beg=obeg+(n>>1);
214 
215  std::fill(nam_beg,oend,n-1);
216 
217  ovalue_type pprev=*obeg;
218  ovalue_type nnam=0;
219  nam_beg[pprev>>1]=nnam;
220  for(Oter i=(obeg+1);i!=lms_end;++i)
221  {
222  const ovalue_type p=*i;
223  bool eq=true;
224  const size_t smax=n-std::max(p,pprev);
225  for(size_t j=0;j<smax;++j)
226  {
227  if((ibeg[p+j]!=ibeg[pprev+j])||(types[p+j]!=types[pprev+j]))
228  {
229  eq=false;
230  break;
231  }
232  if((j>0)&&(types[p+j]>types[p+j-1]))
233  break;
234  }
235  if(!eq)
236  ++nnam;
237  nam_beg[p>>1]=nnam;
238  pprev=p;
239  }
240 
241  // move lms names to back
242  for(Oter i=oend-1,j=oend-1;i!=nam_beg-1;--i)
243  {
244  const ovalue_type v=*i;
245  if(v!=(n-1))
246  {
247  *j=v;
248  --j;
249  }
250  }
251  nam_beg=oend-nlms;
252 
253  // recurse
254  if(nnam<nlms)
255  suffix_array(nam_beg,oend,obeg,true,nnam);
256  else for(size_t i=0;i<nlms;++i)
257  obeg[nam_beg[i]]=i;
258 
259  // induce final result
260  count_elements(ibeg,iend,buckets.begin(),buckets.end(),true);
261  counts_to_intervals(buckets.begin(),buckets.end());
262 
263  // replace names by lms ptrs
264  {
265  s_type=ivalue_type(intpack::int_type(types[0]));
266  Oter j=nam_beg;
267  for(size_t i=1;i<n-1;++i)
268  {
269  const ivalue_type tmp=ivalue_type(intpack::int_type(types[i]));
270  if(tmp>s_type)
271  {
272  (*j)=i;
273  ++j;
274  }
275  s_type=tmp;
276  }
277  }
278 
279  for(Oter i=obeg;i!=lms_end;++i)
280  {
281  (*i)=nam_beg[*i];
282  }
283 
284  std::fill(lms_end,oend,0);
285 
286  {
287  ovalue_type j;
288  for(ROter i(lms_end);i!=r_oend;++i)
289  {
290  j=(*i);
291  (*i)=0;
292  obeg[--buckets[ibeg[j]]]=j;
293  }
294  }
295  }
296 
297  // induced sort
298  count_elements(ibeg,iend,buckets.begin(),buckets.end(),true,1);
299  counts_to_intervals(buckets.begin(),buckets.end());
300 
301  Iter ilast=iend-1;
302  obeg[buckets[*ilast]]=n-1;
303  ++buckets[*ilast];
304 
305  for(Oter i=obeg;i!=oend;++i)
306  {
307  const ovalue_type x=*i;
308  if((x!=0)&&(types[x-1]==0))
309  {
310  const ivalue_type c=ibeg[x-1];
311  obeg[buckets[c]]=x-1;
312  ++buckets[c];
313  }
314  }
315 
316  count_elements(ibeg,iend,buckets.begin(),buckets.end(),true);
317  counts_to_intervals(buckets.begin(),buckets.end());
318 
319  for(ROter i=r_obeg;i!=r_oend;++i)
320  {
321  const ovalue_type x=*i;
322  if((x!=0)&&(types[x-1]==1))
323  {
324  const ivalue_type c=ibeg[x-1];
325  --buckets[c];
326  obeg[buckets[c]]=x-1;
327  }
328  }
329 }
330 
331 template<typename Ater,typename Oter>
332 void suffix_ranks(Ater abeg,Ater aend,Oter rbeg)
333 {
334  typedef typename std::iterator_traits<Oter>::value_type value_type;
335  value_type r=0;
336  for(Ater i=abeg;i!=aend;++i)
337  rbeg[*i]=r++;
338 }
339 
341 template<typename Ster,typename Ater>
342 std::pair<size_t,size_t> suffix_lcpmax(Ster sbeg,Ster send,Ater abeg)
343 {
344  std::pair<size_t,size_t> ret(0,0);
345 
346  const size_t size=send-sbeg;
347  const Ater aend=abeg+size;
348 
349  intpack rnk(size);
350  rnk.repack_for(size-1);
351  suffix_ranks(abeg,aend,rnk.begin());
352 
353  size_t h=0;
354 
355  for(intpack::iterator i=rnk.begin();i!=rnk.end();++i)
356  {
357  const intpack::int_type r=*i;
358 
359  if(r>0)
360  {
361  Ster k=sbeg+(*(abeg+r-1))+h;
362  Ster j=sbeg+(*(abeg+r))+h;
363  while((j!=send)&&(k!=send)&&((*j)==(*k)))
364  {
365  ++h;
366  ++j;
367  ++k;
368  }
369  }
370  if(h>ret.first)
371  {
372  ret.first=h;
373  ret.second=r;
374  }
375  if(h>0)
376  --h;
377  }
378 
379  return ret;
380 }
381 
383 template<typename Ster,typename Ater,typename Oter>
384 void suffix_lcp(Ster sbeg,Ster send,Ater abeg,Oter obeg)
385 {
386  typedef typename std::iterator_traits<Oter>::value_type value_type;
387 
388  const size_t size=send-sbeg;
389  const Ater aend=abeg+size;
390 
391  intpack rnk(size);
392  rnk.repack_for(size-1);
393  suffix_ranks(abeg,aend,rnk.begin());
394 
395  value_type h=0;
396 
397  for(intpack::iterator i=rnk.begin();i!=rnk.end();++i)
398  {
399  const intpack::int_type r=*i;
400  if(r>0)
401  {
402  Ster k=sbeg+(*(abeg+r-1))+h;
403  Ster j=sbeg+(*(abeg+r))+h;
404  while((j!=send)&&(k!=send)&&((*j)==(*k)))
405  {
406  ++h;
407  ++j;
408  ++k;
409  }
410  }
411  obeg[r]=h;
412  if(h>0)
413  --h;
414  }
415 }
416 
418 template<typename Ster,typename Ater,typename Qter>
419 typename std::pair<Ater,Ater> suffix_query(Ster sbeg,Ster send,
420  Ater abeg,Ater aend,
421  Qter qbeg,Qter qend)
422 {
423  typedef typename std::pair<Ater,Ater> range;
424  typedef typename std::iterator_traits<Ster>::value_type svalue_t;
425  typedef typename std::iterator_traits<Ater>::value_type avalue_t;
426  typedef typename std::iterator_traits<Qter>::value_type qvalue_t;
427 
428  if(abeg==aend)
429  return range(abeg,aend);
430  if(qbeg==qend)
431  return range(abeg,abeg);
432 
433  range ret(abeg,aend);
434 
435  const size_t k=send-sbeg;
436 
437  size_t ofs=0;
438  size_t idxL=0;
439  size_t idxR=0;
440 
441  for(Qter q=qbeg;q!=qend;++q)
442  {
443  idxL=(*ret.first)+ofs;
444  idxR=(*(ret.second-1))+ofs;
445 
446  Ater bound=ret.second;
447  size_t step=(bound-ret.first)>>1;
448  while((idxL>=k)||(sbeg[idxL]<static_cast<svalue_t>(*q)))
449  {
450  const Ater itmp=ret.first+step;
451  const size_t idxM=(*itmp)+ofs;
452  if((idxM>=k)||(sbeg[idxM]<static_cast<svalue_t>(*q)))
453  {
454  ret.first=itmp;
455  ++ret.first;
456  }
457  else
458  {
459  bound=itmp;
460  ++ret.first;
461  }
462  step=(bound-ret.first)>>1;
463 
464  if(ret.first==ret.second)
465  return ret;
466 
467  idxL=(*ret.first)+ofs;
468  }
469 
470  bound=ret.first;
471  step=(ret.second-bound)>>1;
472 
473  while((idxR<k)&&(sbeg[idxR]>static_cast<svalue_t>(*q)))
474  {
475  const Ater itmp=ret.second-step;
476  const size_t idxM=(*(itmp-1))+ofs;
477  if((idxM<k)&&(sbeg[idxM]>static_cast<svalue_t>(*q)))
478  {
479  ret.second=itmp;
480  --ret.second;
481  }
482  else
483  {
484  bound=itmp;
485  --ret.second;
486  }
487  step=(ret.second-bound)>>1;
488 
489  if(ret.first==ret.second)
490  return ret;
491 
492  idxR=(*(ret.second-1))+ofs;
493  }
494 
495  ++ofs;
496  }
497 
498  return ret;
499 }
500 
501 template<typename Ater,typename Qter>
502 struct suffix_query_result
503 {
505  Ater array_beg;
507  Ater array_end;
509  Ater array_ins;
511  Qter query_end;
512 
513  suffix_query_result(Ater b,Ater e,Qter q)
514  :array_beg(b),array_end(e),array_ins(b),query_end(q)
515  {}
516 };
517 
519 template<typename Ster,typename Ater,typename Qter>
520 suffix_query_result<Ater,Qter> suffix_query_prefix(Ster sbeg,Ster send,
521  Ater abeg,Ater aend,
522  Qter qbeg,Qter qend)
523 {
524  typedef typename std::pair<Ater,Ater> range;
525  typedef typename std::iterator_traits<Ster>::value_type svalue_t;
526  typedef typename std::iterator_traits<Ater>::value_type avalue_t;
527  typedef typename std::iterator_traits<Qter>::value_type qvalue_t;
528 
529  suffix_query_result<Ater,Qter> ret(abeg,aend,qbeg);
530 
531  if(abeg==aend)
532  return ret;
533  if(qbeg==qend)
534  return ret;
535 
536  const size_t k=send-sbeg;
537 
538  size_t idxL=0;
539  size_t idxR=0;
540 
541  Ater cur_array_end=aend;
542 
543  for(size_t ofs=0;ret.query_end!=qend;++ret.query_end,++ofs)
544  {
545  idxL=(*ret.array_ins)+ofs;
546  idxR=(*(cur_array_end-1))+ofs;
547 
548  Ater bound=cur_array_end;
549  size_t step=(bound-ret.array_ins)>>1;
550  while((idxL>=k)||(sbeg[idxL]<static_cast<svalue_t>(*ret.query_end)))
551  {
552  const Ater itmp=ret.array_ins+step;
553  const size_t idxM=(*itmp)+ofs;
554  if((idxM>=k)||(sbeg[idxM]<static_cast<svalue_t>(*ret.query_end)))
555  {
556  ret.array_ins=itmp;
557  ++ret.array_ins;
558  }
559  else
560  {
561  bound=itmp;
562  ++ret.array_ins;
563  }
564  step=(bound-ret.array_ins)>>1;
565 
566  if(ret.array_ins==cur_array_end)
567  return ret;
568 
569  idxL=(*ret.array_ins)+ofs;
570  }
571 
572  bound=ret.array_ins;
573  step=(cur_array_end-bound)>>1;
574 
575  while((idxR<k)&&(sbeg[idxR]>static_cast<svalue_t>(*ret.query_end)))
576  {
577  const Ater itmp=cur_array_end-step;
578  const size_t idxM=(*(itmp-1))+ofs;
579  if((idxM<k)&&(sbeg[idxM]>static_cast<svalue_t>(*ret.query_end)))
580  {
581  cur_array_end=itmp;
582  --cur_array_end;
583  }
584  else
585  {
586  bound=itmp;
587  --cur_array_end;
588  }
589  step=(cur_array_end-bound)>>1;
590 
591  if(ret.array_ins==cur_array_end)
592  return ret;
593 
594  idxR=(*(cur_array_end-1))+ofs;
595  }
596 
597  ret.array_beg=ret.array_ins;
598  ret.array_end=cur_array_end;
599  }
600 
601  return ret;
602 }
603 
611 template<typename Ster,typename Ater,typename Query>
612 typename std::pair<Ater,Ater> suffix_query_incremental(Ster sbeg,Ster send,
613  Ater abeg,Ater aend,
614  const Query q,
615  const size_t ofs)
616 {
617  typedef typename std::pair<Ater,Ater> range;
618  typedef typename std::iterator_traits<Ster>::value_type svalue_t;
619  typedef typename std::iterator_traits<Ater>::value_type avalue_t;
620  typedef Query qvalue_t;
621 
622  if(abeg==aend)
623  return range(abeg,aend);
624 
625  range ret(abeg,aend);
626 
627  const size_t k=send-sbeg;
628 
629  size_t idxL=0;
630  size_t idxR=0;
631 
632  idxL=(*ret.first)+ofs;
633  idxR=(*(ret.second-1))+ofs;
634 
635  Ater bound=ret.second;
636  size_t step=(bound-ret.first)>>1;
637  while((idxL>=k)||(sbeg[idxL]<static_cast<svalue_t>(q)))
638  {
639  const Ater itmp=ret.first+step;
640  const size_t idxM=(*itmp)+ofs;
641  if((idxM>=k)||(sbeg[idxM]<static_cast<svalue_t>(q)))
642  {
643  ret.first=itmp;
644  ++ret.first;
645  }
646  else
647  {
648  bound=itmp;
649  ++ret.first;
650  }
651  step=(bound-ret.first)>>1;
652 
653  if(ret.first==ret.second)
654  return ret;
655 
656  idxL=(*ret.first)+ofs;
657  }
658 
659  bound=ret.first;
660  step=(ret.second-bound)>>1;
661 
662  while((idxR<k)&&(sbeg[idxR]>static_cast<svalue_t>(q)))
663  {
664  const Ater itmp=ret.second-step;
665  const size_t idxM=(*(itmp-1))+ofs;
666  if((idxM<k)&&(sbeg[idxM]>static_cast<svalue_t>(q)))
667  {
668  ret.second=itmp;
669  --ret.second;
670  }
671  else
672  {
673  bound=itmp;
674  --ret.second;
675  }
676  step=(ret.second-bound)>>1;
677 
678  if(ret.first==ret.second)
679  return ret;
680 
681  idxR=(*(ret.second-1))+ofs;
682  }
683 
684  return ret;
685 }
686 
688 template<typename Ster,typename Ater,typename Qter>
689 typename std::pair<Ater,Ater> suffix_query(
690  Ster sbeg,Ster send,Ater abeg,Qter qbeg,Qter qend)
691 {
692  return suffix_query(sbeg,send,abeg,abeg+(send-sbeg),qbeg,qend);
693 }
694 
696 template<typename Ster,typename Ater,typename Qter>
697 suffix_query_result<Ater,Qter> suffix_query_prefix(
698  Ster sbeg,Ster send,Ater abeg,Qter qbeg,Qter qend)
699 {
700  return suffix_query_prefix(sbeg,send,abeg,abeg+(send-sbeg),qbeg,qend);
701 }
702 
703 template<typename Ster,typename Ater>
704 std::ostream& suffix_dump(std::ostream& os,Ster sbeg,Ster send,Ater abeg)
705 {
706  const Ater aend=abeg+(send-sbeg);
707  for(Ater i=abeg;i!=aend;++i)
708  {
709  os<<(*i)<<'\t';
710  for(Ster j=sbeg+(*i);j!=send;++j)
711  os<<(*j);
712  os<<'\n';
713  }
714  return os;
715 }
716 
718 template<typename Ster,typename Ater>
719 bool suffix_test(Ster sbeg,Ster send,Ater abeg)
720 {
721  bool pass=true;
722  const size_t size=send-sbeg;
723  const Ater aend=abeg+size;
724 
725  if(size>0)
726  {
727  {
728  std::cerr<<"sort"<<std::endl;
729  std::vector<typename std::iterator_traits<Ater>::value_type> v(abeg,aend);
730  std::sort(v.begin(),v.end());
731 
732  std::cerr<<"test uniqueness"<<std::endl;
733  if(std::unique(v.begin(),v.end())!=v.end())
734  {
735  std::cerr<<"sarr uniqueness FAIL"<<std::endl;
736  pass=false;
737  }
738 
739  std::cerr<<"test min"<<std::endl;
740  if((*std::min_element(v.begin(),v.end()))!=0)
741  {
742  std::cerr<<"sarr min FAIL"<<std::endl;
743  pass=false;
744  }
745 
746  std::cerr<<"test max"<<std::endl;
747  if((*std::max_element(v.begin(),v.end()))!=(size-1))
748  {
749  std::cerr<<"sarr max FAIL"<<std::endl;
750  pass=false;
751  }
752  }
753 
754  if(pass&&(abeg!=aend))
755  {
756  std::cerr<<"test ordering"<<std::endl;
757  Ater last=abeg;
758  for(Ater i=(abeg+1);i!=aend;++i)
759  {
760  Ster j=sbeg+std::max(*i,*last);
761  Ster k=sbeg+std::min(*i,*last);
762  if(std::lexicographical_compare(j,send,k,send)!=((*last)>(*i)))
763  {
764  std::cerr<<"sarr sort FAIL"<<std::endl;
765 
766  //VPRN(size);
767  //VPRN(i-abeg);
768  //VPRN(*i);
769  //VPRN(*last);
770  //VPRN(seq.substr(*i));
771  //VPRN(seq.substr(*last));
772 
773  pass=false;
774  }
775  last=i;
776  }
777  }
778  }
779 
780  if(pass)
781  std::cerr<<"sarr OK"<<std::endl;
782 
783  return pass;
784 }
785 
786 } // namespace
787 
788 #endif // SHORE_ALGO_SUFFIX_HPP__
789