# -*- coding: iso-8859-1 -*- # Copyright 20003 - 2008: Julien Bourdaillet (julien.bourdaillet@lip6.fr), Jean-Gabriel Ganascia (jean-gabriel.ganascia@lip6.fr) # This file is part of MEDITE. # # MEDITE is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # MEDITE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with Foobar; if not, write to the Free Software # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA from __future__ import generators #from __future__ import absolute_import import sys, logging, weakref, bisect, time, os, os.path import numpy.oldnumeric as Numeric import psyco import recouvrement import _suffix_tree class SuffixTreeInterface(object): """Fonctions de parcours de l'arbre et de recherche des chaines répétées""" def generate_inner_nodes_2(self): 'Iterator through only MAximal exact matches' for n in self.post_order_nodes_2: if not n.is_leaf: yield n def generate_MEM(self): 'Iterator through only Maximal exact matches' def dfs(n): nb_children_leaf = 0 for c in n.children: if c.is_leaf: nb_children_leaf += 1 for m in dfs(c): yield m if nb_children_leaf>=1: yield n for n in dfs(self.root): yield n def generate_post_order_nodes(self): 'Iterator through all nodes in the tree in post-order.' def dfs(n): for c in n.children: for m in dfs(c): yield m yield n for n in dfs(self.root): yield n def generate_post_order_nodes_max(self): 'Iterator through all nodes in the tree in post-order.' def dfs(n): for c in n.children: for m in dfs(c): if m.nb_occ>1: yield m if n.nb_occ>1: yield n for n in dfs(self.root): if n.nb_occ>1: yield n def generate_pre_order_nodes(self): 'Iterator through all nodes in the tree in post-order.' def dfs(n): yield n for c in n.children: for m in dfs(c): yield m for n in dfs(self.root): yield n def generate_leaves(self): 'Iterator through all leaves in the tree.' for n in self.post_order_nodes: if n.is_leaf: yield n def generate_inner_nodes(self): 'Iterator through all leaves in the tree.' for n in self.post_order_nodes: if not n.is_leaf: yield n # set class properties MEM = property(generate_MEM, None, None, "post_order_nodes2") post_order_nodes = property(generate_post_order_nodes, None, None, "post_order_nodes") post_order_nodes_max = property(generate_post_order_nodes_max, None, None, "post_order_nodes_max") pre_order_nodes = property(generate_pre_order_nodes, None, None, "pre_order_nodes") leaves = property(generate_leaves, None, None, "leaves") inner_nodes = property(generate_inner_nodes, None, None, "inner_nodes") inner_nodes_2 = property(generate_inner_nodes_2, None, None, "inner_nodes_2") def _annotate_nodes(self): """Ajoute à chaque noeud interne(répété) ses positions dans la séquence. """ for n in self.post_order_nodes: if n.is_leaf: n.path_indices = n.path_position else: n.path_indices = [] for c in n.children: try: n.path_indices.extend(c.path_indices) except TypeError: n.path_indices.append(c.path_indices) #print path_indices n.nb_occ = len(n.path_indices) def get_MEM(self, min_size=1): """Renvoie un dico de tous les Maximal Exact Matches de taille min min_size Le dico est indexé par la taille de MEM qui renvoie une liste de toutes les positions de cette taille. Linéaire(nb de MEM) < linéaire(taille de la séquence)""" # à chaque position de seq_repeat doit correspondre la position # de la fin de la répétition la plus longue commençant à cette position #self._annotate_nodes() seq_repeat = {} for n in self.generate_inner_nodes():#generate_MEM(): longueur_MEM = n.edge_label_end - n.path_position if longueur_MEM >= min_size: n.path_indices = [] for c in n.children: if c.is_leaf: n.path_indices.append(c.path_position) elif len(c.path_indices) == 1: n.path_indices.extend(c.path_indices) for pi in n.path_indices: seq_repeat[pi] = pi+longueur_MEM e="""for i in xrange(pi, pi+longueur_chaine): if not seq_repeat.has_key(i): seq_repeat[i] = pi+longueur_chaine elif pi+longueur_chaine < seq_repeat[i]: break else: seq_repeat[i] = pi+longueur_chaine""" #print seq_repeat dic_MEM = {} for debut in seq_repeat: fin = seq_repeat[debut] if (fin > debut): try: dic_MEM[fin-debut].append(debut) except KeyError: dic_MEM[fin-debut] = [debut] #print dic_MEM return dic_MEM class SuffixTree(SuffixTreeInterface,_suffix_tree.SuffixTree): """A higher-level wrapper around the C suffix tree type, _suffix_tree.SuffixTree. This class adds a few methods to the suffix tree, methods that are more easily expressed in Python than in C, and that can be written using the primitives exported from C. """ def __init__(self,s): """Build a suffix tree from the input string s. The string must not contain the special symbol $.""" if '$' in s: raise "The suffix tree string must not contain $!" self.sequence = s #print len(s) #s = unicode(s,'utf-8') #raw_unicode_escape') #print len(s) #if isinstance(s, unicode): # s = s.encode('utf-8') #print "SuffixTree.__init__ type(s) ",type(s) ," len ",len(s) _suffix_tree.SuffixTree.__init__(self,s) class SuffixNodeU(dict): def __init__(self, debut, fin, slink=None, depth=0): self.suffixLink = slink self.depth = depth # représente le début et la fin de l'arc qui pointe sur ce noeud # ces 2 valeurs servent pour le parcours de l'arbre ultérieur à # sa construction permetant la recherche des MEMs self.edge_label_begin = debut # debut de l'arc self.edge_label_end = fin # fin de l'arc # path_position correspond au 1er caractère de la branche, # celui du 1er caractère du noeud fils de root # ainsi texte[path_position:edge_label_end] décrit la chaine depuis root def changeValues(self, debut, fin): self.edge_label_begin = debut self.edge_label_end = fin def getChildren(self): '''Iterator through all the node's children''' for k,p,node in self.itervalues(): yield node children = property(getChildren, None, None, "children of the node") def isLeaf(self): if len(self) == 0: return True else: return False is_leaf = property(isLeaf, None, None, "leaf or not leaf ?") def __str__(self, seq=None, space='', short=False, charOnly=None): res = "deb="+str(self.edge_label_begin)+" / fin="+str(self.edge_label_end) try: res += " / path_position="+str(self.path_position) except AttributeError: pass if seq is not None: res += " / " +seq[self.edge_label_begin:self.edge_label_end] res += "\n" if short and seq is not None: res = seq[self.edge_label_begin:self.edge_label_end]+'\n' space += '|---' if charOnly is not None and self.has_key(charOnly): k,p,node = self[charOnly] res += space + node.__str__(seq, space, short) else: for char,(k,p,node) in self.iteritems(): res += space + node.__str__(seq, space, short) return res class SuffixTree_(SuffixTreeInterface): """Construit un suffix tree avec l'algorithme d'Ukkonen""" def __str__(self): return self.root.__str__(seq=self.sequence) def __init__(self, sequence, generalizedST=False): assert isinstance(sequence,basestring) if not generalizedST: if '¤' in sequence: raise "The suffix tree string must not contain ¤!" self.sequence = sequence+'¤' self.inf = len(sequence) # infini self.bottom = SuffixNodeU(-2,-2) self.root = SuffixNodeU(0,0,self.bottom) s = self.root k = 0 ; i = -1 for t in self.sequence: # pour chaque char t en position i, consturction de ST(sequence[:i]) i += 1 # transition t entre bottom en root self.bottom[t] = -2,-2,self.root # ajout de t en renvoie du end point de ST(sequence[:i]) s,k = self.update(s,k,i) # le end point n'est pas forcément explicite, on le canonize s,k = self.canonize(s,k,i) # end point + 1 char = active point de ST(sequence[:i+1]) #print self self.annotate_path_position() self.sequence = sequence def update(self, s,k,i): """Ajout du caractère en i """ # active point = s,k,i-1 oldroot = self.root t = self.sequence[i] # active point = end point ? renvoie l'active point explicité end_point,r = self.testAndSplit(s,k,i-1,t) while not end_point: r2 = SuffixNodeU(i,self.inf) r[t] = i,self.inf,r2 # crée nouvelle transition # crée lien suffixe du noeud de l'étape précédente if not (oldroot is self.root): oldroot.suffixLink = r oldroot = r s,k = self.canonize(s.suffixLink,k,i-1) # remonte lien suffixe end_point,r = self.testAndSplit(s,k,i-1,t) # end point ? if not (oldroot is self.root): oldroot.suffixLink = r return s,k # renvoie le end point def testAndSplit(self, s,k,p,t): """Test si s,k,p est le end point et sinon explicite et renvoie s,k,p""" if k <= p: # s,k,p implicite tk = self.sequence[k] k2,p2,s2 = s[tk] # le char suivant l'arc s[tk] = t ? => end point if t == self.sequence[k2+p-k+1]: return True,s else:# sinon on split l'arc tk2 = self.sequence[k2] end = k2+p-k r = SuffixNodeU(k2,end+1)#+1 pour avoir un slicing correct en python s[tk2] = k2, end, r s2.changeValues(end+1, p2+1)#MAJ valeur pour le slicing python tk3 = self.sequence[end+1] r[tk3] = end+1, p2, s2 return False,r # s,k,p explicite elif not s.has_key(t): return False,s # pas de tk-transition else: return True,s # tk-transition => end point def canonize(self, s,k,p): """Si s,k,p implicite, renvoie son plus proche ancêtre explicite""" if p < k: return s,k # s,k,p explicite else: tk = self.sequence[k] k2,p2,s2 = s[tk] # si l'arc qui commence en s[tk] est inférieur à w, on avance dessus while p2-k2 <= p-k: k = k+p2-k2+1 s = s2 if k <= p: # essaye un nouvel arc tk = self.sequence[k] k2,p2,s2 = s[tk] return s,k def annotate_path_position(self): """ajoute l'attribut path_position qui correspond à la position du 1er caractère de la branche""" node = self.root self.root.path_position = 0 self.__annotate_path_position(node, 0) def __annotate_path_position(self, node, longueur): for s,k,fils in node.itervalues(): if len(fils) == 0: # feuille # début de l'arc - longueur des arcs traversés depuis root fils.path_position = fils.edge_label_begin - longueur else: fils.path_position = fils.edge_label_begin - longueur self.__annotate_path_position(fils, longueur+fils.edge_label_end-fils.edge_label_begin) class TrueGeneralisedSuffixTree(SuffixTree): """A suffix tree for a set of strings.""" def __init__(self,sequences): '''Build a generalised suffix tree. The strings must not contain the special symbols $ or ascii numbers from 1 to the number of sequences.''' self.sequences = sequences self.start_positions = [0] self.concat_string = '' for i in xrange(len(sequences)): if chr(i+1) in sequences[i]: raise "The suffix tree string must not contain chr(%d)!"%(i+1) self.concat_string += sequences[i]+chr(i+1) self.start_positions += [len(self.concat_string)] self.start_positions += [self.start_positions[-1]+1] # empty string self.sequences += [''] #print self.start_positions #print "GeneralisedSuffixTree.__init__ type(self.concat_string) ",type(self.concat_string) ," len ",len(self.concat_string) SuffixTree.__init__(self,self.concat_string) self._annotate_nodes() def _translate_index(self,idx): 'Translate a concat-string index into a (string_no,idx) pair.' for i in xrange(len(self.start_positions)-1): if self.start_positions[i] <= idx < self.start_positions[i+1]: return (i,idx-self.start_positions[i]) raise IndexError, "Index out of range: "+str(idx) def _annotate_nodes(self): for n in self.post_order_nodes: if n.is_leaf: seq,idx = self._translate_index(n.path_position) n.path_indices = [(seq, idx)] n.sequences = [seq] n.nb_occ = len(n.path_indices) else: path_indices = [] ; sequences = [] for c in n.children: path_indices += c.path_indices sequences += c.sequences seqs_in_subtree = {} for s in sequences: seqs_in_subtree[s] = 1 n.path_indices = path_indices n.sequences = [s for s in seqs_in_subtree] n.nb_occ = len(n.path_indices) def shared_substrings(self,minimum_length=0): '''Iterator through shared sub-strings. Returned as a list of triples (sequence,from,to).''' for n in self.inner_nodes: if len(n.sequences) == len(self.sequences) - 1: l = n.edge_label_end - n.path_position if l >= minimum_length: yield [(seq, idx, idx+l) for (seq,idx) in n.path_indices] # for (seq,idx) in n.path_indices: # print "idx="+str(idx)+" / idx+1="+str(idx+1) # yield (seq, idx, idx+l) def inter_inclus(self,i,d): if (i[0]<=d[0]<=i[1] and i[0]<=d[1]<=i[1]): return True else: return False def shared_substrings2naif(self,long_min_pivots): '''Iterator through shared sub-strings. Returned as a list of triples (sequence,from,to).''' DRes={} maxIdxDres=0 longTexte1=len(self.sequences[0])+1 for n in self.inner_nodes: #for n in self.post_order_nodes: #print "n.pathlabel= *"+n.path_label+"* edge_label= *"+n.edge_label+"*" # print "edge_label= *"+n.edge_label+"*" long_path=n.edge_label_end-n.path_position if long_path1: if DRes.has_key(long_path): DRes[long_path].extend(ladd) else: maxIdxDres=max(maxIdxDres,long_path) DRes[long_path]=ladd LRes=[] for i in xrange(maxIdxDres,0,-1): if DRes.has_key(i): print str(i)+" "+str(DRes[i]) for x in DRes[i]: for y in LRes: if self.inter_inclus(y,x): break else: LRes.append(x) print "text n="+str(len(self.sequences[0])+len(self.sequences[1]))+" / leaves="+str(nleaf)+" / internal nodes="+str(nnode) return LRes def shared_substrings2(self,long_min_pivots): '''Iterator through shared sub-strings. Returned as a list of triples (sequence,from,to).''' lInter=[] longTexte1=len(self.sequences[0])+1 nleaf=nnode=0 #for n in self.inner_nodes: for n in self.post_order_nodes: if n.is_leaf: nleaf+=1 else: nnode+=1 #print "n.pathlabel= *"+n.path_label+"* edge_label= *"+n.edge_label+"*" # print "edge_label= *"+n.edge_label+"*" long_path=n.edge_label_end-n.path_position if long_path1: lInter=MediteAppli.union2(ladd,lInter) # for x in ladd: # lInter = MediteAppli.ajout_intervalle(lInter,x) sys.stderr.write("text n="+str(len(self.sequences[0])+len(self.sequences[1]))+" / leaves="+str(nleaf)+" / internal nodes="+str(nnode)+"\n") return lInter def shared_substrings3(self,long_min_pivots): '''Iterator through shared sub-strings. Returned as a list of triples (sequence,from,to).''' longTexte1=len(self.sequences[0]) # longueur du 1er texte t=self.sequences[0] + self.sequences[1] # texte total longTexte=len(t) # longueur du texte total lMax=[] for i in xrange(longTexte): lMax.append(0) nleaf=nnode=0 for n in self.post_order_nodes: if n.is_leaf: nleaf+=1 else: nnode+=1 long_path=n.edge_label_end-n.path_position if long_path=fin):continue else: ladd.append((deb,fin)) if (len(ladd)>1 and ladd[0][0]=longTexte1): # if (ladd[0][0]lMax[i]: lMax[i]=x[1] #else: break dic={} prevInserted=lastmax=0 nbPos=0 for i in xrange(longTexte): nbPos+=1 pos=lMax[i] long=pos-i if pos!=prevInserted and pos>0 and pos>lastmax: texte=t[i:pos] if not dic.has_key(texte): dic[texte]=[] dic[texte].append(i) prevInserted=pos if pos>lastmax: lastmax=pos #sys.stderr.write("len dic = "+str(len(dic))+"\n") #print dic dic2={} for x in dic: a=dic[x] if (len(x)=longTexte1)): # sys.stderr.write("SUPP "+x+" "+str(a)+"\n") continue else:dic2[x]=a # sys.stderr.write("len dic2 = "+str(len(dic2))+"\n") # sys.stderr.write("text n="+str(len(self.sequences[0])+len(self.sequences[1]))+ # " / leaves="+str(nleaf)+" / internal nodes="+str(nnode)+ # " / nbPos="+str(nbPos)+"\n") #print dic2 return dic2,lMax def get_repeat(self, min_size=1): """Renvoie un dictionnaire de toutes les chaines répétées. Indexé par la taille des chaines. En valeur une liste de positions. min_size: taille mini des chaines répétées Linéaire(nb de noeuds de l'arbre)=linéaire(taille de la séquence)""" longueur_seq1 = len(self.sequences[0]) dic = {} for n in self.generate_inner_nodes():# chaines répétées longueur_chaine = n.edge_label_end - n.path_position if longueur_chaine >= min_size: seq1 = seq2 = False lPos = [] for (seq,pos) in n.path_indices: if seq == 0: seq1 = True lPos.append(pos) else: seq2 = True lPos.append(pos + longueur_seq1) if seq1 and seq2: try: dic[longueur_chaine].extend(lPos) #print n.path_position except KeyError: dic[longueur_chaine] = lPos return dic def get_seq_repeat_list(self,min_size=1): # à chaque position de seq_repeat doit correspondre la position # de la fin de la répétition la plus longue commençant à cette position #print self.root.__str__(self.concat_string, short=True) longueur_seq1 = len(self.sequences[0]) seq_repeat = [] for i in xrange(len(self.concat_string)): seq_repeat.append((i,i)) i=nbpitot=0 for n in self.generate_inner_nodes():#generate_MEM(): i+=1 #print i,n.__str__(self.concat_string) #if n.path_positionlongueur_seq1: # n.edge_label_end= longueur_seq1 #assert (0<=n.path_position<=n.edge_label_begin<=n.edge_label_end<=longueur_seq1 or # longueur_seq1<=n.path_position<=n.edge_label_begin<=n.edge_label_end<=len(self.sequences[0])+len(self.sequences[1])) longueur_chaine = n.edge_label_end - n.path_position n.path_indices = [] for c in n.children: if c.is_leaf: n.path_indices.append(c.path_position) #elif len(c.path_indices) == 1: n.path_indices.extend(c.path_indices) else: n.path_indices.extend(c.path_indices) if longueur_chaine >= min_size and len(n.path_indices) > 1: #print n.path_indices seq1 = seq2 = False for pos in n.path_indices: #print self.concat_string[pos:pos+longueur_chaine] if pos < longueur_seq1: seq1 = True else: seq2 = True if seq1 and seq2: for pi in n.path_indices: #try: # if seq_repeat[pi][1] < pi+longueur_chaine: # seq_repeat[pi] = (pi,pi+longueur_chaine) #except KeyError: # seq_repeat[pi] = (pi,pi+longueur_chaine) for pos in xrange(pi,pi+longueur_chaine): #for pos in xrange(pi+longueur_chaine,pi,-1): #if (seq_repeat[pos][1] < pi+longueur_chaine or # pi < seq_repeat[pos][0] and seq_repeat[pos][1] == pi+longueur_chaine): # seq_repeat[pos] = (pi,pi+longueur_chaine) # si bloc décrit à la pos sourante est inclus dans le bloc qu'on ajoute if (pi <= seq_repeat[pos][0] <= seq_repeat[pos][1] <= pi+longueur_chaine or seq_repeat[pos][0] <= pi <= seq_repeat[pos][1] <= pi+longueur_chaine): seq_repeat[pos] = (pi,pi+longueur_chaine) #for pos in xrange(pi+longueur_chaine-1,pi-1,-1): # if (pi <= seq_repeat[pos][0] <= seq_repeat[pos][1] <= pi+longueur_chaine or #pi <= seq_repeat[pos][0] <= pi+longueur_chaine <= seq_repeat[pos][1]): # seq_repeat[pos][0] <= pi <= seq_repeat[pos][1] < pi+longueur_chaine): # seq_repeat[pos] = (pi,pi+longueur_chaine) #print seq_repeat #if seq_repeat[pi][1] < pi+longueur_chaine: # seq_repeat[pi] = (pi,pi+longueur_chaine) #nbpitot+=1 #seq_repeat[pi] = pi+longueur_chaine #else: print n.path_indices #print seq_repeat return seq_repeat def get_seq_repeat(self,min_size=1): # à chaque position de seq_repeat doit correspondre la position # de la fin de la répétition la plus longue commençant à cette position #print self.root.__str__(self.concat_string, short=True) longueur_seq1 = len(self.sequences[0]) seq_repeat_deb = Numeric.zeros(len(self.concat_string),Numeric.Int) seq_repeat_fin = Numeric.zeros(len(self.concat_string),Numeric.Int) for i in xrange(len(self.concat_string)): seq_repeat_deb[i] = i seq_repeat_fin[i] = i i=nbpitot=0 for n in self.generate_inner_nodes():#generate_MEM(): i+=1 #print i,n.__str__(self.concat_string) #if n.path_positionlongueur_seq1: # n.edge_label_end= longueur_seq1 #assert (0<=n.path_position<=n.edge_label_begin<=n.edge_label_end<=longueur_seq1 or # longueur_seq1<=n.path_position<=n.edge_label_begin<=n.edge_label_end<=len(self.sequences[0])+len(self.sequences[1])) longueur_chaine = n.edge_label_end - n.path_position n.path_indices = [] all_children_leaf = True for c in n.children: if c.is_leaf: n.path_indices.append(c.path_position) #elif len(c.path_indices) == 1: n.path_indices.extend(c.path_indices) else: n.path_indices.extend(c.path_indices) all_children_leaf = False listePos = [i-1 for i in n.path_indices] listeCarac = [] for i in listePos: assert -1 <= i < len(self.concat_string),(i,len(self.concat_string)) if i == -1: x = chr(4) else: x = self.concat_string[i] listeCarac.append(x) #for idx,pos in listePos: # if idx == 0: listeCarac.append(self.sequences[0][pos]) # else: listeCarac.append(self.sequences[1][pos]) setCarac = set(listeCarac) if (longueur_chaine >= min_size and len(n.path_indices) > 1 and all_children_leaf and len(setCarac)==len(listeCarac)): #print n.path_indices seq1 = seq2 = False for pos in n.path_indices: #print self.concat_string[pos:pos+longueur_chaine] if pos < longueur_seq1: seq1 = True else: seq2 = True if seq1 and seq2: for pi in n.path_indices: #try: # if seq_repeat[pi][1] < pi+longueur_chaine: # seq_repeat[pi] = (pi,pi+longueur_chaine) #except KeyError: # seq_repeat[pi] = (pi,pi+longueur_chaine) for pos in xrange(pi,pi+longueur_chaine): #for pos in xrange(pi+longueur_chaine,pi,-1): #if (seq_repeat[pos][1] < pi+longueur_chaine or # pi < seq_repeat[pos][0] and seq_repeat[pos][1] == pi+longueur_chaine): # seq_repeat[pos] = (pi,pi+longueur_chaine) # si bloc décrit à la pos sourante est inclus dans le bloc qu'on ajoute if (pi <= seq_repeat_deb[pos] <= seq_repeat_fin[pos] <= pi+longueur_chaine or seq_repeat_deb[pos] <= pi <= seq_repeat_fin[pos] <= pi+longueur_chaine): seq_repeat_deb[pos] = pi seq_repeat_fin[pos] = pi+longueur_chaine #for pos in xrange(pi+longueur_chaine-1,pi-1,-1): # if (pi <= seq_repeat[pos][0] <= seq_repeat[pos][1] <= pi+longueur_chaine or #pi <= seq_repeat[pos][0] <= pi+longueur_chaine <= seq_repeat[pos][1]): # seq_repeat[pos][0] <= pi <= seq_repeat[pos][1] < pi+longueur_chaine): # seq_repeat[pos] = (pi,pi+longueur_chaine) #print seq_repeat #if seq_repeat[pi][1] < pi+longueur_chaine: # seq_repeat[pi] = (pi,pi+longueur_chaine) #nbpitot+=1 #seq_repeat[pi] = pi+longueur_chaine #else: print n.path_indices #print seq_repeat return seq_repeat_deb,seq_repeat_fin def get_seq_repeat2(self,min_size=1): # à chaque position de seq_repeat doit correspondre la position # de la fin de la répétition la plus longue commençant à cette position #print self.root.__str__(self.concat_string, short=True) texte = self.sequences[0] + '$' + self.sequences[1] longueur_seq1 = len(self.sequences[0]) seq_repeat = [] for i in xrange(len(self.concat_string)): seq_repeat.append((i,i)) i=nbpitot=0 for n in self.generate_inner_nodes():#generate_MEM(): i+=1 longueur_chaine = n.edge_label_end - n.path_position n.path_indices = [] all_children_leaf = True seq1 = seq2 = False car_seq1 = [] ; car_seq2 = [] for c in n.children: if c.is_leaf: n.path_indices.append(c.path_position) dic_car_seq1 = {} ; dic_car_seq2 = {} if c.path_position <= longueur_seq1: seq1 = True try: dic_car_seq1[texte[c.path_position]].append(c.path_position) except KeyError: dic_car_seq1[texte[c.path_position]] = [c.path_position] else: #if c.path_position < len(self.concat_string): assert c.path_position < len(self.concat_string),(c.path_position, len(self.concat_string)) seq2 = True ; car = self.concat_string[c.path_position] try: dic_car_seq2[car].append(c.path_position) except KeyError: dic_car_seq2[car] = [c.path_position] else: all_children_leaf = False if not all_children_leaf: continue if (not seq1) and (not seq2): continue if longueur_chaine >= min_size and len(n.path_indices) > 1: seq1 = seq2 = False ; lOccMEM = [] for car,lOcc in dic_car_seq1.iteritems(): if len(lOcc) == 1: lOccMEM.append(lOcc[0]) seq1 = True for car,lOcc in dic_car_seq2.iteritems(): if len(lOcc) == 1: lOccMEM.append(lOcc[0]) seq2 = True logging.debug(lOccMEM) if seq1 and seq2: for pi in lOccMEM: for pos in xrange(pi,pi+longueur_chaine): if (pi <= seq_repeat[pos][0] <= seq_repeat[pos][1] <= pi+longueur_chaine or seq_repeat[pos][0] <= pi <= seq_repeat[pos][1] <= pi+longueur_chaine): seq_repeat[pos] = (pi,pi+longueur_chaine) return seq_repeat class GeneralisedSuffixTree(object): """Classe proxy qui construit un GeneralisedSuffixTree Celui-ci est construit et on extrait les infos nécessaires pour trouver les MEM Ensuite on supprime l'objet GeneralisedSuffixTree, ce qui permet d'économiser beaucoup de mémoire""" def __init__(self,sequences): self.sequences = sequences self.start_positions = [0] self.concat_string = '' for i in xrange(len(sequences)): if chr(i+1) in sequences[i]: raise "The suffix tree string must not contain chr(%d)!"%(i+1) self.concat_string += sequences[i]+chr(i+1) self.start_positions += [len(self.concat_string)] self.start_positions += [self.start_positions[-1]+1] # empty string self.sequences += [''] self.st = TrueGeneralisedSuffixTree(sequences) def get_MEM_(self, min_size=1): """Renvoie un dico de tous les Maximal Exact Matches de taille min min_size Le dico est indexé par la taille de MEM qui renvoie une liste de toutes les positions de cette taille. Linéaire(nb de MEM) < linéaire(taille de la séquence)""" dic_repeat = self.st.get_repeat(min_size) del self.st #print dic_repeat seq_repeat = {} seq_repeat[-1] = -1 for i in xrange(len(self.sequences[0])+len(self.sequences[1])): seq_repeat[i] = i for length in dic_repeat: for pos in dic_repeat[length]: for pos_i in xrange(pos,pos+length): seq_repeat[pos_i] = max(seq_repeat[pos_i], pos+length) # à chaque position de seq_repeat doit correspondre la position # de la fin de la répétition la plus longue commençant à cette posistion #print seq_repeat assert len(seq_repeat) == len(self.sequences[0])+len(self.sequences[1])+1 dic_MEM = {} for debut in seq_repeat: fin = seq_repeat[debut] # si ça fait plus de 1 carac et si on ne vient pas d'ajouter la répet if (fin > debut and seq_repeat[debut-1] != fin): try: dic_MEM[fin-debut].append(debut) except KeyError: dic_MEM[fin-debut] = [debut] #print dic_MEM return dic_MEM def get_MEM_list(self, min_size=1, sr=1): """Renvoie un dico de tous les Maximal Exact Matches de taille min min_size Le dico est indexé par la taille de MEM qui renvoie une liste de toutes les positions de cette taille. Linéaire(nb de MEM) < linéaire(taille de la séquence)""" seq_repeat = self.get_seq_repeat(min_size) if sr==3: return seq_repeat elif sr == 1: dic_MEM = {} #for (debut,fin) pos = len(seq_repeat)-1 while pos >= 0: (debut,fin) = seq_repeat[pos] if fin > debut: try: dic_MEM[fin-debut].append(debut) except KeyError: dic_MEM[fin-debut] = [debut] pos = debut-1 #print dic_MEM return dic_MEM elif sr==4: dic_MEM = {} longueur_s1 = len(self.sequences[0]) texte = self.sequences[0] + self.sequences[1] pos = len(seq_repeat)-1 while pos >= 0: # attention, ne pas mélanger debut et pos_debut, 2 compteurs différents (debut,fin) = seq_repeat[pos] if fin > debut: longueur = fin - debut if debut > longueur_s1: pos_debut = debut - 1 else: pos_debut = debut t = hash(texte[pos_debut:pos_debut+longueur]) if not dic_MEM.has_key(longueur): dic_MEM[longueur] = {} try: dic_MEM[longueur][t].append(pos_debut) except KeyError: dic_MEM[longueur][t] = [pos_debut] pos = debut-1 #logging.debug( dic_MEM) #for longueur,dic_cle in dic_MEM.iteritems(): # for cle,lOcc in dic_cle.iteritems(): # logging.debug(str(longueur)+' / '+texte[lOcc[0]:lOcc[0]+longueur]+'/ '+str(lOcc)) return dic_MEM def get_MEM(self, min_size=1, sr=1): """Renvoie un dico de tous les Maximal Exact Matches de taille min min_size Le dico est indexé par la taille de MEM qui renvoie une liste de toutes les positions de cette taille. Linéaire(nb de MEM) < linéaire(taille de la séquence)""" seq_repeat_deb,seq_repeat_fin = self.st.get_seq_repeat(min_size) del self.st # permet d'économiser beaucoup de mémoire logging.debug("suffixTree deleted") if sr==3: return seq_repeat elif sr == 1: dic_MEM = {} #for (debut,fin) pos = len(seq_repeat_deb)-1 while pos >= 0: debut = seq_repeat_deb[pos] fin = seq_repeat_fin[pos] if fin > debut: try: dic_MEM[fin-debut].append(debut) except KeyError: dic_MEM[fin-debut] = [debut] pos = debut-1 #print dic_MEM return dic_MEM elif sr==4: dic_MEM = {} longueur_s1 = len(self.sequences[0]) texte = self.sequences[0] + self.sequences[1] pos = len(seq_repeat_deb)-1 while pos >= 0: # attention, ne pas mélanger debut et pos_debut, 2 compteurs différents debut = seq_repeat_deb[pos] fin = seq_repeat_fin[pos] if fin > debut: longueur = fin - debut if debut > longueur_s1: pos_debut = debut - 1 else: pos_debut = debut t = hash(texte[pos_debut:pos_debut+longueur]) if not dic_MEM.has_key(longueur): dic_MEM[longueur] = {} try: dic_MEM[longueur][t].append(pos_debut) except KeyError: dic_MEM[longueur][t] = [pos_debut] pos = debut-1 #logging.debug( dic_MEM) #for longueur,dic_cle in dic_MEM.iteritems(): # for cle,lOcc in dic_cle.iteritems(): # logging.debug(str(longueur)+' / '+texte[lOcc[0]:lOcc[0]+longueur]+'/ '+str(lOcc)) return dic_MEM def get_MEM_index_chaine2(self, min_size=1, sr = 1): """ Modifie les coordonnées de la 2e chaine en retirant 1 à cause du séparateur pour le suffix tree """ if sr == 4: seq_repeat = self.get_MEM(min_size, sr=4) logging.log(5,'get_MEM done') return seq_repeat elif sr == 3: seq_repeat = self.get_MEM(min_size, sr=3) logging.log(5,'get_MEM done') #print seq_repeat[:100] seq2 = [] longueur_s1 = len(self.sequences[0]) pos = len(seq_repeat)-1 for pos in xrange(len(seq_repeat)-1): deb,fin = seq_repeat[pos] if pos > longueur_s1: # 2e chaine seq2.append((deb-1,fin-1)) elif pos == longueur_s1: # position du séparateur des 2 séquences continue else: # première chaine seq2.append((deb,fin)) #print seq2 assert len(seq2) == len(self.sequences[0]) + len(self.sequences[1]) # renvoie une séquence de la longueur des 2 chaines concaténées # dont chaque item est une paire(debut du bloc, fon du bloc) return seq2 else: # sr 1 dic_repeat = self.get_MEM(min_size, sr=1) logging.log(5,'get_MEM done') #print seq_repeat[:100] dic_repeat2 = {} longueur_s1 = len(self.sequences[0]) texte = self.sequences[0] + self.sequences[1] for longueur,liste in dic_repeat.iteritems(): if longueur >= min_size: liste2 = [] for deb in liste: if deb > longueur_s1: liste2.append(deb-1) else: liste2.append(deb) for deb in liste2: t = texte[deb:deb+longueur] try: #dic_repeat2[t].append(deb) bisect.insort_right(dic_repeat2[t],deb) except KeyError: dic_repeat2[t] = [deb] del dic_repeat for t,l in dic_repeat2.iteritems(): #logging.debug((t,l)) assert len(t) >= min_size return dic_repeat2 def get_MEM_index_chaine3(self, min_size=1, just_keep_words=True, sr = 4, eliminRecouv=True): """ just_keep_words: si Vrai, on rogne les homologies de façon à n'avoir que des mots (ou suites de mots) entiers renvoie un dico indexé par (cle,longueur) ou cle représente hash(chaine) la chaine dont on fait référence et la longueur de la chaine, la valeur et la liste d'occurences de la chaine""" seq = self.get_MEM_index_chaine2(min_size, sr) #print seq logging.log(5,'get_MEM_index_chaine2 done') if sr == 4: #logging.debug('Recouv 4 eliminrecouv') texte = self.sequences[0]+self.sequences[1] for longueur, dicoOcc in seq.iteritems(): for cle_hash, lOcc in dicoOcc.iteritems(): #for occ in lOcc: # logging.debug(texte[occ:occ+longueur]) logging.debug(texte[lOcc[0]:lOcc[0]+longueur]) #print '====' if eliminRecouv: a = time.time() recouv = recouvrement.Recouvrement4(self.sequences[0]+self.sequences[1],seq,len(self.sequences[0]),min_size) #recouv = recouvrement.RecouvrementNaif(self.sequences[0]+self.sequences[1],seq,len(self.sequences[0]),min_size) #trace('blocs = recouv.eliminer_recouvrements()',locals()) blocs = recouv.eliminer_recouvrements() b = time.time() else: blocs = {} #blocs.NOSMEM_nb_bloc = 0 texte = self.sequences[0]+self.sequences[1] for longueur, dicoOcc in seq.iteritems(): for cle_hash, lOcc in dicoOcc.iteritems(): for occ in lOcc: cle = hash(texte[occ:occ+longueur]) try: blocs[(cle,longueur)].append(occ) except KeyError: blocs[(cle,longueur)] = [occ] #blocs.NOSMEM_nb_bloc += 1 #for cle,lOcc in blocs.iteritems(): # lOcc.reverse() # on reverse car dans la suite dans l'algo c'est nécessaire #for (cle,longueur),liste_pos in blocs.iteritems(): # for occ in liste_pos: # print texte[occ:occ+longueur] #print '=========================' #logging.debug('Recouv 4 FIN eliminrecouv') if eliminRecouv and recouv.SMEM_nb_bloc >= 100: logFile = os.path.join(os.getcwd(), 'resultat_complex_recouv.csv') if not os.path.exists(logFile): #le fichier csv n'existe pas, on écrit sur la premiere ligne le nom des params csv = open(logFile, 'w') #mode creation csv.write("SMEM_nb_bloc;SMEM_nb_occ;SMEM_tot_size;NOSMEM_nb_bloc;NOSMEM_nb_occ;NOSMEM_tot_size;temps;A;B\n") else: csv = open(logFile, 'a') #mode append csv.write(str(recouv.SMEM_nb_bloc)+';'+str(recouv.SMEM_nb_occ)+';'+str(recouv.SMEM_tot_size)+';'+ str(recouv.NOSMEM_nb_bloc)+';'+str(recouv.NOSMEM_nb_occ)+';'+str(recouv.NOSMEM_tot_size)+';'+ str(b-a)+';'+str(len(self.sequences[0]))+';'+str(len(self.sequences[1]))+'\n') csv.close() elif sr == 3: recouv = recouvrement.Recouvrement3(self.sequences[0]+self.sequences[1],seq,len(self.sequences[0])) blocs = recouv.eliminer_recouvrements() else: # 1 #logging.debug('Recouv 1 eliminrecouv') recouv = recouvrement.Recouvrement(self.sequences[0]+self.sequences[1],seq,len(self.sequences[0]),min_size) blocs1 = recouv.eliminer_recouvrements() #blocs1 = {} #trace('blocs1 = recouv.eliminer_recouvrements()',locals()) #logging.info('Recouv 1 FIN eliminrecouv') blocs = {} for texte,liste in blocs1.iteritems(): cle = hash(texte) longueur = len(texte) # ATTENTION, on doit inverser ici SSI on inverse présédemment # dans get_MEM_index_chaine2(sr=1) lors de l'ajout dans dic_repeat2 liste.reverse() blocs[(cle,longueur)] = liste del blocs1 #print 'blocs:'+str(blocs) logging.log(5,'eliminer_recouvrements done') dic_chaine2 = {} longueur_s1 = len(self.sequences[0]) #print len(blocs),longueur_s1,len(self.sequences[0])+len(self.sequences[1]),min_size # ATTENTION le 1er espace (début de liste) est l'espace habituel # le second (fin de liste) est ALT+0160 qui est un espace insécable inséré, en particulier, # par Word avant les ponctuations faibles : ; ! ? et guillemets sep = """. !\r,\n:\t;-?"'`’() """ #liste des séparateurs # on ne garde que les chaine de taille > min et répétées for (cle,longueur),liste_pos in blocs.iteritems(): if (longueur >= min_size and len(liste_pos) >= 2): # on inverse car les pos ont été ajoutées dans l'ordre décroissant dans get_MEM liste_pos.reverse() if (liste_pos[0] < longueur_s1 <= liste_pos[-1]): if just_keep_words: # basé sur seulement 2 homologies dans la liste, quid du comportement avec + d'homologies ? # car pour couper une homologies, on regarde les caractères précédents et suivants # dans les 2 chaines, la premiere de la liste (texte 1) et la dernière (texte 2) #logging.debug( '1:' +self.sequences[0][liste_pos[0]:liste_pos[0]+longueur]+'/'+str(liste_pos)) # les car de début et fin des chaines doivent êrre égaux assert self.sequences[0][liste_pos[0]] == self.sequences[1][liste_pos[-1]-longueur_s1] assert self.sequences[0][liste_pos[0]+longueur-1] == self.sequences[1][liste_pos[-1]+longueur-1-longueur_s1], chaine+' / '+self.sequences[0][liste_pos[0]:liste_pos[0]+longueur]+' / '+self.sequences[1][liste_pos[-1]-longueur_s1:liste_pos[-1]+longueur-longueur_s1]+' / *'+self.sequences[0][liste_pos[0]+longueur]+'* / *'+self.sequences[1][liste_pos[-1]+longueur-longueur_s1]+'**' # recherche du premier séparateur dans la chaine i = liste_pos[0] ; i2 = liste_pos[-1] if (i == 0 or self.sequences[0][i-1] in sep) and (i2-longueur_s1 == 0 or self.sequences[1][i2-1-longueur_s1] in sep): # soit début de séquence et dans ce cas le car précédent est forcément un sep car on est en mode mot # soit caractère précédent est un séparateur decalage_avant = 0 else: # sinon on cherche le 1er sep dans la chaine courante while (i < liste_pos[0]+longueur and self.sequences[0][i] not in sep and self.sequences[1][i2-longueur_s1] not in sep): i += 1 ; i2 += 1 decalage_avant = i - liste_pos[0] assert i >= liste_pos[0] # recherche du séparateur le + à droite dans la chaine j = liste_pos[0]+longueur-1 ; j2 = liste_pos[-1]+longueur-1 if ((j == len(self.sequences[0])-1 or self.sequences[0][j+1] in sep) and (j2-longueur_s1 == len(self.sequences[1])-1 or self.sequences[1][j2+1-longueur_s1] in sep)): pass # idem que i: car de fin de séquence ou car precedent sep else: while (j >= liste_pos[0] and self.sequences[0][j] not in sep and j2 >= liste_pos[-1] and self.sequences[1][j2-longueur_s1] not in sep): j -= 1 ; j2 -= 1 assert j <= liste_pos[0]+longueur-1 #logging.debug('decalage_avanti='+str(decalage_avant)+', longueur2='+str(j+1-i)) if i < j: # si la nouvelle chaine n'est pas vide longueur2 = j+1 -i ; assert longueur2 <= longueur if (longueur2 >= min_size): cle2 = hash(self.sequences[0][i:j+1]) tt = [x+decalage_avant for x in liste_pos] dic_chaine2[(cle2,longueur2)] = tt #logging.debug('2:'+self.sequences[0][tt[0]:tt[0]+longueur2]+'/'+str(tt)) else: # cas standard dic_chaine2[(cle,longueur)] = liste_pos #print 'dic_chaine2:'+str(dic_chaine2) #print len(dic_chaine2)#,dic_chaine2 logging.log(10,'len(blocs)='+str(len(blocs))+' / len(dic_chaine2)='+str(len(dic_chaine2)) +' / taille chaines dic_chaine2='+str(sum([longueur*len(li) for (c,longueur),li in dic_chaine2.iteritems()]))) texte = self.sequences[0] + self.sequences[1] for (cle,longueur),lOcc in dic_chaine2.iteritems(): logging.debug(texte[lOcc[0]:lOcc[0]+longueur]) return dic_chaine2 def get_MEM_index_chaine(self, min_size=1): """ Converti le dic de sortie pour le module MediteAppli index par la chaine directement qui donne les positions de cette chaine""" dic_MEM = self.get_MEM(min_size) texte = self.sequences[0] + '$'+ self.sequences[1] longueur_s1 = len(self.sequences[0]) dic_chaine = {} for longueur,liste_pos in dic_MEM.iteritems(): if len(liste_pos) < 2: continue for pos in liste_pos: chaine = texte[pos:pos+longueur] # à cause du $ dans la 2e chaine, tous les indices sont décalés de 1, donc on les MAJ correctement if pos > longueur_s1: pos -= 1 try: dic_chaine[chaine].append(pos) except KeyError: dic_chaine[chaine] = [pos] #print dic_chaine dic_chaine2 = {} for chaine,liste_pos in dic_chaine.iteritems(): if (len(chaine) >= min_size and len(liste_pos) >= 2): # on inverse car les pos ont été ajoutées dans l'ordre décroissant dans get_MEM liste_pos.reverse() if (liste_pos[0] < longueur_s1 <= liste_pos[-1]): dic_chaine2[chaine] = liste_pos return dic_chaine2 def trace(commande,dic_locals): #import hotshot,hotshot.stats,sys import profile,pstats,sys profile.runctx(commande,globals(), dic_locals,'c:\workspace\medite\statFile') s = pstats.Stats('c:\workspace\medite\statFile') s.sort_stats('time') s.print_stats() #sys.stderr.flush() #sys.stdout.flush() #prof = hotshot.Profile("c:\workspace\medite\statFile") #benchtime = prof.runctx(commande,globals(), dic_locals) #prof.close() #stats = hotshot.stats.load("c:\workspace\medite\statFile") #stats.strip_dirs() #stats.sort_stats('time', 'calls') #stats.print_stats(50) def simple_test(): print 'SIMPLE TEST' st = SuffixTree('mississippi') for l in st.leaves: print 'leaf:', '"'+l.path_label+'"', ':', '"'+l.edge_label+'" path_position '+ \ str(l.path_position)+' / edge_label_start '+str(l.edge_label_start)+ \ ' / edge_label_end '+str(l.edge_label_end) for n in st.inner_nodes: print 'inner:', '"'+n.edge_label+'" path_position '+ \ str(n.path_position)+' / edge_label_start '+str(n.edge_label_start)+ \ ' / edge_label_end '+str(n.edge_label_end) print 'done.\n\n' def generalised_test(): print 'GENERALISED TEST' sequences = ['xabxa','babxba'] st = GeneralisedSuffixTree(sequences) for shared in st.shared_substrings(): print '-'*70 for seq,start,stop in shared: print seq, '['+str(start)+':'+str(stop)+']', print sequences[seq][start:stop], print sequences[seq][start:] print '='*70 for shared in st.shared_substrings(2): print '-'*70 for seq,start,stop in shared: print seq, '['+str(start)+':'+str(stop)+']', print sequences[seq][start:stop], print sequences[seq][start:] print '='*70 print 'done.\n\n' def test3(): print 'TEST3' #sequences = ['le chat est grand','1 chat est vert'] #['le chat est grand','un grand chat est vert'] sequences = ["Le 7 décembre, à onze heures du matin, et à deux heures et demie du soir, en pinçant la patte dont les vaisseaux sont liés, on y observe des mouvements réflexes. On laisse la grenouille dans les mêmes conditions que la veille.", "Le 7 décembre, à 11 h. du matin et à 2 h. 1/2, en pinçant la patte antérieure et la patte liée, on a des mouvements réflexes dans cette dernière."] st = GeneralisedSuffixTree(sequences) # st.shared_substrings2() texte=sequences[0]+sequences[1] for (start,stop) in st.shared_substrings2(): print '['+str(start)+':'+str(stop)+']', print "*"+texte[start:stop]+ "* \n ", print '='*70 def test(): simple_test() #generalised_test() #test3() if __name__ == '__main__': test()