import os, string, math, re, io
import seq_io, tools, progressbar
try:
    import psyco
    psyco.profile()
except:
    pass
    
#####################################################################################################
class OUP:
    def __init__(self,pattern_type,task,path,query,subject=None):
        # Task to perform
        self.task = task
        # Pattern type
        self.pattern_type = pattern_type
        # Path to source files
        self.path = path
        # Source file or sequence
        self.query = query
        # Source file, or  sequence, or name of OUP database, or None
        self.sbjct = subject
        self.index_file = ""
        # I/O library
        self.IO = seq_io.IO()
        # Report object
        self.oReport = None
    
    def execute(self):
        if re.match("^n[0-6]_[1-7]mer",self.pattern_type):
            self._generate_patterns()
    
    def tostring(self):
        if self.oReport:
            return self.oReport.tostring()
        else:
            return ""
        
    def get_distance_matrix(self,lookup=[]):
        if self.oReport:
            return self.oReport.distance_matrix(lookup)
        else:
            return []
    
    def get_OUV(self,lookup=[]):
        if self.oReport:
            return self.oReport.OUV_matrix(lookup)
        else:
            return []
    
    def svg(self):
        if self.oReport:
            return self.oReport.svg()
        else:
            return ""
    
    def getPattern(self,seq,ptype,title=""): # ptype -> n1_4mer
        ptype,norm,wlength = self._parse_ptype(ptype)
        oPattern = Pattern(wlength,title)
        oPattern.setPattern(seq,norm)
        return oPattern
    
    def _parse_ptype(self,ptype):
        norm = int(ptype[1:ptype.find("_")])
        wlength = int(ptype[ptype.find("_")+1:ptype.find("mer")])
        ptype = ptype[0]
        return ptype,norm,wlength

    def _generate_patterns(self):
        query_list = {}
        if not self.query:
            return
        for query_file in self.query.split(","):
            query,children = self._get_seqlist(os.path.join(self.path,query_file))
            print()
            print(("Calculating query %s patterns in process..." % self.pattern_type))
            bar = progressbar.indicator(len(query),os.path.basename(query_file))
            counter = 1
            for query_name in query:
                oPattern = self.getPattern(query[query_name],self.pattern_type,query_name)
                for child_name in children:
                    oChildPattern = self.getPattern(children[child_name],self.pattern_type,child_name)
                    oPattern.add_child(oChildPattern)
                query_list[query_name] = oPattern
                bar(counter)
                counter += 1
            bar.stop()
        sbjct_list = {}
        if self.sbjct:
            for sbjct_file in self.sbjct.split(","):
                sbjct,children = self._get_seqlist(os.path.join(self.path,sbjct_file))
                print()
                print(("Calculating sbjct %s patterns in process..." % self.pattern_type))
                bar = progressbar.indicator(len(sbjct),os.path.basename(sbjct_file))
                counter = 1
                for sbjct_name in sbjct:
                    oPattern = self.getPattern(sbjct[sbjct_name],self.pattern_type,sbjct_name)
                    for child_name in children:
                        oChildPattern = self.getPattern(children[child_name],self.pattern_type,child_name)
                        oPattern.add_child(oChildPattern)
                    sbjct_list[sbjct_name] = oPattern
                    bar(counter)
                    counter += 1
                bar.stop()
        if self.task == "pattern":
            self.oReport = PlotOUPattern(self.pattern_type,query_list)
        elif self.task == "2D-plot":
            self.oReport = Plot2D(self.pattern_type,query_list,sbjct_list)
        elif self.task == "3D-plot":
            self.oReport = Plot3D(self.pattern_type,query_list,sbjct_list)
        elif self.task == "subtraction":
            self.oReport = WordSubtraction(self.pattern_type,query_list,sbjct_list)
        elif self.task == "d-matrix" or self.task == "s-matrix" :
            self.oReport = DistanceMatrix(self.pattern_type,query_list,sbjct_list)
        
    def _parse_ptype(self,ptype):
        pt = ptype[0]
        norm = int(ptype[1:ptype.find("_")])
        wlength = int(ptype[ptype.find("_")+1:ptype.find("mer")])
        return pt,norm,wlength
        
    def _get_seqlist(self,source):
        seqlist = None
        children = {} # subsequences specified as sequence [start1-stop1;start2-stop2;start3-;-stop4;...]
        if type(source) == type(io.StringIO()):
            source = source.read()
        ch = source.find(" [")
        if ch > -1 and source[-1] == "]":
            for child_name in source[ch+2:source.find("]")].split(";"):
                children[child_name] = ""
            source = source[:ch]
        if source and not os.path.exists(source):
            seqlist = {"Query":source}
        elif source and source[source.rfind("."):].upper() in (".GBK",".GB"):
            seqlist = self.IO.openGBK(source)[0]
        elif source and source[source.rfind("."):].upper() in (".FASTA",".FST",".FAS",".FNA",".FSA"):
            seqlist = self.IO.openFasta(source)[0]
        if children and len(seqlist)==1:
            seqname = list(seqlist.keys())[0]
            sequence = seqlist[seqname]
            for child_name in list(children.keys()):
                del children[child_name]
                d = child_name.find("-")
                if d == -1:
                    continue
                elif d == 0:
                    child_name = "1"+child_name
                elif d == len(list(child_name))-1:
                    child_name += str(len(list(sequence)))
                try:
                    start,stop = [int(s) for s in child_name.split("-")]
                except:
                    continue
                if stop <= start:
                    continue
                if start < 1:
                    start = 1
                if stop > len(list(sequence)):
                    stop = len(list(sequence))
                children["%s [%d-%d]" % (seqname,start,stop)] = sequence[start-1:stop-1]
        else:
            children = {}
                
        return seqlist,children
        
    def _check_dbname(self):
        dbname = os.path.basename(self.sbjct)
        if dbname in self.shortnames:
            dbname = self.shortnames[dbname]
            self.sbjct = os.path.join(self.path,dbname)
        if dbname not in list(self.shortnames.values()):
            return
        try:
            index = self.IO.openDBFile(self.index_file)[1]
        except:
            return
        if not index:
            return
        filelist = os.listdir(self.path)
        tbname,pattern_type = dbname.split(":")
        if (pattern_type not in index['tables'] or 
            tbname not in index['tables'][pattern_type] or
            index['tables'][pattern_type][tbname] not in filelist):
            return
        return index

######################################################################################################
class Report: # The report of a search through OUP_databases
    def __init__(self):
        self.container = []
        
    def __len__(self):
        return len(self.container)
    
    def __getitem__(self,key):
        for record in self.container:
            if record.title == key:
                return record
        return 
    
    def __contains__(self,key):
        for record in self.container:
            if record.title == key:
                return True
        return False
    
    def __iter__(self):
        records = self.get()
        iterator = iter(records)
        return iterator
    
    def __delitem__(self,key):
        for i in range(len(self.container)):
            if self.container[i].title == key:
                del self.container[i]
                return
    
    def add(self,item):
        self.container.append(item)
    
    def get(self,sorting=None):
        records = []
        for record in self.container:
            records.append(record.copy())
        return records

######################################################################################################
class pattern_collection:
    def __init__(self,ptype,first,second=None):
        self._query_array = {}
        if first:
            self._query_array.update(first)
        self._sbjct_array = {}
        if second:
            self._sbjct_array.update(second)
        self.ptype = ptype
        
    def __len__(self):
        return len(self._query_array)+len(self._sbjct_array)
    
    def __getitem__(self,key):
        if key in self._query_array:
            return self._query_array[key]
        elif key in self._sbjct_array:
            return self._sbjct_array[key]
        return 
    
    def __contains__(self,key):
        if key in self._query_array or key in self._sbjct_array:
            return True
        return False
    
    def __iter__(self):
        iterator = iter(list(self._query_array.keys())+list(self._sbjct_array.keys()))
        return iterator
    
    def __delitem__(self,key):
        if key in self._query_array:
            del self._query_array[key]
        if key in self._sbjct_array:
            del self._sbjct_array[key]
    
    def _get_pattern(self,seq,ptype,title=""): # ptype -> n1_4mer
        ptype,norm,wlength = parse_pattern_type(ptype)
        oPattern = Pattern(wlength,title)
        oPattern.setPattern(seq,norm)
        return oPattern
    
    def _find_pattern(self,pattern):
        for key in self:
            if pattern == key:
                return self[key]
            for child in self[key].child_patterns:
                if child.title == pattern:
                    return child
        return
    
    def _svg_arrow(self,x,y,alpha,width,color="black"):
        ax = x + 2.0*width*math.cos(alpha)
        ay = y - 2.0*width*math.sin(alpha)
        bx = x - 5.0*width*math.cos(alpha) - 3.0*width*math.sin(alpha)
        by = y + 5.0*width*math.sin(alpha) - 3.0*width*math.cos(alpha)
        cx = x - 3.0*width*math.cos(alpha)
        cy = y + 3.0*width*math.sin(alpha)
        dx = x - 5.0*width*math.cos(alpha) + 3.0*width*math.sin(alpha)
        dy = y + 5.0*width*math.sin(alpha) + 3.0*width*math.cos(alpha)
        return ("<path style=\"stroke:%s;stroke-width:%f;fill:%s\" d=\"M%f,%fL%f,%fL%f,%fL%f,%fZ\" />" %
                    ("black",1.0,"black",ax,ay,bx,by,cx,cy,dx,dy))
    
    def _format_title(self,title,space): # space - number of characters per line
        words = title.split(" ")
        title = [""]
        for word in words:
            title[-1] += word+" "
            if len(list(title[-1])) >= space:
                title.append("")
        if len(list(title)) > 3:
            title[2] += " ..."
            title = title[:3]
        return title
    
    def do(self,command,ArgList=None):
        if command == "Get pattern":
            return self.get_pattern(ArgList)
        elif command == "Get distance" and len(ArgList)==2:
            key1,key2 = ArgList
            return self.get_distance(key1,key2)
    
    def keys(self,flg=0):
        if flg == 1:
            return list(self._query_array.keys())
        elif flg == 2:
            return list(self._sbjct_array.keys())
        else:
            return list(self._query_array.keys())+list(self._sbjct_array.keys())
    
    def get_pattern(self,key):
        if key in self:
            return self[key]
        
    def get_patterns(self,lookup=[],flg=0):
        patterns = []
        for key in self.keys(flg):
            if lookup and key not in lookup:
                continue
            patterns.append(self[key])
        return patterns
    
    def get_parent_and_child_patterns(self):
        parent_patterns = self.get_patterns()
        child_patterns = []
        for oPattern in parent_patterns:
            if oPattern.child_patterns:
                child_patterns.extend(oPattern.child_patterns)
        return parent_patterns+child_patterns
    
    def get_wlength(self):
        ptype,norm,wlength = parse_pattern_type(self.ptype)
        return wlength
                    
    def get_normalization(self):
        ptype,norm,wlength = parse_pattern_type(self.ptype)
        return norm
                    
    def get_ptype(self):
        ptype,norm,wlength = parse_pattern_type(self.ptype)
        return ptype
                    
    def get_distance(self,first="",second=""):
        if first in self:
            first = self[first]
        else:
            first = self._find_pattern(first)
        if not first:
            return
        if second in self:
            second = self[second]
        else:
            second = self._find_pattern(second)
        if not second:
            return
        return first-second
    
    def get_local_distances(self,key="",child_name=""):
        if (key and key not in self) or (child_name and not self[key].has_child(child_name)):
            return []
        if not key:
            if self._query_array:
                key = self._query_array[list(self._query_array.keys())[0]]
            elif self._sbjct_array:
                key = self._sbjct_array[list(self._sbjct_array.keys())[0]]
            else:
                return []
        distances = []
        for key in self:
            distances.append([key])
            for child in self[key].child_patterns:
                distances[-1].append(self[key]-child)
        return distances
    
    def variance(self,key=""):
        if key and key not in self:
            return
        elif key:
            return self[key].getOUV()
        else:
            return self._query_array[list(self._query_array.keys())[0]].getOUV()
    
    def pattern_skew(self,key):
        if key and key not in self:
            return
        elif key:
            return self[key].getPS()
        else:
            return self._query_array[list(self._query_array.keys())[0]].getPS()
    
    def tostring(self,mode=""): # mode = "distance table"/"pattern table"/"word table"
        return ""
        if mode == "distance table" or (not mode and len(self) > 1):
            return self.distance_table()
        elif mode == "pattern table" or (not mode and len(self)==1):
            return self.pattern_table()
        elif mode == "word table":
            return self.word_table()
    
    def distance_table(self,lookup=[]):
        if not len(self):
            return ""
        print()
        print("Calculating OU distance matrix in process...")
        abbreviations = {}
        all_patterns = self.get_parent_and_child_patterns()
        if lookup:
            n = len(lookup)
            titles = [p.title for p in all_patterns]
            for i in range(len(lookup)-1,-1,-1):
                if lookup[i] not in titles:
                    del lookup[i]
            if not lookup:
                return
            lookup = [filter(lambda p: p.title==title, all_patterns)[0] for title in lookup]
        else:
            n = len(all_patterns)
            lookup = all_patterns
        output = '\t%d\n' % n
        lines = []
        bar = progressbar.indicator(n**2,"OUP Dist.")
        counter = 1
        for query in lookup:
            line = []
            key = abbreviate(query.title,list(abbreviations.keys()))
            abbreviations[key] = query.title
            line.append(key + " "*(10-len(list(key))))
            for subject in lookup:
                if lookup and subject.title not in lookup:
                    continue
                if subject.title == query.title:
                    line.append(0)
                else:
                    line.append(query-subject)
            bar(counter*n)
            counter += 1
            lines.append(line)
        bar.stop()
        output += "\n".join([" ".join([format_number(n,2) for n in s]) for s in lines])+"\n\nAbbreviations:\n"
        for key in abbreviations:
            output += "\t%s:\t%s\n" % (key,abbreviations[key])
        return output

    def distance_matrix(self,lookup=[]):
        if not len(self):
            return []
        print()
        print("Calculating OU distance matrix in process...")
        all_patterns = self.get_parent_and_child_patterns()
        if lookup:
            n = len(lookup)
            titles = [p.title for p in all_patterns]
            for i in range(len(lookup)-1,-1,-1):
                if lookup[i] not in titles:
                    del lookup[i]
            if not lookup:
                return
            lookup = [filter(lambda p: p.title==title, all_patterns)[0] for title in lookup]
        else:
            n = len(all_patterns)
            lookup = all_patterns
        seqlist = []
        output = tools.matrix([n,n])
        bar = progressbar.indicator(n**2,"OUP Dist.")
        counter = 1
        for i in range(len(lookup)):
            query = lookup[i]
            seqlist.append(query.title)
            for j in range(len(lookup)):
                subject = lookup[j]
                if subject.title == query.title:
                    output[i][j] = 0
                else:
                    output[i][j] = (query-subject)
            print (output,query.ouv,subject.ouv)
            bar(counter*n)
            counter += 1
        bar.stop()
        return output,seqlist

    def OUV_matrix(self,lookup=[]):
        if not len(self):
            return []
        all_patterns = self.get_parent_and_child_patterns()
        if lookup:
            n = len(lookup)
            titles = [p.title for p in all_patterns]
            for i in range(len(lookup)-1,-1,-1):
                if lookup[i] not in titles:
                    del lookup[i]
            if not lookup:
                return
            lookup = [filter(lambda p: p.title==title, all_patterns)[0] for title in lookup]
        else:
            n = len(all_patterns)
            lookup = all_patterns
        seqlist = []
        variances = []
        for query in lookup:
            variances.append(query.getOUV())
            seqlist.append(query.title)
        return variances,seqlist

    def sort_dist_table(self,a,b):
        return int(sum(b[1:])-sum(a[1:]))

######################################################################################################
class PlotOUPattern(pattern_collection):
    def __init__(self,ptype,query):
        pattern_collection.__init__(self,ptype,query)
    
    def tostring(self):
        output = []
        for seqname in self._query_array:
            if not self._query_array[seqname].child_patterns:
                output.append(self._query_array[seqname].tostring())
            else:
                for child in self._query_array[seqname].child_patterns:
                    output.append(child.tostring())
        return "\n\n".join(output)
    
    def svg(self):
        if len(self) == 0:
            svg = ""
        else:
            svg,width,height = self.svg_plot()
        return "\n".join(["<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % (width,height,width,height),
                svg,"</svg>"])
        
    def svg_plot(self,width=0,height=0,flg_finish=False):
        series = 0
        svg = [] 
        if not width or width < 400:
            width = 400
        if height != width:
            height = width
        for seqname in self._query_array:
            if not self._query_array[seqname].child_patterns:
                svg += self._query_array[seqname].svg(width,height,0,series*height,False)
                series += 1
            else:
                for child in self._query_array[seqname].child_patterns:
                    svg += child.svg(width,height,0,series*height,False)
                    series += 1
        if flg_finish:
            svg.insert(0,"<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % (width+300,height*series,width+300,height*series))
            svg.append("</svg>")
        return "\n".join(svg),width,height*series

######################################################################################################
class Plot2D(pattern_collection):
    def __init__(self,ptype,first,second=None):
        pattern_collection.__init__(self,ptype,first,second)
        self.references = []
        
    def get_cross_reference(self):
        # cross_reference = [['query_name','sbjct_name',dist],
        #                    [query_name,query - query.child_1,...],
        #                    [query_name,query - sbjct.child_1,...],
        #                    [sbjct_name,sbjct - sbjct.child_1,...],
        #                    [sbjct_name,sbjct - query.child_1,...],
        #                    [query_child_1 - sbjct_child_1, query_child_1 - sbjct_child_2,...]]
        references = []
        if not self._query_array:
            return references
        query_name = sbjct_name = ""
        for query_name in self._query_array:
            cross_reference = []
            query = self._query_array[query_name]
            cross_reference.append([query_name])
            for child in query.child_patterns:
                cross_reference[-1].append(query-child)
            for sbjct_name in self._sbjct_array:
                sbjct = self._sbjct_array[sbjct_name]
                cross_reference.insert(0,[query_name,sbjct_name,query-sbjct])
                cross_reference.append([query_name])
                for child in sbjct.child_patterns:
                    cross_reference[-1].append(query-child)
                cross_reference.append([sbjct_name])
                for child in sbjct.child_patterns:
                    cross_reference[-1].append(sbjct-child)
                cross_reference.append([sbjct_name])
                for child in query.child_patterns:
                    cross_reference[-1].append(sbjct-child)
                cross_reference.append([])
                for query_child in query.child_patterns:
                    for sbjct_child in sbjct.child_patterns:
                        cross_reference[-1].append(query_child-sbjct_child)
            if len(cross_reference) == 1:
                cross_reference.insert(0,[query_name,"",0])
                cross_reference.append([query_name])
                for i in range(2):
                    cross_reference.append([""])
                cross_reference.append([])
            references.append(cross_reference)
        return references
    
    def tostring(self):
        if len(self) == 0:
            output = ""
        else:
            output = self.get_table()
        return output
    
    def get_table(self):
        if not self.references:
            self.references = self.get_cross_reference()
        if not self.references:
            return ""
        output = []
        for data in self.references:
            query_name = data[0][0]
            sbjct_name = data[0][1]
            D = data[0][2]
            output.append("\t".join(["",query_name,sbjct_name]))
            output.append("\t".join([query_name,"0",str(D)]))
            output.append("\t".join([sbjct_name,str(D),"0"]))
            query_internal_distances = data[1][1:]
            query_external_distances = data[2][1:]
            sbjct_internal_distances = data[3][1:]
            sbjct_external_distances = data[4][1:]
            for i in range(len(sbjct_internal_distances)):
                d1 = query_external_distances[i]
                d2 = sbjct_internal_distances[i]
                output.append("\t".join([query_name+":"+str(i+1),str(d1),str(d2)]))
            for i in range(len(query_internal_distances)):
                d1 = query_internal_distances[i]
                d2 = sbjct_external_distances[i]
                output.append("\t".join([sbjct_name+":"+str(i+1),str(d1),str(d2)]))
            output.append("")
        return "\n".join(output)
    
    def svg(self,mode="",width=0,height=0):
        if len(self) == 0:
            svg = ""
        else:
            svg,width,height = self.svg_plot()
        return "\n".join(["<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % (width,height,width,height),
                svg,"</svg>"])
        
    def svg_plot(self,width=0,height=0,flg_finish=False):
        if not self.references:
            self.references = self.get_cross_reference()
        if not self.references:
            return "",width,height
        if not width:
            width = 750
        if not height:
            height = 500
        svg = []
        x_indend = 30
        y_indend = 30
        if width <= height:
            span100percent = width-x_indend-30
        else:
            span100percent = height-y_indend-30
        X = span100percent/2 + x_indend + 5
        Y = 30
        title_zone = 200
        series = 1
        for data in self.references:
            query_name = data[0][0]
            sbjct_name = data[0][1]
            D = data[0][2]
            query_internal_distances = data[1][1:]
            query_external_distances = data[2][1:]
            sbjct_internal_distances = data[3][1:]
            sbjct_external_distances = data[4][1:]
            x0 = X+x_indend
            y0 = Y+height-y_indend-title_zone
            x1 = X+x_indend+span100percent*D/100.0
            y1 = Y+height-y_indend-title_zone
            r = span100percent*D/200.0
            # add discriminating stream-line
            svg.append("<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" fill=\"none\" stroke=\"%s\" stroke-width=\"%f\" />" %
                (X+x_indend+r,Y+20,X+x_indend+r,Y+height-title_zone,"grey",.5))
            if D and D < 25:
                opacity = 0.6
                r = span100percent*.1
            else:
                opacity = 1.0
                r = span100percent*D/200.0
            # add query shadow areas
            svg.append("<circle style=\"stroke:%s;fill:%s;stroke-width:%f\" cx=\"%f\" cy=\"%f\" r=\"%f\" />" % 
                ("grey","#CCFFCC",1.0,x0,y0,r))
            svg.append("<circle style=\"stroke:%s;fill:%s;stroke-width:%f\" cx=\"%f\" cy=\"%f\" r=\"%f\" />" % 
                ("none","#66FF66",1.0,x0,y0,span100percent*.075))
            # add subject shadow areas
            svg.append("<circle style=\"stroke:%s;fill:%s;stroke-width:%f;opacity:%f\" cx=\"%f\" cy=\"%f\" r=\"%f\" />" % 
                ("grey","#CCFFCC",1.0,opacity,x1,y1,r))
            svg.append("<circle style=\"stroke:%s;fill:%s;stroke-width:%f;opacity:%f\" cx=\"%f\" cy=\"%f\" r=\"%f\" />" % 
                ("none","#66FF66",1.0,opacity,x1,y1,span100percent*.075))
            # add axes
            svg.append(self._svg_cartesian(X,Y,width,height-title_zone,x_indend,y_indend))
            # add 100% mark
            svg.append("<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" fill=\"none\" stroke=\"%s\" stroke-width=\"%f\" />" %
                (X+x_indend+span100percent,Y+height-y_indend-title_zone,
                X+x_indend+span100percent,Y+height-y_indend-title_zone+6,"black",3.0))
            svg.append("<text x=\"%d\" y=\"%d\" text-anchor=\"middle\">%s</text>" % 
                (X+x_indend+span100percent,Y+height-y_indend-title_zone+22,"100%"))
            # add spots
            title = "D = %f;" % D
            svg.append(self._svg_spot(x0,y0,"#009900",title))
            svg.append(self._svg_spot(x1,y1,"#009900",title))
            for i in range(len(sbjct_internal_distances)):
                d1 = query_external_distances[i]
                d2 = sbjct_internal_distances[i]
                x = (d1**2 - d2**2 + D**2)/2.0/D
                y = y0-(math.sqrt(d1**2 - x**2)*span100percent/100.0)
                x = x0+(x*span100percent/100.0)
                title = "GI #%d: d1 = %f; d2 = %f;" % (i+1,d1,d2)
                svg.append(self._svg_spot(x,y,"blue",title,.8))
            for i in range(len(query_internal_distances)):
                d1 = query_internal_distances[i]
                d2 = sbjct_external_distances[i]
                x = (d1**2 - d2**2 + D**2)/2.0/D
                y = y0-(math.sqrt(d1**2 - x**2)*span100percent/100.0)
                x = x0+(x*span100percent/100.0)
                title = "GI #%d: d1 = %f; d2 = %f;" % (i+1,d1,d2)
                svg.append(self._svg_spot(x,y,"red",title,.8))
            # add titles
            title = self._format_title(query_name,(width-X)/10)
            for i in range(len(list(title))):
                svg.append("<text x=\"%d\" y=\"%d\">%s</text>" % 
                    (X,Y+height-title_zone/2+(i+1)*22,title[i]))
            title = self._format_title(sbjct_name,(width-X)/10)
            for i in range(len(list(title))):
                svg.append("<text x=\"%d\" y=\"%d\">%s</text>" % 
                    (X+width/2+5,Y+height-title_zone/2+(i+1)*22,title[i]))
                
            Y += height
            series += 1
            
        if flg_finish:
            svg.insert(0,"<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % (width+300,height*series,width+200,height*series))
            svg.append("</svg>")
        return "\n".join(svg),width+300,height*series
    
    def _svg_spot(self,X,Y,color,title="",opacity=1.0,r=10):
        svg = ""
        if title:
            svg = "<a title=\"%s\">\n" % title
        svg += ("<circle style=\"stroke:%s;fill:%s;stroke-width:%f;opacity:%f\" cx=\"%f\" cy=\"%f\" r=\"%f\" />" % 
                ("black",color,1.0,opacity,X,Y,r))
        if title:
            svg += "\n</a>"
        return svg
            
    def _svg_cartesian(self,X,Y,width,height,x_indend=30,y_indend=30,flg_finish=False):
        svg = []
        # vertical axis
        svg.append("<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" fill=\"none\" stroke=\"%s\" stroke-width=\"%f\" />" %
            (X+x_indend,Y+20,X+x_indend,Y+height,"black",3.0))
        # vertical arrow
        svg.append(self._svg_arrow(X+30,Y+20,math.pi/2.0,2.0))
        # horizontal axis
        svg.append("<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" fill=\"none\" stroke=\"%s\" stroke-width=\"%f\" />" %
            (X,Y+height-y_indend,X+width-5,Y+height-y_indend,"black",3.0))
        # horizontal arrow
        svg.append(self._svg_arrow(X+width-2,Y+height-30,0,2.0))
        if flg_finish:
            svg.insert(0,"<svg xmlns=\"http://www.w3.org/2000/svg\" viewbox=\"0 0 %d %d\">" % (width,height))
            svg.append("</svg>")
        return "\n".join(svg)

######################################################################################################
class DistanceMatrix(pattern_collection):
    def __init__(self,ptype,first,second=None):
        pattern_collection.__init__(self,ptype,first,second)
        
    def _getColor(self,oFirstWord,oSecondWord):
        n = 4**self.get_wlength()
        p = (oFirstWord.rank+oSecondWord.rank)/2.0/n
        d = 1.5*n*p*(1.0-p)
        if abs(oFirstWord.rank-oSecondWord.rank) > d:
            return "red"
        else:
            return "blue"
        
    def tostring(self):
        return self.distance_table()
    
    def svg(self):
        return

######################################################################################################
class WordSubtraction(pattern_collection):
    def __init__(self,ptype,first,second=None):
        pattern_collection.__init__(self,ptype,first,second)
        
    def _getColor(self,oFirstWord,oSecondWord):
        n = 4**self.get_wlength()
        p = (oFirstWord.rank+oSecondWord.rank)/2.0/n
        d = 1.5*n*p*(1.0-p)
        if abs(oFirstWord.rank-oSecondWord.rank) > d:
            return "red"
        else:
            return "blue"
        
    def tostring(self):
        data = []
        pattern_names = list(self.keys())
        for i in range(len(self)-1):
            first = self[pattern_names[i]]
            first_child = first.getFirstChild()
            for j in range(i+1,len(self)):
                second = self[pattern_names[j]]
                second_child = second.getFirstChild()
                if first_child and second_child:
                    data.extend(self.subtract(first_child,second_child))
                elif first_child:
                    data.extend(self.subtract(first_child,second))
                elif second_child:
                    data.extend(self.subtract(first,second_child))
                else:
                    data.extend(self.subtract(first,second))
                return "\n".join(data)
            
    def subtract(self,first,second):
        data = []
        data.append("%s | %s" % (first.title,second.title))
        data.append("Distance: %s%%" % format_number(first-second,2))
        data.append("\t".join(["Word","dRanks","Rank_1","Freq_1","Rank_2","Freq_2"]))
        subtraction = first | second
        subtraction.sort()
        subtraction.reverse()
        data.extend(["\t".join([s[1],str(s[0]),"\t".join(s[2:])]) for s in subtraction])
        data.append("")
        return data
    
    def svg(self,mode="",width=0,height=0):
        svg = []
        x_indend = 60
        y_indend = 60
        if not width:
            width = 572
        if not height or height != width+y_indend:
            height = width+y_indend
        if len(self):
            pattern_names = list(self.keys())
            for i in range(len(self)-1):
                first = self[pattern_names[i]]
                first_child = first.getFirstChild()
                for j in range(i+1,len(self)):
                    second = self[pattern_names[j]]
                    second_child = second.getFirstChild()
                    if first_child and second_child:
                        svg.append(self.svg_plot(first_child,second_child,width,height,x_indend,y_indend))
                    elif first_child:
                        svg.append(self.svg_plot(first_child,second,width,height,x_indend,y_indend))
                    elif second_child:
                        svg.append(self.svg_plot(first,second_child,width,height,x_indend,y_indend))
                    else:
                        svg.append(self.svg_plot(first,second,width,height,x_indend,y_indend))
                    return "\n".join(["<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % 
                            (width,height,width,height),
                            "\n".join(svg),"</svg>"])
        
    def svg_plot(self,first,second,width=0,height=0,x_indend=30,y_indend=30,flg_finish=False):
        step = float(width-x_indend)/(4**self.get_wlength())
        box_side = step
        if box_side < 2:
            box_side = 2
        svg = []
        svg.append("<path d=\"M%d,%dL%d,%dL%d,%dZ\" style=\"fill:none;stroke:grey;stroke_width:.5\"/>" % 
            (x_indend-1,y_indend,x_indend-1,height-y_indend,width,height-y_indend))
        for oFirstWord in first:
            oSecondWord = second[oFirstWord.word]
            svg.append("<a title=\"%s\"><rect x=\"%d\" y=\"%d\" width=\"%d\" height=\"%d\" style=\"fill:%s\"/></a>" %
                    ("%s, %d, %d" % (oFirstWord.word,oFirstWord.rank,oSecondWord.rank),
                    x_indend+step*(oFirstWord.rank-1),
                    y_indend+step*(oSecondWord.rank-1),
                    box_side,box_side,self._getColor(oFirstWord,oSecondWord)))
        
        svg.append("<text x=\"%f\" y=\"%f\" style=\"text-anchor:middle\" font-size=\"12\">%s</text>" % 
            (x_indend+(width-x_indend)/2.0,y_indend/2.0,("Distance = %f" % float(first-second))+"%"))
        svg.append("<text x=\"%f\" y=\"%f\" style=\"text-anchor:middle\" font-size=\"12\">%s</text>" % 
            (x_indend+(width-x_indend)/2.0,height-y_indend/2.0,first.title))
        svg.append("<text x=\"%f\" y=\"%f\" style=\"text-anchor:middle\" font-size=\"12\" transform=\"rotate(-90)\">%s</text>" % 
            (-(width-x_indend)/2.0 - x_indend,y_indend/2.0,second.title))
        
        if flg_finish:
            svg.insert(0,"<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % (width+300,height*series,width+200,height*series))
            svg.append("</svg>") 
        return "\n".join(svg)
    
######################################################################################################
class Plot3D(pattern_collection):
    def __init__(self,ptype,first,second=None):
        pattern_collection.__init__(self,ptype,first,second)
        self.container = Container(self.do)
        for query_name in self._query_array:
            query = self._query_array[query_name]
            self._addPattern(query_name,"query")
            for child in query.child_patterns:
                self._addPattern(child.title,"query_child")
        for sbjct_name in self._sbjct_array:
            sbjct = self._sbjct_array[sbjct_name]
            self._addPattern(sbjct_name,"sbjct")
            for child in sbjct.child_patterns:
                self._addPattern(child.title,"sbjct_child")
        self.container.set_space(3)
        
    #### METHODS
    def tostring(self):
        return self.container.tostring()
    
    def svg(self,width=0,height=0):
        return self.container.svg(width,height)
        
    def _addPattern(self, oPattern, info=""):
        self.container.add_pattern(oPattern,info)

    def _resolution(self):
        return self.container.get_resolution()
    
    def _dimension(self,ArgList=None):
        return self.container.get_dimension()
    
######################################################################################################
class Search_Report(Report): # The report of a search through OUP_databases
    def __init__(self,dMax=0,stdMax=0):
        Report.__init__(self)
        self.dMax = dMax
        self.stdMax = stdMax
        
        self.oFilter = None
        self.flg_clarified = False
    
    def append(self,seqname):
        self.add(Query_Record(seqname,self.do)) 
    
    def tostring(self):
        headers = []
        reports = []
        for record in self:
            reports.append(record.tostring())
        headers.append("Thresholds: MaxD = %s; Max stDev = %s;" % (format_number(self.dMax,2),
                        format_number(self.stdMax,2)))
        if self.oFilter:
            headers.append("Filter: %s;" % self.oFilter.tostring())
        else:
            headers.append("Filter: None;")
        if self.flg_clarified:
            headers.append("Clarified: True;")
        else:
            headers.append("Clarified: False;")
        headers.append("")
        return "\n".join(headers+reports)
    
    def svg(self,X=5,Y=25,width=800,height=0,flg_finish=True):
        headers = []
        headers.append("Thresholds: MaxD = %s; Max stDev = %s;" % (format_number(self.dMax,2),
                        format_number(self.stdMax,2)))
        if self.oFilter:
            headers.append("Filter: %s;" % self.oFilter.tostring())
        else:
            headers.append("Filter: None;")
        if self.flg_clarified:
            headers.append("Clarified: True;")
        else:
            headers.append("Clarified: False;")
        svg = []
        step = 20
        for line in headers:
            svg.append("<text x=\"%d\" y=\"%d\" style=\"fill:%s;stroke:none;font-weight:%s\">%s</text>" % 
                    (X,Y,"black","bold",line))
            Y += step
        
        for record in self:
            report = record.tostring().split("\n")
            for i in range(len(list(report))):
                line = report[i].split("\t")
                if i == 0:
                    Y += step
                    svg.append("<text x=\"%d\" y=\"%d\" style=\"fill:%s;stroke:none;font-weight:%s\">%s</text>" % 
                        (X,Y,"black","bold",line[0]))
                    svg.append("<text x=\"%d\" y=\"%d\" style=\"fill:%s;stroke:none;font-weight:%s\">%s</text>" % 
                        (self.tab(X,1,line[0]),Y,"red","normal",line[2]))
                    Y += step
                    continue
                if len(list(line)) == 2:
                    svg.append("<text x=\"%d\" y=\"%d\" style=\"fill:%s;stroke:none;font-weight:%s\">%s</text>" % 
                        (self.tab(X,1),Y,"black","normal",line[1]))
                    Y += step
                    continue
                elif len(list(line)) == 4:
                    svg.append("<text x=\"%d\" y=\"%d\" style=\"fill:%s;stroke:none;font-weight:%s\">%s</text>" % 
                        (self.tab(X,2),Y,"red","normal",line[2]))
                    svg.append("<text x=\"%d\" y=\"%d\" style=\"fill:%s;stroke:none;font-weight:%s\">%s</text>" % 
                        (self.tab(X,3,line[2]),Y,"blue","normal",line[3]))
                    Y += step
                    continue
            Y += step
        if flg_finish:
            svg.insert(0,"<svg xmlns=\"http://www.w3.org/2000/svg\" viewbox=\"0 0 %d %d\">" % (width,Y))
            svg.append("</svg>")
        return "\n".join(svg)
    
    def tab(self,x,n,before=""):
        space = 50
        letter = 7
        return space*(n+(x+len(list(before))*letter)/space) 
        

######################################################################################################
class Query_Record(Report): # Sequence record
    def __init__(self,seqname,trigger=None):
        Report.__init__(self)
        self.title = seqname
        self.trigger = trigger
        self.patterns = {}

    def append(self,tbname,stDev,Flag=True):
        self.add(Table_Record(tbname,stDev,Flag,self.do))
    
    def add_pattern(self,ptype,oPattern):
        self.patterns[ptype] = oPattern
    
    def has_pattern(self,ptype):
        return ptype in self.patterns
    
    def get_pattern(self,ptype):
        if ptype not in self.patterns:
            return None
        return self.patterns[ptype]
    
    def get_best_hit(self):
        top_record = self.get_top_record() # top_record = [distance,stdDev,title]
        if not top_record:
            return ""
        a,b,title = top_record
        return "%s [%s] %s" % (format_number(a,2),format_number(b,2),title)
    
    def copy(self):
        oNew = Query_Record(self.title,self.trigger)
        for ptype in self.patterns:
            oNew.add_pattern(ptype,self.patterns[ptype].copy())
        for record in self:
            oNew.add(record)
        return oNew
    
    def tostring(self):
        headers = []
        reports = []
        for record in self:
            reports.append(record.tostring())
        headers.append("%s\t\tBest hit: %s;" % (self.title,self.get_best_hit()))
        return "\n".join(headers+reports)
    
    def SVG(self):
        pass

######################################################################################################
class Table_Record(Report): # Sequence record
    def __init__(self,tbname,stDev,Flag=True,trigger=None):
        Report.__init__(self)
        self.title = tbname
        self.stDev = stDev
        self.Flag = Flag
        self.trigger = trigger
        
    def do(self,command=""):
        if command == "Get stDev":
            return self.stDev

    def append(self,cluster_id):
        self.add(Cluster_Record(cluster_id,self.do))
    
    def copy(self):
        oNew = Table_Record(self.title,self.stDev,self.Flag,self.trigger)
        for record in self:
            oNew.add(record)
        return oNew
    
    def tostring(self):
        headers = []
        reports = []
        for record in self:
            reports.append(record.tostring())
        headers.append("\t%s [%s]" % (self.title,format_number(self.stDev,2)))
        return "\n".join(headers+reports)
    
    def SVG(self):
        pass

######################################################################################################
class Cluster_Record(Report): # Sequence record
    def __init__(self,cluster_id,trigger=None):
        Report.__init__(self)
        self.title = cluster_id
        self.trigger = trigger

    def append(self,seqname,euclidian_distance,corrected_distance):
        self.add(Subject_Record(seqname,euclidian_distance,corrected_distance,self.do))
    
    def copy(self):
        oNew = Cluster_Record(self.title,self.trigger)
        for record in self:
            oNew.add(record)
        return oNew
    
    def get_best_hit(self):
        top_record = self.get_top_record() # top_record = [distance,stdDev,title]
        return top_record
    
    def get_top_record(self):
        hits = []
        for record in self:
            hits.append([record.maxD,record.minD,self.trigger("Get stDev"),record.title])
        if not hits:
            return hits
        if len(hits) > 1:
            for hit in hits:
                hit.insert(0,self.get_hit_score(hit))
            hits.sort()
        return hits[0][1:]
    
    def get_hit_score(self,hit):
        return (hit[0]+hit[1])*hit[2]

    def tostring(self):
        reports = []
        for record in self:
            reports.append(record.tostring())
        return "\n".join(reports)
    
    def SVG(self):
        pass

######################################################################################################
class Subject_Record(Report): # Sequence record
    def __init__(self,seqname,minD,maxD,trigger=None):
        Report.__init__(self)
        self.title = seqname
        self.minD = minD
        self.maxD = maxD
        self.trigger = trigger

    def copy(self):
        oNew = Subject_Record(self.title,self.minD,self.maxD,self.trigger)
        for record in self:
            oNew.add(record)
        return oNew
    
    def tostring(self):
        return "\t\t%s-%s\t%s" % (format_number(self.minD,2),format_number(self.maxD,2),self.title)
    
    def SVG(self):
        pass

###############################################################################
class Container:
    def __init__(self,trigger,outgroups=[]):
        self.trigger = trigger
        self.D12 = 0.0
        # outgroups
        self.outgroups = []
        # collection of reference pattern objects
        self.references = {}
        # coordinates
        self.XYZF = {}
        # hidden patterns
        self.hidden_patterns = []
        # hidden referenses
        self.hidden_referenses = []
        # distances between patterns and the first outgroup
        self.residues = {}
        # matrix of distances between outgroups and all patterns
        self.matrix = {}
        # set outgroups if provided
        if outgroups:
            self.set_outgroups(outgroups)
        
    #### ACCESSIONS
    def do(self,command,ArgList=None):
        if command == "Get distance" and len(ArgList)==2:
            key1,key2 = ArgList
            if key1 in self.matrix and key2 in self.matrix[key1]:
                return self.matrix[key1][key2]
            elif key2 in self.matrix and key1 in self.matrix[key2]:
                return self.matrix[key2][key1]
            else:
                return self.trigger(command,ArgList)
        else:
            return self.trigger(command,ArgList)
        
    def tostring(self):
        axes = ["X","Y","Z","F"]
        dim = self.get_dimension()
        if dim < 4:
            axes = axes[:dim]
        elif dim > 4:
            axes.append("\t"*(dim-4))
        output = ["\t".join(["Sequence","\t".join(axes),"Coverage"])]
        for seqname in self.outgroups:
            output.append(seqname)
            for val in self.XYZF[seqname]:
                output[-1] += "\t%f" % val
            output[-1] += "\t%f" % self.get_coverage(seqname)
        for seqname in self.get_patternList():
            if seqname in self.outgroups:
                continue
            output.append(seqname)
            for val in self.XYZF[seqname]:
                output[-1] += "\t%f" % val
            output[-1] += "\t%f" % self.get_coverage(seqname)
        return "\n".join(output)
    
    def svg(self,width=600,height=400):
        width = 600
        height = 400
        shift_y = 150
        # coordinates
        self.x_A = 5
        self.y_A = 100+shift_y
        self.x_G = 405
        self.y_G = 100+shift_y
        self.x_T = 405
        self.y_T = 350+shift_y
        self.x_C = 5
        self.y_C = 350+shift_y
        self.x_a = 225
        self.y_a = 50+shift_y
        self.x_g = 585
        self.y_g = 50+shift_y
        self.x_t = 585
        self.y_t = 300+shift_y
        self.x_c = 225
        self.y_c = 300+shift_y

        # changeable geometry attributes
        self.ziroX = 0
        self.maxX = 0
        self.ziroY = 0
        self.maxY = 0
        self.middleX = 0

        self.shift_X = 0
        self.shift_Y = 0
        self.shift_Z = 0
        self.scale = 0
        self.min_Z = "Auto"

        self.obliquePegs = None

        self.alpha = (self.y_C - self.y_a)/math.sqrt((self.y_C - self.y_a)*(self.y_C - self.y_a) +
                                                     (self.x_a - self.x_C)*(self.x_a - self.x_C))
        self.beta = math.sqrt(1-self.alpha*self.alpha)
        self.backline_color = "#%02X%02X%02X" % (240,240,240)
        
        svg = self.svg_plot(shift_y)
        svg += self.svg_graph(shift_y)
        return "\n".join(["<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % 
                (width,height+shift_y,width,height+shift_y),svg,"</svg>"])
    
    def svg_plot(self,shift_y=0):
        svg = []
        
        MIN_Y = 0
        D12 = self.get_basicDistance()
            
        # set geometry
        self.ziroX = 180 + self.shift_X - self.shift_Y*self.beta + self.scale
        self.maxX = 500+self.shift_X - self.shift_Y*self.beta - self.scale
        self.ziroY = 150 + self.shift_Y*self.alpha
        self.maxY = self.ziroY + (self.maxX - self.ziroX)*self.beta/self.alpha
        self.middleX = float(self.maxX+self.ziroX)/2-220*float(self.maxY-self.ziroY)/255
        
        svg.append("<path d=\"M%f,%fL%f,%fL%f,%fL%f,%fZ\" style=\"stroke:none;fill:%s;opacity:%f\" />" %
                (self.x_a,self.y_a,self.x_g,self.y_g,self.x_T,self.y_T,self.x_C,self.y_C,'grey',0.8))
        # Back lines
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (435,self.y_t,self.x_c,self.y_c,self.backline_color,1.0))
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (self.x_t,self.y_t,435,self.y_c,"black",1.5))
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (self.x_a,self.y_a,self.x_c,self.y_c,self.backline_color,1.0))
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (self.x_C,self.y_C,self.x_c,self.y_c,self.backline_color,1.0))
        # Front lines
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (self.x_C,self.y_C,self.x_T,self.y_T,"black",1.5))
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (self.x_C,self.y_C,self.x_A,self.y_A,"black",1.5))
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (self.x_A,self.y_A,self.x_a,self.y_a,"black",1.5))
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (self.x_T,self.y_T,self.x_t,self.y_t,"black",1.5))
        svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                (self.x_g,self.y_g,self.x_t,self.y_t,"black",1.5))
        
        # Plot dots
        self.perspective = (self.x_A-self.x_G)/(self.x_a-self.x_g)/2
        self.totalSlope = math.sqrt((self.x_a-self.perspective-self.x_C)**2 + (self.y_C-self.y_a)**2)
        self.YSlope = self.totalSlope*(self.maxY-self.ziroY)/(self.y_C-self.y_a)
        self.midX = self.ziroX + (self.maxX + self.ziroX)/2

        # Add marking values
        text_Y = 370
        xyzf = list(self.XYZF.values())
        x1 = self.ziroX + min(xyzf)[0]*(self.maxX-self.ziroX)/D12
        x1 = x1 - (((self.midX - x1)/(self.midX - self.ziroX))*
                   (((text_Y + self.ziroY)/self.alpha)/self.totalSlope)*self.perspective +
                   (self.x_a-self.x_C)*(text_Y-self.ziroY)/self.alpha/self.totalSlope)
        x2 = self.ziroX + (max(xyzf)[0]-min(xyzf)[0])*(self.maxX-self.ziroX)/2.0/D12
        x2 = x2 - (((self.midX - x2)/(self.midX - self.ziroX))*
                   (((text_Y + self.ziroY)/self.alpha)/self.totalSlope)*self.perspective +
                   (self.x_a-self.x_C)*(text_Y-self.ziroY)/self.alpha/self.totalSlope)
        x3 = self.ziroX + max(xyzf)[0]*(self.maxX-self.ziroX)/D12
        x3 = x3 + (((self.midX - x2)/(self.midX - self.ziroX))*
                   (((text_Y + self.ziroY)/self.alpha)/self.totalSlope)*self.perspective -
                   (self.x_a-self.x_C)*(text_Y-self.ziroY)/self.alpha/self.totalSlope)
        svg.append("<text x=\"%f\" y=\"%f\" text-anchor=\"left\">%s</text>" % (x1,text_Y+shift_y,"0"))
        svg.append("<text x=\"%f\" y=\"%f\" text-anchor=\"left\">%s</text>" % (x2,text_Y+shift_y,str(int(D12/2.0))+"%"))
        svg.append("<text x=\"%f\" y=\"%f\" text-anchor=\"left\">%s</text>" % (x3,text_Y+shift_y,str(int(D12))+"%"))

        # Show statistics
        svg.append("<text x=\"%f\" y=\"%f\" text-anchor=\"left\">%s</text>" % 
                (50,45,"%u dimensional space" % self.get_dimension()))
        svg.append("<text x=\"%f\" y=\"%f\" text-anchor=\"left\">%s</text>" % 
                (50,60,"resolution %g%%" % (100.0*self.get_resolution())))
        return "\n".join(svg)
    
    def svg_graph(self,y_shift):
        svg = []
        for pattern in self.get_patternList():
            svg.append(self.drawElement(pattern,self.get_basicDistance(),y_shift))
        return "\n".join(svg)
    
    def drawElement(self,pattern,D12,y_shift):
        svg = []
        # basic x
        x = self.ziroX + self.XYZF[pattern][0]*(self.maxX-self.ziroX)/D12
        # basic y
        if len(self.XYZF[pattern]) > 1:
            y = self.ziroY + self.XYZF[pattern][1]*self.YSlope*self.alpha/D12
        else:
            y = self.ziroY
        # correction for perspective
        if x < self.midX:
            x = x - (((self.midX - x)/(self.midX - self.ziroX))*(((y + self.ziroY)/self.alpha)/self.totalSlope)*self.perspective +
                     (self.x_a-self.x_C)*(y-self.ziroY)/self.alpha/self.totalSlope)
        else:
            x = x + (((self.midX - x)/(self.midX - self.ziroX))*(((y + self.ziroY)/self.alpha)/self.totalSlope)*self.perspective -
                     (self.x_a-self.x_C)*(y-self.ziroY)/self.alpha/self.totalSlope)
        # adjust Z coefficients
        Z = []
        xyzf = list(self.XYZF.values())
        for item in xyzf:
            if len(list(item)) > 2:
                Z.append(item[2]/D12)
            else:
                Z.append(0.0)
        if self.min_Z == 'Auto':
            if self.shift_Z:
                for item in xyzf:
                    item[2] -= self.shift_Z
                sef.shift_Z = 0
        elif self.min_Z == 'Min':
            self.shift_Z = -min(Z)
            for item in xyzf:
                item[2] += self.shift_Z
        elif self.min_Z == 'Middle':
            minVal = min(Z)
            maxVal = max(Z)
            self.shift_Z = -(minVal+maxVal)/2.0
            for item in xyzf:
                item[2] += self.shift_Z
        elif self.min_Z == 'Max':
            self.shift_Z = -max(Z)
            for item in xyzf:
                item[2] += self.shift_Z
        else:
            for item in xyzf:
                item[2] += self.shift_Z
        # correction for z
        if self.obliquePegs:
            if len(self.XYZF[pattern]) > 2:
                peg = (self.maxX - self.ziroX)*self.XYZF[pattern][2]*(1-self.alpha*self.alpha)/D12
            else:
                peg = 0
            xz = x - peg*self.beta
            yz = y - peg*self.alpha
        else:
            if len(self.XYZF[pattern]) > 2:
                peg = (self.maxX - self.ziroX)*self.XYZF[pattern][2]/D12
            else:
                peg = 0
            xz = x
            yz = y - peg
        # draw dot
        svg.append("<a title=\"%s\">" % pattern)
        if not self.references[pattern] or self.references[pattern]=="query":
            svg.append("<rect x=\"%f\" y=\"%f\" width=\"%f\" height=\"%f\" style=\"stroke:none;fill:%s\" /></a>" %
                (xz-4,yz-4+y_shift,7,7,"red"))
            svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                    (x,y+y_shift,xz,yz+y_shift,"red",.5))
        if self.references[pattern]=="sbjct":
            svg.append("<rect x=\"%f\" y=\"%f\" width=\"%f\" height=\"%f\" style=\"stroke:none;fill:%s\" /></a>" %
                (xz-4,yz-4+y_shift,7,7,"blue"))
            svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                    (x,y+y_shift,xz,yz+y_shift,"blue",.5))
        if self.references[pattern]=="query_child":
            svg.append("<circle cx=\"%f\" cy=\"%f\" r=\"%f\" style=\"stroke:none;fill:%s\" /></a>" % 
                (xz,yz+y_shift,3,"red"))
            svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                    (x,y+y_shift,xz,yz+y_shift,"red",.5))
        if self.references[pattern]=="sbjct_child":
            svg.append("<circle cx=\"%f\" cy=\"%f\" r=\"%f\" style=\"stroke:none;fill:%s\" /></a>" % 
                (xz,yz+y_shift,3,"blue"))
            svg.append("<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:%s;stroke-width:%f\" />" %
                    (x,y+y_shift,xz,yz+y_shift,"blue",.5))
        return "\n".join(svg)
    
    def add_pattern(self,pattern,info=""):
        if pattern in self.XYZF:
            return
        # Fill in the matrix
        for seqname in self.outgroups:
            if pattern not in self.matrix[seqname]:
                self.matrix[seqname][pattern] = self.trigger("Get distance",[pattern,seqname])
        if pattern not in self.XYZF:
            self.XYZF[pattern] = []
        if self.outgroups and pattern not in self.residues:
            self.residues[pattern] = self.matrix[self.outgroups[0]][pattern]**2
        if not self.outgroups:
            self.residues[pattern] = 0.0
        for dim in range(1,len(self.outgroups)):
            if len(self.XYZF[pattern]) >= dim:
                continue
            self.next_coordinate(pattern,dim)
        self.references[pattern] = info
    
    def add_outgroup(self,pattern):
        if pattern not in self.XYZF:
            self.XYZF[pattern] = []
        # Set residues table
        if len(self.outgroups)==0 and self.residues:
            for seqname in self.residues:
                if seqname == pattern:
                    self.residues[seqname] = 0.0
                else:
                    self.residues[seqname] = self.trigger("Get distance",[pattern,seqname])**2
        if self.outgroups:
            basic_distance = self.trigger("Get distance",[self.outgroups[0],pattern])
        else:
            basic_distance = 0.0
        if len(self.outgroups)==1:
            self.set_basicDistance(basic_distance)
        if pattern not in self.residues:
            self.residues[pattern] = basic_distance**2
        # Fill in the matrix
        self.matrix[pattern] = {}
        for seqname in self.XYZF:
            if pattern == seqname:
                self.matrix[pattern][seqname] = 0.0
            else:
                distance = self.trigger("Get distance",[pattern,seqname])
                self.matrix[pattern][seqname] = distance
        for seqname in self.outgroups:
            if pattern not in self.matrix[seqname]:
                if pattern == seqname:
                    self.matrix[seqname][pattern] = 0
                else:
                    self.matrix[seqname][pattern] = self.trigger("Get distance",[seqname,pattern])
        # Fill in coordinates up to existing dimensions
        if len(self.outgroups) > 1 and len(self.XYZF[pattern]) < len(self.outgroups)-1:
            for dim in range(1,len(self.outgroups)):
                self.next_coordinate(pattern,dim)
        # add new outgroup
        self.outgroups.append(pattern)
        # calculate  for all patterns next coordinates for the dimension created by new outgroup
        dim = len(self.outgroups)-1
        if dim > 0:
            for seqname in self.XYZF:
                self.next_coordinate(seqname,dim)
   
    def get_visible(self):
        visible = []
        for pattern in self.XYZF:
            if pattern not in self.hidden_patterns:
                visible.append(pattern)
        return visible
    
    def hide(self,pattern):
        if pattern not in self.XYZF:
            return
        self.hidden_patterns.append(pattern)
        
    def disclose(self,pattern):
        if pattern in self.hidden_patterns:
            self.hidden_patterns.remove(pattern)
        
    def get_outgroup(self,index):
        if index >= 0 and index < len(self.outgroups):
            return self.outgroups[index]
        else:
            return None
        
    def set_outgroups(self,outgroups):
        dim = len(outgroups)
        if dim > 5:
            dim = 5
        for index in range(dim):
            if index >= len(self.outgroups):
                break
            if self.outgroups[index] != outgroups[index]:
                self.delslice(index)
                break
        if index==dim-1:
            return
        for k in range(index,dim):
            self.add_outgroup(outgroups[k])
        
    def keys(self):
        outgroups = []
        outgroups.extend(self.outgroups)
        return outgroups
    
    def values(self):
        return list(self.residues.values())
    
    def index(self,pattern):
        if pattern in self.outgroups:
            return self.outgroups.index(pattern)
        else:
            return None
    
    def delete_pattern(self,pattern_name):
        self.disclose(pattern_name)
        if pattern_name in self.outgroups:
            alert("Pattern " + pattern_name + " is an outgroup and cannot be removed")
            return
        if pattern_name in self.XYZF:
            del self.XYZF[pattern_name]
        if pattern_name in self.residues:
            del self.residues[pattern_name]
        if pattern_name in self.references:
            del self.references[pattern_name]
        return 1
        
    def rename_pattern(self,oldname,newname):
        if oldname in self.outgroups:
            i = self.outgroups.index(oldname)
            self.outgroups[i] = newname
        if oldname in self.references:
            self.references[newname] = self.references[oldname]
            del self.references[oldname]
        if oldname in self.residues:
            self.residues[newname] = self.residues[oldname]
            del self.residues[oldname]
        if oldname in self.XYZF:
            self.XYZF[newname] = []
            self.XYZF[newname].extend(self.XYZF[oldname])
            del self.XYZF[oldname]
        if oldname in self.hidden_patterns:
            i = self.hidden_patterns.index(oldname)
            self.hidden_patterns[i] = newname
        if oldname in self.hidden_referenses:
            i = self.hidden_referenses.index(oldname)
            self.hidden_referenses[i] = newname
        for seqname in self.matrix:
            if oldname in self.matrix[seqname]:
                self.matrix[seqname][newname] = self.matrix[seqname][oldname]
                del self.matrix[seqname][oldname]
            if seqname == oldname:
                self.matrix[newname] = {}
                self.matrix[newname].update(self.matrix[oldname])
                del self.matrix[oldname]
    
    def get_coverage(self,pattern):
        if not self.outgroups:
            return 0.0
        if pattern in self.residues and pattern in self.matrix[self.outgroups[0]]:
            if self.matrix[self.outgroups[0]][pattern] <= 1 or self.residues[pattern]==0:
                return 1.0
            basis = abs(self.matrix[self.outgroups[0]][pattern]-50)
            basis = math.exp(basis-35)
            sign = abs(self.residues[pattern])/self.residues[pattern]
            coverage = (self.matrix[self.outgroups[0]][pattern]-sign*math.sqrt(abs(self.residues[pattern]))+basis)/(self.matrix[self.outgroups[0]][pattern]+basis)
            return coverage
        else:
            return None
        
    def get_resolution(self):
        resolution = 0.0
        if not self.XYZF:
            return resolution
        for pattern in self.XYZF:
            if self.matrix[self.outgroups[0]][pattern] == 0:
                resolution += 1.0
            else:
                resolution += 1.0 - abs(1.0-self.get_coverage(pattern))
                #resolution += 1.0 - math.sqrt(abs(self.residues[pattern]))/self.matrix[self.outgroups[0]][pattern]
        resolution /= float(len(self.XYZF))
        return resolution
    
    def calculate_coverage(self,pattern):
        if self.matrix[self.outgroups[0]][pattern] == 0:
            return 0.0
        basis = abs(self.matrix[self.outgroups[0]][pattern]-50)
        return (self.matrix[self.outgroups[0]][pattern]-math.sqrt(abs(self.residues[pattern]))+basis)/(self.matrix[self.outgroups[0]][pattern]+basis)
    
    def set_coordinates(self,seqname,values):
        if seqname not in self.XYZF:
            print(("auxiliaries:1150, Error setting the coordinate for the pattern " + seqname))
            return None
        self.XYZF[seqname] = []
        self.XYZF[seqname].extend(values)
        
    def get_RelativeCoordinates(self,pattern):
        if pattern not in self.XYZF:
            return None
        coordinates = []
        coordinates.extend(self.XYZF[pattern])
        if self.D12:
            D12 = float(self.D12)
        else:
            D12 = 1.0
        for i in range(len(coordinates)):
            coordinates[i] /= D12
        return coordinates
    
    def get_AbsoluteCoordinates(self,pattern):
        if pattern not in self.XYZF:
            return None
        coordinates = []
        coordinates.extend(self.XYZF[pattern])
        return coordinates
    
    def get_coordinates(self):
        xyzf = {}
        xyzf.update(self.XYZF)
        return xyzf
    
    def get_distance(self,first,second):
        if first == second:
            return 0.0
        elif first in self.matrix and second in self.matrix[first]:
            return self.matrix[first][second]
        elif second in self.matrix and first in self.matrix[second]:
            return self.matrix[second][first]
        elif first in self.references and second in self.references:
            return self.references[first] - self.references[second]
        else:
            return None
    
    def get_patternList(self):
        return list(self.XYZF.keys())
    
    def set_basicDistance(self,distance=None):
        if distance:
            self.D12 = distance
            return
        if len(self.outgroups) < 2:
            return
        self.D12 = self.trigger("Get distance",[self.outgroups[0],self.outgroups[1]])
        
    def get_basicDistance(self):
        return self.D12
    
    def get_dimension(self):
        if self.outgroups:
            return len(self.outgroups)-1
        else:
            return 0
    
    def get_dataset(self):
        dataset = {"All":[],"References":[],"Main":[]}
        for pattern in self.XYZF:
            if pattern not in self.hidden_patterns:
                dataset["All"].append(pattern)
        dataset["References"].extend(list(self.references.keys()))
        dataset["Main"].extend(self.outgroups)
        return dataset
        
    def delslice(self,dim):
        if dim > len(self.outgroups):
            return
        for i in range(len(self.outgroups)-1,dim-1,-1):
            del self.outgroups[i]
        for pattern in self.XYZF:
            if dim:
                self.XYZF[pattern] = self.XYZF[pattern][:dim-1]
            else:
                self.XYZF[pattern] = []
            if self.outgroups:
                self.residues[pattern] = self.matrix[self.outgroups[0]][pattern]**2
            for i in range(len(self.XYZF[pattern])):
                self.residues[pattern] -= self.XYZF[pattern][i]**2
        items_todelete = []
        for seqname in self.matrix:
            if seqname not in self.outgroups:
                items_todelete.append(seqname)
        for seqname in items_todelete:
            del self.matrix[seqname]
    
    def size(self):
        return len(self.XYZF)
    
    def wlength(self):
        if self.outgroups:
            try:
                return self.references[self.outgroups[0]].getWordLength()
            except:
                return None
        return None

    def copy(self):
        container = Container()
        container.outgroups.extend(self.outgroups)
        container.XYZF.update(self.XYZF)
        container.residues.update(self.residues)
        container.hidden_patterns.extend(self.hidden_patterns)
        container.hidden_referenses.extend(self.hidden_referenses)
        container.D12 = self.D12
        return container
    
    def get_dictionary(self):
        container = {"outgroups":[],"coordinates":{},"referenses":{},"residues":{},
            "matrix":{},"hidden patterns":[],"hidden referenses":[]}
        container["outgroups"].extend(self.outgroups)
        container["coordinates"].update(self.XYZF)
        container["referenses"].update(self.references)
        container["residues"].update(self.residues)
        container["matrix"].update(self.matrix)
        container["hidden patterns"].extend(self.hidden_patterns)
        container["hidden referenses"].extend(self.hidden_referenses)
        return container
    
    def load(self,container):
        self.outgroups = []
        self.outgroups.extend(container["outgroups"])
        self.matrix = {}
        self.matrix.update(container["matrix"])
        if len(self.outgroups) > 1:
            self.set_basicDistance(self.matrix[self.outgroups[0]][self.outgroups[1]])
        self.XYZF = {}
        self.XYZF.update(container["coordinates"])
        self.references = {}
        self.references.update(container["referenses"])
        self.residues = {}
        self.residues.update(container["residues"])
        self.hidden_patterns = []
        self.hidden_patterns.extend(container["hidden patterns"])
        self.hidden_referenses = []
        self.hidden_referenses.extend(container["hidden referenses"])
    
    def update(self,other):
        self.outgroups = []
        self.outgroups.extend(other.outgroups)
        self.matrix = {}
        self.matrix.update(other.matrix)
        self.references = {}
        self.references.update(other.references)
        self.XYZF = {}
        self.XYZF.update(other.XYZF)
        self.residues = {}
        self.residues.update(other.residues)
        self.hidden_patterns = []
        self.hidden_patterns.extend(other.hidden_patterns)
        self.hidden_referenses = []
        self.hidden_referenses.extend(other.hidden_referenses)
        self.D12 = other.D12
    
    #### METHODS
    def next_coordinate(self,pattern,dim=None):
        # First the confidence must be set for this pattern
        if dim == None:
            dim = len(self.outgroups)-1
        x = self.get_coordinate(pattern,dim)
        self.XYZF[pattern].append(x)
        self.residues[pattern] -= x**2
        if abs(self.residues[pattern]) < 1E-010:
            self.residues[pattern] = 0.0
        
    def get_coordinate(self,pattern,dim):
        if dim == 0:
            return
        first_pattern = self.outgroups[0]
        last_pattern = self.outgroups[dim]
        try:
            a = self.matrix[first_pattern][pattern]**2
        except:
            print(("auxiliaries:1361",first_pattern,pattern,self.matrix[first_pattern][pattern]))
            5/0
        try:
            b = self.matrix[last_pattern][pattern]**2
        except:
            print(("auxiliaries:1366",last_pattern,pattern,self.matrix[last_pattern][pattern]))
            5/0
        try:
            c = self.matrix[first_pattern][last_pattern]**2
        except:
            print(("auxiliaries:1371",first_pattern,last_pattern,self.matrix[first_pattern][last_pattern]))
            5/0
        if c==0:
            return self.XYZF[first_pattern][dim-1]
        d1 = d2 = 0
        for i in range(dim-1):
            c -= self.XYZF[last_pattern][i]**2
            a += (self.XYZF[pattern][i] - self.XYZF[last_pattern][i])**2 - self.XYZF[pattern][i]**2
        if a == 0 or c == 0:
            return 0.0
        if abs(c) < 1E-010:
            return 0.0
        try:
            x = (a - b + c)/math.sqrt(c)/2
        except:
            print(("auxiliaries:1393",c,last_pattern))
            5/0
        if abs(x) < 1E-010:
            return 0.0
        else:
            return x
        
    def calculate_distance(self,first,second):
        if not self.outgroups:
            return None
        if first in self.XYZF:
            first_coordinates = self.XYZF[first]
            first_confidence = self.residues[first]
        else:
            result = self.calculate_coordinates(first)
            if not result:
                return None
            first_coordinates,first_confidence = result
        if second in self.XYZF:
            second_coordinates = self.XYZF[second]
            second_confidence = self.residues[second]
        else:
            result = self.calculate_coordinates(second)
            if not result:
                return None
            second_coordinates,second_confidence = result
        distance = get_distanceByCoordinates(first_coordinates,second_coordinates)
        if distance:
            return distance*self.D12
        else:
            return distance
        
    def set_space(self,dNum=0,maxResolution=0.75):
        resolution = 0
        for pattern in self.references:
            self.add_pattern(pattern)
        # if outgroups were not set
        if not self.outgroups:
            oNode = Node(0,self.do)
            counter = 0
            for categories in (("query","sbjct"),("query_child","sbjct_child"),("")):
                if dNum and counter > dNum+1:
                    break
                for pattern in self.get_patternList():
                    if self.references[pattern] not in categories:
                        continue
                    oNode.add(pattern)
                    counter += 1
            outgroups = oNode.getOutermosts()
            if outgroups:
                self.set_outgroups(outgroups)
            del oNode
        if len(self.outgroups)==1:
            self.XYZF[self.outgroups[0]] = [0.0]
        if dNum and dNum <= self.get_dimension():
            self.delslice(dNum+1)
            return
        dim = 1
        while dim > 0:
            if len(self.XYZF) < dim + 1:
                break
            # set next outgroup pattern
            if len(self.outgroups) < dim + 1:
                values = list(self.values())
                maxVal = max(values)
                for seqname in self.XYZF:
                    if self.residues[seqname] == maxVal:
                        self.add_outgroup(seqname)
                        break
            dim += 1
            new_resolution = self.get_resolution()
            if not dNum and new_resolution <= resolution:
                self.delslice(dim-1)
                break
            else:
                resolution = new_resolution
            if not dNum and dim > 2 and resolution >= maxResolution:
                break
            if dNum and dim == dNum+1:
                break
            elif dim > 20:
                break
            else:
                pass
                
    def get_StdDeviation(self,seqname,oPattern=None):
        if seqname not in self.XYZF:
            return
        if not self.references:
            return 0.0
        sum = 0
        for pattern in self.references:
            if pattern in self.hidden_referenses:
                continue
            if oPattern:
                dist = self.references[pattern]-oPattern
            else:
                dist = self.get_distance(seqname,pattern)
            if dist==None:
                dist = self.trigger("Get distance",[seqname,pattern])
            sum += (dist - get_distanceByCoordinates(self.XYZF[seqname],self.XYZF[pattern],
                        self.get_coverage(seqname),self.get_coverage(pattern)))**2
        return math.sqrt(sum)/len(self.references)

    def report(self):
        output = []
        for seqname in self.XYZF:
            for pattern in self.XYZF:
                if seqname == pattern:
                    continue
                result = [str(math.log(math.sqrt(2.0/(self.get_coverage(seqname)**2+self.get_coverage(pattern)**2)))),
                    str((self.trigger("Get distance",[seqname,pattern])-get_EuclidianDistance(self.XYZF[seqname],self.XYZF[pattern]))/self.D12)]
                output.append("\t".join(result))
        return "\n".join(output)
        
    # Container emulating methods    
    def __call__(self,pattern,basic_distance=None):
        return self.calculate_coordinates(pattern,basic_distance)
        
    def __len__(self):
        return len(self.outgroups)
    
    def __iter__(self):
        return list(self.keys())
    
    def __getitem__(self,index):
        if type(index) == type(1):
            return self.outgroups[index]
        else:
            return self.residues[index]
    
    def __setitem__(self,key,value):
        #self.confidence[key] = value
        pass

###################################################################################################################
class Node:
    def __init__(self, level=0, trigger=None):
        
        # CONSTANTS and flags
        # Value of this attribute - the method 'do' of the parent node ia
        # assigned by the method 'insert' only.
        self.trigger = trigger
        # Mark a node that requires updating of own neighbours and neighbours
        # in all child nodes.
        # This flag may be set by 'setOutermosts', 'fixOutermosts', 'insert', 'resetUpdateNeed' and 'copy' methods
        self.needUpdate = 0

        # ATTRIBUTES
        # collection of the subordinate nodes
        # new nodes are added by the 'insert' method only.
        # the container may is empty only in leaf nodes.
        self.container = []
        # List of names of outermost elements of the node.
        # It may have 2 elements, 1 element if the node has only 1 child and
        # in leaf nodes. This attribute is logically linked to the 'distance'
        # and 'neighbours' attributes
        self.outermosts = []
        # Temporary variable to keep previous outermost elements after new one was set
        # until updating of the node. It is set by 'fixOutermosts', copied by 'copy', and reset by 'update' methods
        self.previousOutermosts = None
        # List of neighbours in the format [[value,internalOutermost,externalOutermost],...].
        # It can not be empty but it is None in the top level node.
        # It may be changed by 'setNeighbourhood' and 'copy' methods.
        self.neighbours = None
        # The distance between outermost elements. 
        # It may be changed by 'setOutermosts','fixOutermosts' and 'copy' methods.
        self.distance = 0
        # The level in the node hierarchy.
        # If not passed in arguments, it is assigned by the 'insert' or 'copy'
        # It is passed in arguments only in the method 'split'
        self.level = level
        # The index of a child node in the collection of the parent node.
        # Only methods 'insert' and 'copy' may change it.
        self.index = 0
        # The length of the collection, corresponds to len(self.container).
        # May be changed by the 'insert' and 'split' methods only.
        self.length = 0
        # The name of the element. Is assigned only in leaf nodes by the method 'insert'
        self.name = ""
        self.thresholds = [50,30,20,10]
        
    #### METHODS
    
    # This method is passed from the parent bode to child nodes to establich the
    # event communication between them
    def do(self,mode,ArgList=None):
        if mode == "Split":
            self.split(ArgList)
        elif mode == "Feedback":
            self.update()
        elif mode == "Get outermost":
            return self.getOutermosts()
        elif mode == "Path":
            if self.isTopLevelNode():
                return self.getPath(ArgList)
            else:
                return self.trigger("Path",ArgList)
        elif mode == "Size":
            if self.isTopLevelNode():
                return self.size()
            else:
                return self.trigger("Size")
        elif mode == "SnapShot":
            if self.isTopLevelNode():
                self.takeSnapShot(ArgList)
            else:
                self.trigger("SnapShot",ArgList)
        elif mode == "Get thresholds":
            thresholds = self.trigger(mode,ArgList)
            if not thresholds:
                return self.thresholds
        else:
            return self.trigger(mode,ArgList)
    
    # This method calls 'split' method in the parent node if conditions require 
    # fixing the attributes of the parent node, and it triggers updating of the 
    # parent node until the chain of callsreachs the top level node.
    # This method is callable by 'add' method when a new leaf element was added 
    # and then by 'do' method for the parent nodes until the top level node will be reached. 
    def update(self):
        thresholds = self.trigger("Get thresholds")
        if not thresholds:
            thresholds = self.thresholds
        # check if child outermosts are bigger than own
        for node in self.container:
            if self.distance < node.distance:
                self.setOutermosts(node.getOutermosts(),node.distance)
        # check if neighbours of this node need update
        if self.needUpdate:
            self.setNeighbourhood()
        if not self.isTopLevelNode():
            # check if the node needs splitting
            if self.length > 1 and int(self.distance) > thresholds[self.level]:
                self.trigger("Split",self.index)
            self.previousOutermosts = None
            self.trigger("Feedback")
        
    # INPUT METHODS
    
    # Methods to manipulate with the nodes
    
    # This is only one publick method available to change the database.
    # Add a new element to the database
    # evoke the method add for subordinate nodes
    # or the method insert to add a new node on the same level.
    # This method may change the outermost elements of the
    # current node that requires updating of the neighbours of this
    # node and all subordinate nodes
    def add(self,name,topNeighbour=None):
        thresholds = self.trigger("Get thresholds")
        if not thresholds:
            thresholds = self.thresholds
        # check if the name is unique
        if self.level == 0:
            if self.getPath(name):
                #print "The entry " + name + " already exists!"
                return
        if self.length:
            # Find the closest and the most distant elements in the child nodes
            neighbour = [100,""]
            outer = [0,""]
            for j in range(self.length):
                local_distances = self.container[j].getLocalDistances(name)
                if min(local_distances) < neighbour:
                    neighbour = min(local_distances)
                if max(local_distances) > outer:
                    outer = max(local_distances)

            # check if the outermost element need fixing
            if outer[0] > self.distance:
                self.fixOutermosts(name,outer)

            # if the patterns are the same but names are different
            if self.distance == 0 and outer[0] == 0 and name != outer[1]:
                self.setOutermosts([name,outer[1]])
                
            # Set the outermost nodes
            if outer[0] > self.distance:
                # change the outermost elements
                self.fixOutermosts(name,outer)
            
            # Insert a new leaf node
            if self.isLeavesContainer():
                self.insert(Node(),name)
                # update itself and all parent nodes
                self.update()
                return
            # Add new element to the closest child node    
            elif neighbour[0] < thresholds[self.level+1]:
                index = self.indexOf(neighbour[1])
                self.container[index].add(name,neighbour)
            # Insert a new empty node into the container and add
            # the new element to this node
            else:
                self.insert(Node())
                self.container[-1].add(name,neighbour)                
            return
        else:
            if self.isLeavesContainer():
                self.insert(Node(),name)
                self.setOutermosts([name])
                # update itself and all parent nodes
                self.update()
            else:
                self.setOutermosts([name])
                self.insert(Node())
                self.container[-1].add(name)
            return
    
    # This method inserts a new node into collection of nodes of the current node.
    # The argument 'name' is submitted for inserting leaf nodes only.
    # Inserting of a node itself can not influence the set of outermost and neighbour
    # elements of the current node, but the index of the current node and neighbours 
    # of the inserted node must be updated!!!
    # The link between the parent and child nodes must be established as well.
    # This method is used by the methods 'add', 'split', 'propagateNodes',
    # 'updateOutermosts' and 'copy'.
    def insert(self,node,name=None):
        node.trigger = self.do
        node.level = self.level+1
        if node.isLeafNode():
            if name:
                node.name = name
                node.setOutermosts([name])
            else:
                print(("Error in insert",node.level,5/0))
        self.container.append(node)
        node.index = len(self.container)-1
        self.length += 1
        # inserted node needs update
        node.needUpdate = 1
        # container of leaf nodes needs update as well, other wise its child nodes
        # never would be updated
        if self.isLeavesContainer():
            self.needUpdate = 1
        return node
    
    # This method is called from one of child nodes that passes its own index in
    # the parent node's collection. Method 'update' triggers the method 'split' in
    # the parent node.
    # The method splits the child node into two nodes and insert them into collection.
    # Following things must be kept in memory:
    # 1) number of nodes in hierarchi shouldn't be changed. Method first inserts two
    #   empty nodes and than inserts into these nodes the new split nodes;
    # 2) splitting a child node doesn't influence the set of outermost and neighbour
    #   elements of the current node;
    # 3) while the method 'insert' takes care for updating the neighbours of the inserted
    #   nodes, inserting of nodes into new empty node is prone of errors. The method must 
    #   ensure that sets of outermost and neighbour elements in all split nodes are
    #   fixed properly.
    def split(self, index):
        # temporary container for child nodes
        tmp_container = []
        # fill the container with all child nodes but one that is to be split
        for item in self.container:
            if item.index != index:
                tmp_container.append(item.copy())
            else:
                node_to_split = item.copy()
            del item
        # compare own outermosts with node_to_split's outermosts
        flg_ChildNodeNeedUpdate = None
        if self.distance < node_to_split.distance:
            self.setOutermosts(node_to_split.getOutermosts(),node_to_split.distance)
            flg_ChildNodeNeedUpdate = 1
        # empty the permanent container
        self.container = []
        self.length = 0
        # insert child nodes back to the permanent container
        for item in tmp_container:
            self.insert(item)
            if not flg_ChildNodeNeedUpdate:
                # no need to update these nodes
                item.resetUpdateNeed()
        # Create and insert 2 empty nodes and assign their levels
        first_node = self.insert(Node(self.level+1))
        second_node = self.insert(Node(self.level+1))
        self.needUpdate = 1
        # copy child nodes from the node_to_split into 2 new empty nodes
        self.propagateNodes(node_to_split,first_node,second_node)
        # set neighbours
        self.setNeighbourhood()
        return first_node,second_node
    
    # This method is part of the process started by 'split' method. 
    # Any other methods should not use it.
    # It throws the error when:
    # 1) the list of outermosts in the node_to_split contains only 1 element
    # 2) child node doesn't relate to any of outermost elements
    def propagateNodes(self,node_to_split,first_node,second_node):
        # Find indices of the outermost elements
        outermosts = node_to_split.getOutermosts()
        if len(outermosts) < 2:
            print(("Error in propagating the nodes", node_to_split.toString(), 5/0))
        i = node_to_split.indexOf(outermosts[0])
        j = node_to_split.indexOf(outermosts[1])
        if i == None or j == None:
            print(("Error propagating nodes", i,j))
        # Both outermost elements came from the same child node. It means that this node
        # must be split as well
        if i == j:
            print((4778, node_to_split.trigger))
            node1,node2 = node_to_split.split(i)
            first_node.insert(node1,node1.name)
            second_node.insert(node2,node2.name)
        else:
            # Insert outermost nodes to the corresponding empty nodes of higher level
            first_node.insert(node_to_split.container[i].copy(),node_to_split.container[i].name)
            second_node.insert(node_to_split.container[j].copy(),node_to_split.container[j].name)
            # Distribute other children among outermost nodes
            if node_to_split.length > 2:
                for k in range(node_to_split.length):
                    if k not in (i,j):
                        if node_to_split.container[k].getNeighbour() == outermosts[0]:
                            first_node.insert(node_to_split.container[k].copy(),
                            node_to_split.container[k].name)
                        elif node_to_split.container[k].getNeighbour() == outermosts[1]:
                            second_node.insert(node_to_split.container[k].copy(),
                            node_to_split.container[k].name)
                        else:
                            print(("Error in propagating nodes",5/0))
        if outermosts[0] in first_node.getLeaves():
            first_node.updateOutermosts(outermosts[0],node_to_split.previousOutermosts)
            second_node.updateOutermosts(outermosts[1],node_to_split.previousOutermosts)
        else:
            first_node.updateOutermosts(outermosts[1],node_to_split.previousOutermosts)
            second_node.updateOutermosts(outermosts[0],node_to_split.previousOutermosts)
            
        del node_to_split
        
    def rename_pattern(self,oldname,newname):
        if self.name == oldname:
            self.name = newname
        for i in range(len(self.outermosts)):
            if self.outermosts[i] == oldname:
                self.outermosts[i] = newname
        if self.neighbours:
            for i in range(len(self.neighbours)):
                if oldname in self.neighbours[i]:
                    j = self.neighbours[i].index(oldname)
                    self.neighbours[i][j] = newname
        for node in self.container:
            node.rename_pattern(oldname,newname)
        if self.previousOutermosts == oldname:
            self.previousOutermosts = newname
       
    # Methods to manupulate with outermost elements
    # Each node keeps names of one or two outermost elements:
    # self.outermosts = [first,second - may be omitted]
    # self.distance is the distance between first and second outermosts
    # If there is only one outermost element, self.distance = 0
    # self.outermosts = None only in the node of top level
    # In the leaf nodes self.outermosts always contains only one
    # element - the host of the node
    
    # This method is callable by 'propagateNodes' method only.
    # It has to fix names of outermost elements in new empty nodes
    # after they were inserted and populated with split child nodes.
    # Logic of splitting is: split-propagate-update.
    # If the node is empty, this method throws an error.
    def updateOutermosts(self,first_outermost,previousOutermosts):
        size = self.size()
        if size == 1:
            self.setOutermosts([first_outermost])
        elif size > 1:
            # Select second outermost element
            distantNeighbour = self.findDistantNeighbour(first_outermost)
            second_outermost = distantNeighbour[1]
            distance = distantNeighbour[0]
            better_candidate_index = None
            
            for node in self.container:
                if distance < node.distance:
                    better_candidate_index = node.index
                    distance = node.distance
            if better_candidate_index != None:
                outermosts = self.container[better_candidate_index].getOutermosts()
                if len(outermosts) != 2:
                    print(("Error1 updating outermosts",5/0))
                self.setOutermosts(outermosts,distance)
            else:
                if second_outermost == first_outermost:
                    self.setOutermosts([first_outermost])
                else:
                    self.setOutermosts([first_outermost,second_outermost],distance)
            # restore the outermosts that this node had before adding a new element followed by node splitting,
            # if the previous outermosts were better and still available
            if (previousOutermosts and
                len(previousOutermosts) == 3 and
                previousOutermosts[2] > self.distance and
                self.indexOf(previousOutermosts[0]) != None and
                self.indexOf(previousOutermosts[1]) != None):
                self.setOutermosts([previousOutermosts[0],previousOutermosts[1]],previousOutermosts[2])
            else:
                pass
                
        else:
            print(("Error2 updating outermosts",5/0))
        for node in self.container:
            if self.distance < node.distance:
                print(("Error3 updating outermosts",self.distance,node.distance))
                5/0

    # This method replaces the existing outermosts with the elements submitted in
    # the list of arguments and fix self.distance accordingly.
    # This method is used in the methods 'add', 'insert', 'split' and 'updateOutermosts'.
    # Errors:
    # 1) second supposed outermos elements are not present in this node
    # 3) First and second outermost elements can not be the same;
    # 4) Carring the method out can not result in emptying of the list of outermosts;
    # 5) If the node contains more than 1 leaf element, there must be 2 outermost elements.
    def setOutermosts(self,outermosts,distance=0):
        self.outermosts = []
        if outermosts and len(outermosts):
            for name in outermosts:
                self.outermosts.append(name)
        self.distance = distance
        if len(self.outermosts) == 2:
            if self.indexOf(self.outermosts[1]) == None:
                print(("Error1 setting outermosts", self.outermosts[1], 5/0))
            if self.outermosts[0]==self.outermosts[1]:
                print(("Error2  setting outermosts",5/0))
        # mark that this node and its child nodes must to be updated
        self.needUpdate = 1
        for node in self.container:
            node.needUpdate = 1
        if len(self.outermosts) == 0:
            print(("Error3 setting outermosts",5/0))
        if self.size() >= 2 and len(self.outermosts) < 2:
            print(("Error4 setting outermosts",5/0))
        for node in self.container:
            if self.distance < node.distance:
                print(("Error5 setting outermosts",self.distance,node.distance))
                5/0
    
    # This methos is called when the node has got a new element which is more distant to
    # the existing outermost elements than they are to each other. The new element adds
    # itself to the list of outermosts if this list contains only 1 element, or one of the
    # outermost elements is replaced with the new one.
    # The method takes arguments: (name of new element,name of the element to be replaced,
    #       [distance,name of the distant neighbour])
    # Updating of neighbours in in this node and all subordinate nodes must follow!!!
    # Only 'add' method calls this method.
    # Errors:
    # 1) The supposed distant neighbours is not present in this node
    # 2) First and second outermost elements can not be the same;
    # 3) Carring the method out can not result in emptying of the list of outermosts;
    # 4) If the node contains more than 1 leaf element, there must be 2 outermost elements.
    def fixOutermosts(self,name,distantNeighbour):
        self.previousOutermosts = []
        if len(self.outermosts):
            for val in self.outermosts:
                self.previousOutermosts.append(val)
            self.previousOutermosts.append(self.distance)
            
        if self.indexOf(distantNeighbour[1]) == None:
            print(("Error1 fixing outermosts", distantNeighbour, 5/0))
        # The closest outermost element to the given one
        closeNeighbour = min(self.getLocalDistances(name))[1]
        # fix the distance
        self.distance = distantNeighbour[0]
        # replace outermosts
        self.replaceOutermosts(closeNeighbour,name)
        # if distantNeighbour isn't already the second outermost, replace it
        if distantNeighbour[1] not in self.outermosts:
            index = self.outermosts.index(name)
            second_outermost = self.outermosts[index-1]
            self.replaceOutermosts(second_outermost,distantNeighbour[1])
        # mark that this node and its child nodes must to be updated
        self.needUpdate = 1
        for node in self.container:
            node.needUpdate = 1
        if len(self.outermosts) > 1:
            if self.outermosts[0]==self.outermosts[1]:
                print(("Error2  fixing outermosts", self.outermosts))
                5/0
        if len(self.outermosts) == 0:
            print(("Error3 fixing outermosts",5/0))
        if self.size() >= 2 and len(self.outermosts) < 2:
            print(("Error4 fixing outermosts",5/0))
        for node in self.container:
            if self.distance < node.distance:
                print(("Error5 fixing outermosts",self.distance,node.distance))
                5/0

    # This method is callable by the method 'fixOutermosts' only. It replaces in the 
    # list self.outermosts an oldName with a newName. It does't care for updating 
    # any other attributes of the current and subordinate nodes.
    # It throws errors if:
    # 1) oldName is not in the list self.outermosts or self.outermosts is empty;
    # 2) self.outermosts has got empty after replacing;
    # 3) after replacement two names of outermost elements in self.outermosts appeared
    #   to be the same.
    def replaceOutermosts(self,oldName,newName):
        if len(self.outermosts) == 1:
            self.outermosts.append(newName)
            return
        elif oldName in self.outermosts:
            index = self.outermosts.index(oldName)
            self.outermosts[index] = newName
        else:
            print(("Error replacing outermosts",oldName,newName,self.outermosts))
        if self.length:
            for item in self.container:
                item.setNeighbourhood()
        if len(self.outermosts) == 0:
            print(("Error1 replacing outermosts",5/0))
        if len(self.outermosts)==2 and self.outermosts[0]==self.outermosts[1]:
            print(("Error2 replacing outermosts",5/0))
        return
    
    # The methods of manipulation with neighbours
    # self.neighbours = [[distance,internalNeighbour,externalNeighbour],...]
    
    # This method set neighbours for the node. It calculate all possible distances
    # between outermost elements of the parent node and own outermost elements and
    # put this information to the list self.neighbours.
    # The method throws errors if:
    # 1) The list of outermost elements of the parent node is empty;
    # 2) The method ends up with an empty list self.neighbours
    # This method must be called every time when own outermosts or the outermost
    # elements where changed includin situation when the node was inserted into
    # the collection of another node.
    # The methods 'insert', 'setOutermosts', 'replaceOutermosts' and 'updateOutermosts'
    # use this method.
    def setNeighbourhood(self):
        if not self.isTopLevelNode():
            Neighbourhood = self.trigger("Get outermost")
            if len(Neighbourhood) == 0:
                print("Error1 setting neighbourhood")
            neighbours = []
            for member in Neighbourhood:
                elements = self.checkNeighbours(member)
                if elements == None:
                    localDistances = self.getLocalDistances(member)
                    for value in localDistances:
                        neighbours.append([value[0],value[1],member])
                else:
                    for value in elements:
                        neighbours.append([value[0],value[1],member])
            self.neighbours = []
            for item in neighbours:
                if item not in self.neighbours:
                    self.neighbours.append([item[0],item[1],item[2]])
            if len(self.neighbours) == 0:
                print(("Error2 setting neighbourhood",5/0))
        # update neighbours in all child nodes
        for node in self.container:
            # check if the child node needs update
            if node.needUpdate:
                node.setNeighbourhood()
        # the node was succesfully updated
        self.needUpdate = 0

    # An ancillary function to help the method 'setNeighbourhood' avoid unnecessary calculations.
    # If the entries [distance,internalOutermost,externalOutermost] are already in self.neighbours,
    # no need to calculate the distance againe. If there are any mismatches, the function returns None. 
    def checkNeighbours(self,externalNeighbour):
        elements = []
        if self.neighbours == None:
            return None
        for item in self.neighbours:
            if externalNeighbour in item:
                members = [item[1],item[2]]
                index = members.index(externalNeighbour)
                if members[index-1] in self.outermosts:
                    elements.append([item[0],members[index-1]])
                else:
                    return None
        if len(elements):
            return elements
        else:
            return None
    
    # Methods to get new information
    
    # This method calculate distances between a submitted element 'name'
    # and the outermost elements of the current node
    # It returns nested lists [[distance,first_outermost],[distance,second_outermos]]
    # This method is callable from the methods 'add', 'setNeighbourhood' and 'insertNeighbour'
    def getLocalDistances(self,name):
        outermosts = self.getOutermosts()
        if outermosts:
            first_name = outermosts[0]
        else:
            first_name = name
            outermosts = []
        #pattern = self.trigger("Get pattern",name)
        if first_name == name:
            first_distance = 0
        else:
            first_distance = self.trigger("Get distance",[name,first_name])
        if len(outermosts)==2:
            second_name = outermosts[1]
            if second_name == name:
                second_distance = 0
            else:
                second_distance = self.trigger("Get distance",[name,second_name])
            return [[first_distance,first_name],[second_distance,second_name]]
        else:
            return [[first_distance,first_name]]

    # Set the attribute 'needUpdate' to 0
    def resetUpdateNeed(self):
        self.needUpdate = 0
        for node in self.container:
            node.resetUpdateNeed()
    
    # ACCESSION METHODS
    
    # PRIVATE ACCESSION METHODS
    
    # This method is to make a copy of a node preserving node specific information
    # but cutting links with parent nodes. It creates an empty node and copies all 
    # attributes and child nodes in a recursive maner.
    # This method is used by the methods 'split', 'propagateNodes' and by itself.
    def copy(self):
        node = Node()
        # ATTRIBUTES
        outermosts = self.getOutermosts()
        if len(outermosts):
            node.outermosts = []
            for name in outermosts:
                node.outermosts.append(name)
        else:
            node.outermosts = []
        if self.previousOutermosts:
            node.previousOutermosts = []
            for value in self.previousOutermosts:
                node.previousOutermosts.append(value)
        node.neighbours = []
        if len(self.neighbours):
            for item in self.neighbours:
                node.neighbours.append([item[0],item[1],item[2]])
        node.level = self.level
        node.trigger = self.trigger
        node.index = self.index
        node.distance = self.distance
        node.name = self.name
        node.needUpdate = self.needUpdate
        node.container = []
        if self.length:
            for item in self.container:
                node.insert(item.copy(),item.name)
        return node
    
    # Return a list of outermost elements [name] or [name,name],
    # or None, but never []
    def getOutermosts(self):
        if self.outermosts == None:
            return None
        outermosts = []
        if len(self.outermosts):
            for name in self.outermosts:
                outermosts.append(name)
        if self.level and len(outermosts) == 0:
            #self.takeSnapShot()
            print(("Error in getOutermosts",self.level))
            5/0
        return outermosts
    
    def getDistance(self):
        return self.distance
    
    # In the lists of neighbours [[distance,internalOutermost,externalOutermost],...] 
    # of all subordinate nodes this method looks for the entry with given externalOutermost
    # and the bigest 'distance', then returns the distance and the name of the internalOutermost.
    def findDistantNeighbour(self,extOutermost):
        if self.length == 1:
            return self.container[0].findDistantNeighbour(extOutermost)
        neighbours = []
        for node in self.container:
            for item in node.neighbours:
                candidates = [item[1],item[2]]
                # check neighbours with the given externar outermost and a distance not 0
                if item[0] and extOutermost in candidates:
                    index = candidates.index(extOutermost)
                    # check if the found neighbour is still in this node
                    if self.indexOf(candidates[index-1]) == None:
                        continue
                    else:
                        neighbours.append([item[0],candidates[index-1]])
        if len(neighbours):
            neighbours.sort()
            return neighbours[-1]
        else:
            index = self.indexOf(extOutermost)
            return self.container[index].findDistantNeighbour(extOutermost)

    # Return list of neighbours [[distance,internalNeighbour,externalNeighbour],...]
    def getListOfNeighbours(self):
        neighbours = []
        if self.neighbours == None:
            return neighbours
        for item in self.neighbours:
            neighbours.append(item[2])
        return neighbours
    
    # This method returns the name of an outermost element of the parent node
    # which is the closest one among two to this node
    def getNeighbour(self):
        return min(self.neighbours)[2]
    
    # Return the neighbour [distance,internalNeighbour,externalNeighbour] with the smallest distance
    def getClosestNeighbour(self):
        return min(self.neighbours)
    
    # Return the high-level index of a leaf element in the current node
    def indexOf(self,name):
        if self.isLeafNode():
            if name==self.name:
                return 0
            else:
                return None
        else:
            for child in self.container:
                index = child.indexOf(name)
                if index != None:
                    return child.index
            return None
    
    # Return true if the node is last but one in the node hierarchy.
    # This node contains the collection of leaf nodes.
    def isLeavesContainer(self):
        thresholds = self.trigger("Get thresholds")
        if not thresholds:
            thresholds = self.thresholds
        if self.level == len(thresholds)-1:
            return 1
        else:
            return 0
    # Return true for leaf nodes
    def isLeafNode(self):
        thresholds = self.trigger("Get thresholds")
        if not thresholds:
            thresholds = self.thresholds
        if self.level == len(thresholds):
            return 1
        else:
            return 0
        
    # Return true if this node is on the top level, including the top level
    # of a copied branch of nodes.
    # In all top level nodes self.trigger == None
    def isTopLevelNode(self):
        if self.level:
            return 0
        else:
            return 1
        
    # PUBLIC ACCESSION METHODS
    
    # Return the number of leaf elements of the node
    def size(self):
        if self.isLeafNode():
            return 1
        if self.isLeavesContainer():
            return len(self.container)
        else:
            Size = 0
            for child in self.container:
                Size += child.size()
            return Size
    
    # Return the path to the element from the current node
    def getPath(self,name):
        thresholds = self.trigger("Get thresholds")
        if not thresholds:
            thresholds = self.thresholds
        if self.level < len(thresholds):
            for child in self.container:
                path = child.getPath(name)
                if path:
                    return str(child.index)+"."+path
            return None
        else:
            if len(list(name)) == len(list(self.name)) and name==self.name:
                return tools.basename(self.name)
            elif len(list(name)) < len(list(self.name)) and string.rfind(self.name,name)>-1:
                return name
            else:
                return None
    
    # Return the list of all leaf elements belonging to this node
    def getLeaves(self):
        if self.isLeafNode():
            return [self.name]
        else:
            leaves = []
            for child in self.container:
                child_leaves = child.getLeaves()
                for leaf in child_leaves:
                    leaves.append(leaf)
            return leaves
        
    # Return a list of unique outermost elements from all child nodes
    # Return a list of unique outermost elements from all child nodes
    def getListOfOutermosts(self, threshold=0):
        l = self.getOutermosts()
        if l == None:
            l = []
        if self.isLeafNode() or self.distance < threshold:
            return l
        else:
            for child in self.container:
                child_list = child.getListOfOutermosts(threshold)
                if not child_list:
                    continue
                for item in child_list:
                    if item not in l:
                        l.append(item)
            return l
        
    # This method is used by the graphical tree viewer. It returns the shift
    # of the node from left side of the canvas.
    def getShift(self,mode='threshold'):
        if mode == 'distance':
            return 100 - self.distance
        else:
            thresholds = self.trigger("Get thresholds")
            if not thresholds:
                thresholds = self.thresholds
            return 100 - thresholds[self.level]
    
    # Return number of brunches on the given level
    def getBranching(self,level):
        if self.level == level:
            return self.length-1
        elif self.level < level:
            branching = 0
            for node in self.container:
                branching += node.getBranching(level)
            return branching
        else:
            return 0

    # Transform the node tree to a Python dictionary and return it
    def valueOf(self):
        DB = {"level":self.level,
            "length":self.length,
            "name":self.name,
            "index":self.index,
            "distance":self.distance,
            "size":self.size(),
            "outermosts":self.getOutermosts(),
            "neighbours":self.neighbours,
            "container":[]}
        if self.length:
            for item in self.container:
                DB["container"].append(item.valueOf())
            return DB
        else:
            return DB
    
    # Importing the tree from a dictionary
    def importTree(self,db):
        self.level =  db["level"]
        self.length = 0
        self.name = db["name"]
        self.index = db["index"]
        self.distance = db["distance"]
        self.outermosts = []
        for name in db["outermosts"]:
            self.outermosts.append(name)
        if db["neighbours"]:
            self.neighbours = []
            for item in db["neighbours"]:
                self.neighbours.append([item[0],item[1],item[2]])
        for item in db["container"]:
            node = Node(0,self.do)
            node.importTree(item)
            self.insert(node,node.name)
        self.resetUpdateNeed()

    # Convert the database structure to a string and return it
    def toString(self):
        thresholds = self.trigger("Get thresholds")
        if not thresholds:
            thresholds = self.thresholds
        report = ("level: " + str(self.level) + "\n")
        if self.level == len(thresholds):
            report = report + "\t"*self.level + "name: " + self.name + "\n"
        report = report + (
                "\t"*self.level + "length: " + str(self.length) + "\n" +
                "\t"*self.level + "index: " + str(self.index) + "\n" +
                "\t"*self.level + "needUpdate: " + str(self.needUpdate) + "\n" +
                "\t"*self.level + "distance: " + str(self.distance) + "\n" +
                "\t"*self.level + "size: " + str(self.size()) + "\n" +
                "\t"*self.level + "outermosts: " + str(self.getOutermosts()) + "\n" +
                "\t"*self.level + "neighbour: " + str(self.neighbours) + "\n" +
                "\t"*self.level + "container: " + "\n")
        if self.length:
            for item in self.container:
                report = report + "\t"*self.level + "{" + item.toString() + "\n" + "\t"*self.level + "}\n"
        return report
    
    # Save the string presentation of the database structure
    def takeSnapShot(self):
        ouplib.TextEditor("Node snapshort",self.toString())
        
    # Temporary function to check the database structure and attributes.
    # It throws the errors when:
    # 1) attribute 'length' doesn't correspond the real length of the container;
    # 2) attribute 'needUpdate' is still on;
    # 3) there are more than 1 leaf element in the node but self.outermosts contains only 1 elements;
    # 4) an outermost element isn't presented in neighbours
    # 5) a neighbour element isn't presented in parent outermosts
    def check(self):
        if self.length != len(self.container):
            print(("Check error 1", self.toString(), 5/0))
        if self.needUpdate:
            print(("Check error 2", self.toString(), 5/0))
        if self.size() >= 2 and len(self.outermosts) < 2:
            print(("Check error 3", self.toString(), 5/0))
        for item in self.container:
            neighbours = item.getListOfNeighbours()
            for name in self.outermosts:
                if name not in neighbours:
                    print(("Check error 4", self.level, name, 5/0))
            internal_outermosts = item.getOutermosts()
            for val in item.neighbours:
                if val[1] not in internal_outermosts or val[2] not in self.outermosts:
                    print(("Check error 5", item.level))
                    print(val)
                    print((item.neighbours))
                    print((self.outermosts))
                    print((item.outermosts, 5/0))
            if self.distance < item.distance:
                print("Check error 6")
                print((self.distance,item.distance))
                print((self.outermosts))
                print((item.outermosts))
                5/0
            if item.index >= self.size():
                print("Check error 7")
                print((self.size()))
                print((item.index))
                5/0
            if self.level and not self.trigger:
                print("Check error 8")
                5/0
                    
            item.check()
        
###############################################################################
class Mapper:
    def __init__(self):
        # neighbour pattern
        self.weights = {"T":[2,0],"C":[3,1],"A":[0,2],"G":[1,3]}

        # reverse complement pattern
        # self.weights = {"T":[0,2],"C":[1,3],"A":[2,0],"G":[3,1]}
        
        self.code = {"X":["T","C","A","G"],"Y":["A","G","T","C"]}
        
    #### PUBLIC METHODS
    def __call__(self,word,x=None,y=None):
        # If the arguments are wlength,x,y - return the corresponding word
        if x and y and type(word)==type(1):
            return self._restore(word,x,y)
        elif type(word) == type([]) or type(word) == type(()):
            wlength,x,y = word
            return self._restore(wlength,x,y)
        wlength = len(list(word))
        if not wlength:
            return
        word = word.upper()
        return self._translate(word)
    
    # return list of permutations that differ from the given word from k1 to k2 nucleotides
    # w is either list [wlength,x,y] or a sequence 'ATGC'
    def get_permutations(self,w,k2=1,k1=1):
        wlength,x,y = self._parse(w)
        permutations = _get_permutations(wlength,x,y,k2,k1)
        for i in range(len(permutations)):
            wlength,x,y = permutations[i]
            x,y = self._restore_xy(wlength,x,y)
            permutations[i] = [wlength,x,y]
        return permutations

    # arguments w1 and w2 are either lists [wlength,x,y] or sequences 'ATGC'
    def compare_words(self,w1,w2):
        wlength,x,y = w1
        x,y = self._adjust_xy(wlength,x,y)
        w1 = [wlength,x,y]
        wlength,x,y = w2
        x,y = self._adjust_xy(wlength,x,y)
        w2 = [wlength,x,y]
        return self._count_mismatches(w1,w2)

    def right_subword(self,wlength,x,y,target_wl):
        if target_wl >= wlength:
            return wlength,x,y
        x,y = self._adjust_xy(wlength,x,y)
        self._right_subword(wlength,x,y,target_wl)
        x,y = self.restore_xy(wlength,x,y)
        return wlength,x,y
        
    def left_subword(self,wlength,x,y,target_wl):
        if target_wl >= wlength:
            return wlength,x,y
        x,y = self._adjust_xy(wlength,x,y)
        self._left_subword(wlength,x,y,target_wl)
        x,y = self.restore_xy(wlength,x,y)
        return wlength,x,y
    
    def move_right(self,w,shift=1):
        wlength,x,y = self._parse(w)
        if shift == 0:
            return [wlength,x,y]
        if shift > wlength-1:
            return None
        x,y = self._adjust_xy(wlength,x,y)
        words = self._move_right(wlength,x,y,shift)
        for i in range(len(words)):
            wlength,x,y = words[i]
            x,y = self._restore_xy(wlength,x,y)
            words[i] = [wlength,x,y]
        return words

    def move_left(self,w,shift=1):
        wlength,x,y = self._parse(w)
        if shift == 0:
            return [wlength,x,y]
        if shift > wlength-1:
            return None
        x,y = self._adjust_xy(wlength,x,y)
        words = self._move_left(wlength,x,y,shift)
        for i in range(len(words)):
            wlength,x,y = words[i]
            x,y = self._restore_xy(wlength,x,y)
            words[i] = [wlength,x,y]
        return words
    
    # Incrementation of words
    def right_increment(self,w,n):
        words = _right_increment(w,n)
        for i in range(len(words)):
            wlength,x,y = words[i]
            x,y = self._restore_xy(wlength,x,y)
            words[i] = [wlength,x,y]
        return words
    
    def left_increment(self,w,n):
        words = _left_increment(w,n)
        for i in range(len(words)):
            wlength,x,y = words[i]
            x,y = self._restore_xy(wlength,x,y)
            words[i] = [wlength,x,y]
        return words
    
    # return list of subwords of the length wl
    def right_shift(self,w,letters=[0,1,2,3]):
        wlength,x,y = self._parse(w)
        x,y = self._adjust_xy(wlength,x,y)
        words = self._right_shift(wlength,x,y,letters)
        for i in range(len(words)):
            wlength,x,y = words[i]
            x,y = self._restore_xy(wlength,x,y)
            words[i] = [wlength,x,y]
        return words
    
    def left_shift(self,w,letters=[0,1,2,3]):
        wlength,x,y = self._parse(w)
        x,y = self._adjust_xy(wlength,x,y)
        words = _left_shift(wlength,x,y,letters)
        for i in range(len(words)):
            wlength,x,y = words[i]
            x,y = self._restore_xy(wlength,x,y)
            words[i] = [wlength,x,y]
        return words

    def constituents(self,w,wl):
        wlength,x,y = self._parse(w)
        if wl == wlength:
            return [wlength,x,y]
        if wl > wlength:
            return None
        x,y = self._adjust_xy(wlength,x,y)
        words = _constituents(w,wl)
        for i in range(len(words)):
            wlength,x,y = words[i]
            x,y = self._restore_xy(wlength,x,y)
            words[i] = [wlength,x,y]
        return words
    
    # return 1 letter shorter word with 1 deletion
    def deletion(self,wlength,x,y,ind):
        if ind >= wlength:
            return None
        if ind == 0:
            return self.right_subword(wlength,x,y,wlength-1)
        if ind == wlength-1:
            return self.left_subword(wlength,x,y,wlength-1)
        x,y = self._adjust_xy(wlength,x,y)
        wlength,x,y = _deletion(wlength,x,y)
        x,y = self.restore_xy(wlength,x,y)
        return wlength,x,y
        
    # concatenate several words
    def concatenate(self,words):
        if not words:
            return words
        if len(words)==1:
            return words[0]
        wlength,x,y = words[0]
        x,y = self._adjust_xy(wlength,x,y)
        wlength,x,y = self._concatenate(wlength,x,y)
        x,y = self.restore_xy(wlength,x,y)
        return [wlength,x,y]
    
    def complement(self,w):
        wlength,x,y = self._parse(w)
        x,y = self._adjust_xy(wlength,x,y)
        wlength,x,y = self._complement(wlength,x,y)
        x,y = self.restore_xy(wlength,x,y)
        return [wlength,rX,rY]

    def revcomplement(self,w):
        wlength,x,y = self._parse(w)
        if wlength%2==0:
            return wlength,y,x
        else:
            wlength,x,y = self._revcomplement(wlength,x,y)
            x,y = self._restore_xy(wlength,x,y)
            return [wlength,x,y]        
        
    #### PRIVAT METHODS
    def _parse(self,w):
        if type(w)==type('ATGC'):
            wlength,x,y = list(self.__call__(w))
        else:
            wlength,x,y = w
        return [wlength,x,y]
    
    def _translate(self,word):
        # Calculate the list of [wlength,x,y] for a word with a variable N
        wlength = len(list(word))
        pos = word.find("N")
        if pos > -1:
            result = []
            words = [word]
            while pos > -1:
                for i in range(len(words)):
                    words[i] = words[i][:pos]+"A"+words[i][pos+1:]
                    for l in ("T","G","C"):
                        words.append(words[i][:pos]+l+words[i][pos+1:])
                pos = words[0].find("N")
            for w in words:
                result.append(self.__call__(w))
            return result
        # if word length is an odd number, the first letter in the word is added to the end
        if wlength%2==1:
            word += word[0]
        p = int(len(list(word))/2)
        word1 = list(word)[:p]
        word2 = list(word)[p:]
        x = self._accumulate_y(word1)
        y = self._accumulate_x(word2)
        x,y = self._restore_xy(wlength,x,y)
        return [wlength,x,y]
    
    # This functing transfers continuous word coordinates to parted coordinates
    def _adjust_xy(self,wlength,x,y):
        if wlength%2==1:
            if y%2 == 0:
                d = 0
            else:
                d = 2
            x += d+2*int((x-1)/2)
            if x%2 == 0:
                d = 1
            else:
                d = 0
            y += d+2*int((y-1)/2)+(y-1)%2
            return x,y
        else:
            return x,y
        
    # This functing restores parted word coordinates to continuous coordinates
    def _restore_xy(self,wlength,x,y):
        if wlength%2==1:
            x -= (int((x-1)%4)/2)*2 + 2*int((x-1)/4)
            if y%2 == 0:
                d = 1
            else:
                d = 0
            y -= int(((y-1)%4)/2) + d + 2*int((y-1)/4)
            return x,y
        else:
            return x,y
    
    def _restore(self,wlength,x,y):
        x,y = self._adjust_xy(wlength,x,y)
        MAXVAL = 2**(wlength+wlength%2)
        if x > MAXVAL or y > MAXVAL:
            return None
        if wlength%2==1 and abs(x%4 - y%4) != 2:
            return None
        x -= 1
        y -= 1
        if wlength%2==1:
            W = (wlength+1)//2
        else:
            W = wlength//2
        word1 = word2 = ""
        while W:
            p = int(4**(W-1))
            i = int(round(x//p))
            j = int(round(y//p))
            word1 = self.code["X"][i] + word1
            word2 += self.code["Y"][j]
            x -= i*p
            y -= j*p
            W -= 1
        if wlength%2==1:
            word2 = word2[:-1]
        return word1+word2
    
    # return list of permutations that differ from the given word from k1 to k2 nucleotides
    # w is either list [wlength,x,y] or a sequence 'ATGC'
    def _get_permutations(self,wlength,x,y,k2,k1):
        if wlength < k2:
            k2 = wlength
        permutations = []
        p = (wlength+wlength%2)/2-1
        x_permutations = self._nrange(x,p,k2)
        x_permutations.insert(0,[x])
        y_permutations = self._nrange(y,p,k2)
        y_permutations.insert(0,[y])
        nbX = nbY = None
        board = {}
        for m in range(len(x_permutations)):
            for u in range(len(y_permutations)):
                if m+u < k1 or m+u > k2:
                    continue
                for i in range(len(x_permutations[m])):
                    for j in range(len(y_permutations[u])):
                        rX = x_permutations[m][i]
                        rY = y_permutations[u][j]
                        # correction for odd words
                        if wlength%2==1:
                            t = (rX-1)%4
                            rY = (rY-(rY-1)%4) + abs(t-2-2*(t%2))
                        # check for uniqueness
                        if rX not in board:
                            board[rX] = []
                        if rY not in board[rX]:
                            board[rX].append(rY)
                        else:
                            continue
                        permutations.append([wlength,rX,rY])
        return permutations

    def _right_subword(self,wlength,x,y,target_wl):
        q = int((wlength+wlength%2)/2)-1
        t = int((y-1)/(4**q))
        wlength,x,y = self._right_intermediate(wlength,x,y)
        if wlength%2==0:
            # Add last letter to X - high level X component
            x += abs(t-2-2*(t%2))*(4**q)
            # Add last letter to Y - low level Y component
            y += abs((x-1)%4-2-2*(x-1)%2)
        if target_wl < wlength-1:
            return self._right_subword(wlength-1,x,y,target_wl)
        return wlength-1,x,y
        
    def _left_subword(self,wlength,x,y,target_wl):
        q = int((wlength+wlength%2)/2)-1
        t = int((x-1)/(4**q))
        wlength,x,y = self._left_intermediate(wlength,x,y)
        if wlength%2==0:
            # Add last letter to Y - low level Y component
            y += abs((x-1)%4-2-2*((x-1)%2))
        else:
            # Remove last letter - low level Y
            y = self._hdecrement(wlength-1,y)
            # Add firts letter to Y - high level Y component
            y += abs(t-2-2*(t%2))*(4**(q-1))
        if target_wl < wlength-1:
            return self._left_subword(wlength-1,x,y,target_wl)
        return wlength-1,x,y
    
    def _move_right(self,wlength,x,y,shift):
        words_toProcess = [[wlength,x,y]]
        for i in range(shift):
            words = []
            for W,X,Y in words_toProcess:
                words.extend(self._right_shift(W,X,Y))
            words_toProcess = []
            words_toProcess.extend(self._ridof_redundancy(words))
        return words

    def _move_left(self,wlength,x,y,shift):
        words_toProcess = [[wlength,x,y]]
        for i in range(shift):
            words = []
            for W,X,Y in words_toProcess:
                words.extend(self._left_shift(W,X,Y))
            words_toProcess = []
            words_toProcess.extend(self._ridof_redundancy(words))
        return words
    
    # Incrementation of words
    def _right_increment(self,w,n):
        words = []
        for word in self._generate(n):
            wlength,x,y = self._add(w,word)
            words.append([wlength,x,y])
        return words
    
    def _left_increment(self,w,n):
        words = []
        for word in self._generate(n):
            wlength,x,y = self._add(word,w)
            words.append([wlength,x,y])
        return words
    
    def _constituents(self,w,wl):
        words = [list(self.left_subword(wlength,x,y,wl))]
        for i in range(wl,wlength):
            W,X,Y = words[-1]
            next_subword = self._right_shift(W,X,Y,[self._y_index(wlength,x,y,i)])[0]
            W,X,Y = next_subword
            X,Y = self.restore_xy(W,X,Y)
            words.append([W,X,Y])
        return words
    
    # return 1 letter shorter word with 1 deletion
    def _deletion(self,wlength,x,y,ind):
        q = int((wlength+wlength%2)/2)-1
        if wlength%2==0:
            if ind <= int((wlength+wlength%2)/2)-1:
                t = int((y-1)/(4**q))
                # Delete one letter
                x = self._hslice(wlength,x,ind)+self._hdecrement(wlength,x,ind)-1
                # Add last letter to X - high level X component
                x += abs(t-2-2*(t%2))*(4**q)
                
                # Remove first after middle letter - high level Y component
                y -= t*(4**q)
                # Increase levels of Y components
                y = self._hincrement(wlength,y)
            else:
                # Delete one letter
                ind = wlength-ind-1
                y = int(y/(4**(ind+1)))*(4**(ind+1))+self._hincrement(wlength,y,ind)
                
            # Add last letter to Y - low level Y component
            y += abs((x-1)%4-2-2*(x-1)%2)
        else:
            if ind <= int((wlength+wlength%2)/2)-1:
                # Delete one letter
                x = self._hslice(wlength,x,ind)+self._hdecrement(wlength,x,ind)-1
                # Decrease levels of Y components
                y = self._hdecrement(wlength,y)
            else:
                t = int((x-1)/(4**q))
                # Remove first before middle letter - high level X component
                x -= t*(4**q)

                # Decrease levels of Y components
                y = self._hdecrement(wlength,y)
                # Delete one letter
                ind = wlength-ind-1
                y = self._hslice(wlength,y,ind)+self._hdecrement(wlength,y,ind)-1
                # Add first letter to Y - high level Y component
                y += abs(t-2-2*(t%2))*(4**(q-1))
        wlength -= 1
        return wlength,x,y
        
    # concatenate several words
    def _concatenate(self,wlength,x,y):
        word = [wlength,x,y]
        for i in range(1,len(words)):
            wlength,x,y = words[i]
            x,y = self._adjust_xy(wlength,x,y)
            word = self._add(word,[wlength,x,y])
        wlength,x,y = word
        return [wlength,x,y]
    
    def _complement(self,wlength,x,y):
        rX = rY = 0
        for p in range(int((wlength+wlength%2)/2)-1,-1,-1):
            step = 4**p
            rX += step*abs(int(x/step) - 2*(1+(int(x/step))%2))
            rY += step*abs(int(y/step) - 2*(1+(int(y/step))%2))
            x -= step*int(x/step)
            y -= step*int(y/step)
        return [wlength,rX,rY]

    def _revcomplement(self,wlength,x,y):
        if wlength%2==0:
            return wlength,y,x
        else:
            x,y = self._adjust_xy(wlength,x,y)
            W = wlength+1
            words = self._right_shift(W,y,x)
            for wl,wx,wy in words:
                if abs(wx%4 - wy%4) == 2:
                    return wlength,wx,wy

    # Create a list of all possible words of the given length
    def _generate(self,wlength):
        words = []
        p = 2**(wlength+wlength%2)
        for x in range(p):
            y = 0
            while y < p:
                if wlength%2==1:
                    t = x%4
                    rY = y - y%4 + abs(t-2-2*(t%2))
                    y = y - y%4 + 4
                else:
                    rY = y
                    y += 1
                words.append([wlength,x+1,rY+1])
        return words
    
    def _count_mismatches(self,first,second):
        w,x,y = self._parse(first)
        W,X,Y = self._parse(second)
        mismatches = []
        if W == w:
            mismatches.append(self._substract(x,X)+self._substract(y,Y))
        elif W > w:
            words = self.constituents([W,X,Y],w)
            for word in words:
                d,s = self.count_mismatches(first,word)
                mismatches.append(d)
        else:
            words = self.constituents([w,x,y],W)
            for word in words:
                d,s = self.count_mismatches(word,second)
                mismatches.append(d)
        match = min(mismatches)
        indices = []
        for i in range(len(mismatches)):
            if mismatches[i]==match:
                indices.append(i)
        return match,indices
    
    def _substract(self,first,second):
        if first==second:
            return 0
        count = 0
        p = 1
        d = first-second
        while d:
            p = 4.0**int(math.floor(math.log(abs(d),4)))
            k = round(abs(d)/p)
            first -= int(k*p*abs(d)/d)
            count += 1
            d = first-second
        return count

    def _right_shift(self,wlength,x,y,letters=[0,1,2,3]):
        words = []
        if wlength == 1:
            for y in letters:
                words.append([1,abs(y-2-2*(y%2))+1,y+1])
            return words
        q = int((wlength+wlength%2)/2)-1
        t = int((y-1)/(4**q))
        W,x,y = self._right_intermediate(wlength+wlength%2,x,y)
        # Add last letter to X - high level X component
        x += abs(t-2-2*(t%2))*(4**q)
        if wlength%2==1:
            y = 16*int(y/16)+1
            for i in range(3):
                if  abs(x%4 - y%4) == 2:
                    break
                y += 1
            m = 4
        else:
            m = 1
        for i in letters:
            words.append([wlength,x,y+m*i])
        return words
    
    def _left_shift(self,wlength,x,y,letters=[0,1,2,3]):
        words = []
        if wlength == 1:
            for y in letters:
                words.append([1,abs(y-2-2*(y%2))+1,y+1])
            return words
        q = int((wlength+wlength%2)/2)-1
        t = int((x-1)/(4**q))
        W,x,y = self._left_intermediate(wlength+wlength%2,x,y)
        # Increase levels of X components and normalize
        x = 4*int(self._hincrement(wlength,x)/4)+1
        # Remove last letter
        y = self._hdecrement(wlength,(y-1))
        # Add firts letter to Y - high level Y component
        y += abs(t-2-2*(t%2))*(4**q)
        if wlength%2==1:
            # normalize
            y = 4*int(y/4)+1
            for i in letters:
                words.append([wlength,x+i,y+abs(i-2-2*(i%2))])
        else:
            y += 1
            for i in letters:
                words.append([wlength,x+i,y])
        return words

    def _accumulate_x(self,word): # word must be list
        wlength = len(word)
        x = 1
        for i in range(len(list(word))):
            x += self._get_x(word[i],wlength)
            wlength -= 1
        return x

    def _accumulate_y(self,word): # word must be list
        wlength = len(word)
        y = 1
        for i in range(len(list(word))-1,-1,-1):
            y += self._get_y(word[i],wlength)
            wlength -= 1
        return y
    
    def _get_x(self,letter,wlength):
        k = 4**(wlength-1)
        if letter not in self.weights:
            return
        return self.weights[letter][0]*k

    def _get_y(self,letter,wlength):
        k = 4**(wlength-1)
        if letter not in self.weights:
            return
        return self.weights[letter][1]*k
    
    def _ridof_redundancy(self,items):
        items.sort()
        for i in range(len(items)-1,0,-1):
            if items[i] == items[i-1]:
                del items[i]
        return items

    def _right_intermediate(self,wlength,x,y):
        #GTGA 4,4,5
        # First after middle letter become last befor middle
        # For example, in GTGA -> TG-A* 
        # high level Y component become high leve X component
        q = int((wlength+wlength%2)/2)-1
        t = int((y-1)/(4**q))
        # Remove first after middle letter - high level Y component
        if wlength%2==0:
            y -= t*(4**q)
        # Remove first letter - low level X component
        x -= (x-1)%4
        # Decrease levels of letters in X component
        x = self._hdecrement(wlength,x)
        # and recalculate leveles of Y components
        if wlength%2==0:
            y = self._hincrement(wlength,y)
        else:
            y = self._hdecrement(wlength,y)
        return wlength,x,y
    
    # private method to support other public methods
    def _left_intermediate(self,wlength,x,y):
        q = int((wlength+wlength%2)/2)-1
        t = int((x-1)/(4**q))
        # Remove first before middle letter - high level X component
        if wlength%2==1:
            x -= t*(4**q)
        # Remove last letter - low level Y component
        y -= (y-1)%4
        # Double decrease leveles of Y components
        if wlength%2==1:
            y = self._hdecrement(wlength,y)-1
        return wlength,x,y

    def _hincrement(self,wlength,val,IND=0):
        H = 1
        if not IND:
            IND = int((wlength+wlength%2)/2)-1
        for p in range(IND,0,-1):
            d = (int(val-1)/(4**(p-1)))
            val -= d*(4**(p-1))
            H += d*(4**p)
        H = H%(2**(wlength+wlength%2))
        return H
        
    # hyperdecrement until the index specified by IND
    def _hdecrement(self,wlength,val,IND=0):
        H = 1
        for p in range(int((wlength+wlength%2)/2)-1,IND,-1):
            d = int((val-1)/(4**p))
            val -= d*(4**p)
            H += d*(4**(p-1))
        if H==0:
            H=1
        return H
    
    # hyperslice from IND to zero level
    def _hslice(self,wlength,val,IND):
        H = 1
        for p in range(int((wlength+wlength%2)/2)-1,-1,-1):
            d = int((val-1)/(4**p))
            val -= d*(4**p)
            if p < IND:
                H += d*(4**p)
        if H==0:
            H=1
        return H
    
    # return a letter code by y-component
    def _y_index(self,wlength,x,y,ind):
        if ind >= wlength:
            return None
        p = int((wlength+wlength%2)/2) - ind
        if p <= 0:
            return ((y-1)/4**(int((wlength+wlength%2)/2)-1+p))%4
        else:
            t = ((x-1)/(4**(int((wlength+wlength%2)/2)-p)))%4
            return abs(t-2-2*(t%2))
        
    # concatenate 2 words
    def _add(self,first,second):
        if type(first)==type('ATGC'):
            wl1,x,y = list(self.__call__(first))
            first = [wl1,x,y]
        else:
            wl1,x,y = first
        if type(second)==type('ATGC'):
            wl2,x,y = list(self.__call__(second))
            second = [wl2,x,y]
        else:
            wl2,x,y = second
        wlength = wl1+wl2
        indx = (wl2-wl1)/2
        indy = indx+1
        x = y = 1
        for p in range(int((wlength+wlength%2)/2)-1,-1,-1):
            # calculate X
            if indx <= 0:
                W,X,Y = first
                t = self._y_index(W,X,Y,wl1+indx-1)
            else:
                W,X,Y = second
                t = self._y_index(W,X,Y,indx-1)
            x += int(abs(t-2-2*(t%2))*(4**p))
            indx -= 1
            # calculate Y
            if indy <= 0:
                W,X,Y = first
                t = self._y_index(W,X,Y,wl2+indy-1)
            else:
                W,X,Y = second
                t = self._y_index(W,X,Y,indy-1)
            y += int(t*(4**p))
            indy += 1
        return wlength,x,y
    
    # Next two functions are used together with the get_permutations function to 
    # create a list of permutations
    def _nrange(self,n,p,k):
        if k==0 or p<0:
            return [[n]]
        dN = 0
        rP = p
        index_before = index_after = 0
        permutations = [[]]
        orders = []
        while p >= 0:
            step = 4**p
            orders.append(int((n-dN-1)/step))
            first = n - step*orders[-1]
            added_after = 0
            for i in range(4):
                rN = first + step*i
                if rN != n:
                    if rN > n:
                        permutations[-1].insert(len(permutations[-1])-index_after,rN)
                        added_after += 1
                    else:
                        permutations[-1].insert(index_before,rN)
                        index_before += 1
            index_after += added_after
            dN += step*orders[-1]
            p -= 1
        if k>1:
            permutations.extend(self._populate(permutations[-1],orders,rP,k))
        return permutations
            
    def _populate(self,seeds,orders,p,k):
        indices = []
        elements = []
        indices.extend(seeds)
        for m in range(len(orders)-1):
            order = orders[m]
            step = 4**p
            n_elements = [0]
            for i in range(3,-1,-1):
                if order==i:
                    continue
                elif i < order:
                    index = 0
                else:
                    index = -1
                shift = i*step
                if not n_elements[0]:
                    n_elements[0] = indices.pop(index)-shift
                n_elements.append(shift)
            n_permutations = self._nrange(n_elements[0],p-m-1,k-1)
            n_permutations[0] = n_permutations[0]*3
            for j in range(len(n_permutations[0])):
                if j < int(len(n_permutations[0])/3):
                    shift = n_elements[1]
                elif j >= 2*int(len(n_permutations[0])/3):
                    shift = n_elements[3]
                else:
                    shift = n_elements[2]
                n_permutations[0][j] += shift
            elements.append(n_permutations[0])
            p -= 1
        return elements

###############################################################################
class Pattern:
    def __init__(self,wlength,title=''):
        self.db = word_db()
        self.oNormalizationTable = None
        self.wlength = wlength
        self.seqLength = 0
        self.title = title
        self.pattern_type = "n"
        self.normalization = 0
        self.ouv = 0.0
        self.ps = None
        self.colorscale = 12
            # child patterns
        self.child_patterns = []
        self._initiate()
        
    def __iter__(self):
        # Return list of word objects
        iterator = iter(self.getWordList())
        return iterator
    
    def __getitem__(self,key):
        return self.getWord(key)

    def __sub__(self,other):
        return self._compare(other)

    def __ror__(self, other):
        words = []
        my_wordlist = self.getWordList()
        other_wordlist = other.getWordList()
        for i in range(len(my_wordlist)):
            words.append([my_wordlist[i]-other_wordlist[i], # rank difference
                my_wordlist[i].word,                        # word
                str(my_wordlist[i].rank),                   # self rank
                format_number(float(my_wordlist[i].count)/self.seqLength,1,4),     # self expectation per 10 kbp
                str(my_wordlist[i].rank),                   # other rank
                format_number(float(other_wordlist[i].count)/other.seqLength,1,4), # other expectation per 10 kbp
                ])
        return words

    # Deviation between patterns
    def _initiate(self):
        # Fill in the generic pattern with words
        for i in range(1,2**self.wlength+1):
            for j in range(1,2**self.wlength+1):
                self.db.add([self.wlength,i,j],0,Word(self.db.convertor(self.wlength,i,j)))
    
    def _populate_db(self,strSeq):
        def get_word_deviation(oWord):
            return oWord.deviation
        
        i = 0
        self.seqLength = len(list(strSeq))
        words = []
        word_objects = []
        while i+self.wlength <= self.seqLength:
            words.append(strSeq[i:i+self.wlength].upper())
            i += 1
        words.sort()
        word = ""
        count = 0
        for i in range(len(words)-1,-1,-1):
            if words[i] != word:
                flg_added = self.db.add(word,count)
                if flg_added:
                    word_objects.append(Word(word,count,self.oNormalizationTable.calculate_expectation(word),
                                                float(self.seqLength)/4**self.wlength,flg_added))
                else:
                    pass
                count = 1
                word = words[i]
            else:
                count += 1
        flg_added = self.db.add(word,count)
        if flg_added:
            word_objects.append(Word(word,count,self.oNormalizationTable.calculate_expectation(word),
                                        float(self.seqLength)/4**self.wlength,flg_added))
        if len(word_objects) > 2:
            word_objects.sort(key=get_word_deviation,reverse=True)
        sum = 0
        sum_sq = 0
        N = len(word_objects)
        for i in range(N):
            word_objects[i].set_rank(i+1)
            self.db.set_data(word_objects[i].info,word_objects[i])
            sum += word_objects[i].deviation
            sum_sq += word_objects[i].deviation**2
        if N:
            self.ouv = math.sqrt((sum_sq - sum**2/float(N))/float(N - 1))

    def _compare(self,other,flg_ps=False,flg_mindist=True):
        if self.wlength != other.wlength:
            return
        distances = [self._substract(self,other,flg_ps)]
        if flg_mindist:
            self_reverse = self.reverse_complement()
            other_reverse = other.reverse_complement()
            distances.append(self._substract(self,other_reverse))
            distances.append(self._substract(self_reverse,other_reverse))
            distances.append(self._substract(self_reverse,other))
        if not distances:
            return
        return min(distances)
    
    def _substract(self,first,second,flg_ps=False):
        D = 0
        Dmax = float(4**first.wlength)*(4**first.wlength - 1)/2.0
        Dmin = 0
        if flg_ps:
            if self.wlength % 2 == 0:
                Dmin = (4**first.wlength) - (2**first.wlength)
            else:
                Dmin = 0
        for i in range(1,2**first.wlength):
            for j in range(1,2**first.wlength):
                D += first.getWord([first.wlength,i,j])-second.getWord([second.wlength,i,j])
        return 100.0*(D-Dmin)/(Dmax-Dmin)
                
    def _getColor(self, value, colorNum='All', colorMode='Red-Grey-Blue'):
        if value == None:
            return "black"
        if value > 0:
            value *= 2.0
        if type(value) == type(''):
            AGCT_colors = {"A":[0xFF,0xFF,0],"G":[0xFF,0,0xFF],"C":[0,0xFF,0],"T":[0,0xFF,0xFF]}
            for lett in 'AGCT':
                for i in range(3):
                    AGCT_colors[lett][i] = int(AGCT_colors[lett][i])/len(list(value))
            RGB = [0,0,0]
            for l in value:
                for i in range(3):
                    RGB[i] = RGB[i] + AGCT_colors[l][i]
        elif colorMode=='Grey-Red':
            k = float(value-1)/float(self.colorscale)
            if k > 1:
                RGB = [178,0,0]
            else:
                RGB = [230+int(k*25),
                       230-int(k*230),
                       230-int(k*230)]
        else:
            if colorNum != 'All':
                colorNum = int(colorNum)
                if colorNum == 3:
                    if value > 3:
                        value = self.colorscale
                    elif value < -3:
                        value = -self.colorscale
                    else:
                        value = 0
                else:
                    step = 2*self.colorscale/colorNum
                    value = float(value/step)
                    if value > 0:
                        value = math.ceil(value)*step
                    else:
                        value = math.floor(value)*step
                    
            value = value*.5/self.colorscale
            RGB = [0,0,0]
            if colorMode == 'Red-Green-Blue':
                RGB[0] = int((.5 + value)*255)
                RGB[1] = int((1.0 - 2*abs(value))*255)
                RGB[2] = int((.5 - value)*255)
                for i in range(3):
                    if RGB[i] > 255:
                        RGB[i] = 255
                    elif RGB[i] < 0:
                        RGB[i] = 0
            elif colorMode == 'Red-Yellow-Blue':
                if value <=-0.05:
                    RGB[0] = int((1.0 + 4.0*value)*255)
                else:
                    RGB[0] = 255
                RGB[1] = int((1.0 - 2*abs(value))*255)
                RGB[2] = int(510.0*abs(value))
                for i in range(3):
                    if RGB[i] > 255:
                        RGB[i] = 255
                    elif RGB[i] < 0:
                        RGB[i] = 0
            elif colorMode == 'Red-Grey-Blue':
                if value == 0:
                    RGB[0] = RGB[1] = RGB[2] = 255
                elif value > 0.15 and value < .5:
                    RGB[0] = 222
                    RGB[1] = int((1.0 - 2*(value-0.15))*222)
                    RGB[2] = int((1.0 - 2*(value-0.15))*222)
                elif value < -0.15 and value > -.5:
                    RGB[0] = int((1.0 + 2.0*(value+0.15))*222)
                    RGB[1] = int((1.0 + 2.0*(value+0.15))*222)
                    RGB[2] = 222
                elif value >= .5:
                    RGB[0] = 160
                    RGB[1] = 10
                    RGB[2] = 10
                elif value <= -.5:
                    RGB[0] = 10
                    RGB[1] = 10
                    RGB[2] = 160
                else:
                    RGB[0] = 222
                    RGB[1] = 222
                    RGB[2] = 222
            elif colorMode == 'Grey-Green':
                if value < 0:
                    value = -value
                if value > 0.15 and value < .5:
                    RGB[0] = 222
                    RGB[1] = int((1.0 - 2*(value-0.15))*222)
                    RGB[2] = 222
                elif value >= .5:
                    RGB[0] = 10
                    RGB[1] = 160
                    RGB[2] = 10
                else:
                    RGB[0] = 222
                    RGB[1] = 222
                    RGB[2] = 222
            else:
                print("Error in color.getColor")
                return None
            
        color = "#%02X%02X%02X" % (RGB[0],RGB[1],RGB[2])
        return color

    def setPattern(self, strSeq, normalization=0, pattern_type="n"):
        self.normalization = normalization
        self.pattern_type = pattern_type
        self.oNormalizationTable = Pattern(self.normalization)
        if self.wlength:
            self.oNormalizationTable.setPattern(strSeq)
        self._populate_db(strSeq)
        
    def calculate_expectation(self,word):
        if len(list(word)) < self.wlength:
            return 0
        if self.wlength == 0:
            return float(self.seqLength)/(4**len(list(word)))
        subword = word[:self.wlength]
        f = self.db.get(subword)
        i = 1
        while i+self.wlength <= len(list(word)):
            if self.wlength:
                denominator = 0
                for w in self.db.convertor.right_shift(subword):
                    denominator += self.db.get(w)
            else:
                denominator = self.seqLength
            subword = word[i:i+self.wlength]
            f *= float(self.db.get(subword))/denominator
            i += 1
        return f
    
    def reverse_complement(self):
        oReversePattern = Pattern(self.wlength,self.title)
        for w in self.db.get_allWords():
            rw = self.db.convertor.revcomplement(w)
            oWord = self.db.get_data(w)
            flg_added = oReversePattern.db.add(rw,oWord.count)
            if flg_added:
                oRWord = Word(self.db.convertor(rw),oWord.count,oWord.expectation,oWord.stdev)
                oRWord.set_rank(oWord.rank)
                oReversePattern.db.set_data(rw,oRWord)
        return oReversePattern

    def getWord(self,w):
        return self.db.get_data(w)
    
    def getOUV(self):
        return self.ouv
    
    def getPS(self):
        if self.ps == None:
            self.ps = self._substract(self,self.reverse_complement(),True)
        return self.ps
    
    def getPatternType(self):
        return "%s%i_%imer" % (self.pattern_type,self.normalization,self.wlength)
    
    def add_child(self,oPattern):
        self.child_patterns.append(oPattern)
    
    def has_child(self,name):
        for child in self.child_patterns:
            if child.title == name:
                return True
        return False
        
    def get_child(self,name):
        for child in self.child_patterns:
            if child.title == name:
                return child
        return None
    
    def getFirstChild(self):
        if not self.child_patterns:
            return
        return self.child_patterns[0]
    
    def getWordList(self):
        wlist = []
        for i in range(1,2**self.wlength+1):
            for j in range(1,2**self.wlength+1):
                w = [self.wlength,i,j]
                oWord = self.db.get_data(w)
                wlist.append(oWord)
        return wlist
    
    def set_colorscale(self,value):
        self.colorscale = value
        
    def tostring(self):
        def get_word_deviation(oWord):
            return oWord.deviation
        
        words = self.getWordList()
        if not words:
            return
        words.sort(key=get_word_deviation,reverse=True)
        output = ["%s\n%s; OUP = %f; PS = %f\n" % (self.title,self.getPatternType(),self.getOUV(),self.getPS())]
        output.append("\t".join(["#","Word","Frequ.","Expect.","Deviation"]))
        for oWord in words:
            output.append("\t".join([str(oWord.rank),oWord.word,
                                        str(oWord.count),
                                        str(int(oWord.expectation)),
                                        format_number(oWord.deviation,4)]))
        return "\n".join(output)
    
    def svg(self,width=0,height=0,X0=0,Y0=0,flg_finish=True):
        svg = [] 
        if not width or width < 400:
            width = 400
        if height != width:
            height = width
        n = 2**self.wlength
        x0 = X0+10
        y0 = Y0+height/6.0
        cell_side = float(width-width/6-x0)/n
        words = self.getWordList()
        deviations = []
        for oWord in words:
            deviations.append(oWord.deviation)
        self.colorscale = 2.5*getMeanAndStDev(deviations)[1]
        if not self.colorscale:
            return svg
        for oWord in words:
            wlength,i,j = oWord.info
            i -= 1
            j -= 1
            val = oWord.deviation
            if val == None:
                val = -3*self.colorscale
            x1 = x0+i*cell_side
            x2 = x1+cell_side
            y1 = y0+j*cell_side
            y2 = y1+cell_side
            fillcolor = self._getColor(oWord.deviation)
            if wlength < 4:
                strokecolor = "green"
            else:
                strokecolor = "none"
            svg.append("<a title=\"%s f=%d, dev=%f, rank=%d;\">" % (oWord.word,oWord.count,val,oWord.rank))
            svg.append("<path  style=\"fill:%s;stroke:%s\" d=\"M%d,%dL%d,%dL%d,%dL%d,%dZ\" />" % 
                (fillcolor,strokecolor,x1,y1,x2,y1,x2,y2,x1,y2))
            svg.append("</a>")

        svg += self.svg_colorscale(width/6,n*cell_side,x0+n*cell_side,y0)
        y = Y0+height/18.0
        svg.append("<text x=\"%f\" y=\"%f\" style=\"text-anchor:middle\" font-size=\"10\" font-weight=\"bold\">%s</text>" % 
            (width/2.0,y,self.title))
        y = Y0+2.0*height/18.0
        svg.append("<text x=\"%f\" y=\"%f\" style=\"text-anchor:middle\" font-size=\"12\" font-weight=\"bold\">%s; OUP = %f; PS = %f</text>" % 
            (width/2.0,y,self.getPatternType(),self.getOUV(),self.getPS()))
        if flg_finish:
            svg.insert(0,"<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % (width+305,height*series,width+305,height*series))
            svg.append("</svg>")
            return "\n".join(svg)
        return svg
    
    def svg_colorscale(self,width,height,X0=0,Y0=0):
        svg = []
        y_shift = height/30.0
        step = 2.0*self.colorscale/24.0
        val0 = (self.colorscale+3.0*step)/2.0
        x1 = X0+width/4.0
        x2 = X0+2.5*width/4.0
        y1 = Y0
        for i in range(30):
            val = val0-step*i
            if val < 0:
                val /= 1.5
            fillcolor = self._getColor(val,self.colorscale)
            y2 = y1+y_shift
            svg.append("<a title=\"%d\">" % (int(val)))
            svg.append("<path  style=\"fill:%s;stroke:none\" d=\"M%d,%dL%d,%dL%d,%dL%d,%dZ\" />" % 
                (fillcolor,x1,y1,x2,y1,x2,y2,x1,y2))
            svg.append("</a>")
            y1 = y2
        y = Y0+1.5*y_shift
        svg.append("<text x=\"%f\" y=\"%f\" style=\"text-anchor:left\" font-size=\"12\">%d</text>" % (x2+2,y,int(self.colorscale/2.0)))
        y = Y0+9.5*y_shift
        svg.append("<text x=\"%f\" y=\"%f\" style=\"text-anchor:middle\" font-size=\"12\">%d</text>" % (x2+10,y,0))
        y = Y0+29.5*y_shift
        svg.append("<text x=\"%f\" y=\"%f\" style=\"text-anchor:left\" font-size=\"12\">-%d</text>" % (x2+2,y,int(self.colorscale)))
        return svg
    
###############################################################################

class Word:
    def __init__(self,word,count=0,expectation=0,stdev=0,info=[]):
        self.word = word
        self.count = count
        self.expectation = expectation
        self.stdev = stdev
        self.deviation = self.get_deviation()
        self.rank = 0
        self.info = info
        
    def __sub__(self,other):
        return abs(self.rank-other.rank)
        
    def get_deviation(self):
        # normstdiv = 1.0+10.0*self.getWordLength()/math.sqrt(self.getSeqLength())
        if not self.expectation:
            return None
        normstdiv = 1.0
        value = self.count
        if value == 0:
            value = 1
        numerator = math.log((value**2)*math.sqrt(self.stdev**2 + self.expectation**2)/(self.expectation**2)/math.sqrt(self.stdev**2 + value**2))
        denumerator = math.sqrt(math.log((self.stdev/self.expectation)**2 + 1))
        return 6.0*numerator/denumerator/normstdiv
    
    def set_rank(self,rank):
        self.rank = rank

###############################################################################
class word_db:
    def __init__(self):
        # {wlength:{'x':{y1:{},y2,...},'y':{x1:{},x2,...}}
        self.db = {}
        # Version of the database
        self.version = 1
        # Data of the database creation
        self.date = ""
        # Convertor
        self.convertor = Mapper()

    #### ACCESSIONS
    def parse(self,w):
        if type(w)==type('ATGC'):
            try:
                wlength,x,y = list(self.convertor(w))
            except:
                return
        else:
            wlength,x,y = w
        return [wlength,x,y]
        
    
    def add(self,w,count=1,data=None):
        try:
            wlength,x,y = self.parse(w)
        except:
            return
        new_entry = {'counter':0,'data':1,'info':[]}
        if wlength not in self.db:
            self.db[wlength] = {'x':{},'y':{}}
        if x not in self.db[wlength]['x']:
            self.db[wlength]['x'][x] = {y:new_entry}
        if y not in self.db[wlength]['y']:
            self.db[wlength]['y'][y] = {x:new_entry}
        if y not in self.db[wlength]['x'][x] or x not in self.db[wlength]['y'][y]:
            self.db[wlength]['x'][x][y] = self.db[wlength]['y'][y][x] = new_entry
        self.db[wlength]['x'][x][y]['counter'] += count
        self.db[wlength]['x'][x][y]['data'] = data
        return [wlength,x,y]
    
    def get(self,w):
        try:
            wlength,x,y = self.parse(w)
        except:
            return None
        if wlength not in self.db:
            return None
        if x not in self.db[wlength]['x']:
            return None
        if y not in self.db[wlength]['x'][x]:
            return None
        return self.db[wlength]['x'][x][y]['counter']
    
    def pop(self,w,flg_complement=False):
        wlength,x,y = self.parse(w)
        if self.has(wlength,x,y):
            count = self.getvalue(wlength,x,y,flg_complement)
            self.delete(wlength,x,y,self.db[wlength]['x'][x][y]['counter'])
        else:
            count = 0
        if flg_complement and self.has(wlength,y,x):
            wlength,cx,cy = self.convertor.revcomplement([wlength,x,y])
            if self.has(wlength,cx,cy):
                self.delete(wlength,cx,cy,self.db[wlength]['x'][cx][cy]['counter'])
        return count
    
    def set_data(self,w,val):
        wlength,x,y = self.parse(w)
        if not self.has([wlength,x,y]):
            return
        self.db[wlength]['x'][x][y]['data'] = val
        
    def get_data(self,w,flg_complement=False):
        wlength,x,y = self.parse(w)
        if wlength in self.db and x in self.db[wlength]['x'] and y in self.db[wlength]['x'][x] and not flg_complement:
            return self.db[wlength]['x'][x][y]['data']
        if flg_complement:
            wlength,x,y = convertor.revcomplement([wlength,x,y])
            return self.get_data(wlength,x,y)
        return None
    
    def append_info(self,w,info):
        wlength,x,y = self.parse(w)
        if not self.has([wlength,x,y]):
            return
        self.db[wlength]['x'][x][y]['info'].append(info)
        
    def get_info(self,w,flg_revcomp=False):
        wlength,x,y = self.parse(w)
        if not self.has([wlength,x,y]):
            return
        info = self.db[wlength]['x'][x][y]['info']
        rinfo = []
        if flg_revcomp:
            rw = self.convertor.revcomplement([wlength,x,y])
            if rw != [wlength,x,y]:
                wlength,x,y = rw
                rinfo = self.db[wlength]['x'][x][y]['info']
                for item in rinfo:
                    item[3] = -item[3]
        return info + rinfo

    def get_allWords(self,wl=None):
        output = []
        for wlength in self.db:
            if wl and wl != wlength:
                continue
            for x in self.db[wlength]['x']:
                for y in self.db[wlength]['x'][x]:
                    output.append([wlength,x,y])
        return output

    def delete(self,w,how_many=1):
        wlength,x,y = self.parse(w)
        if wlength not in self.db:
            return None
        if x not in self.db[wlength]['x']:
            return None
        if y not in self.db[wlength]['x'][x]:
            return None
        self.db[wlength]['x'][x][y]['counter'] -= how_many
        val_to_return = self.db[wlength]['x'][x][y]['counter']
        if val_to_return <= 0:
            val_to_return = 0
            del self.db[wlength]['x'][x][y]
            try:
                del self.db[wlength]['y'][y][x]
            except:
                pass
        if x in self.db[wlength]['x'] and not self.db[wlength]['x'][x]:
            del self.db[wlength]['x'][x]
        if y in self.db[wlength]['y'] and not self.db[wlength]['y'][y]:
            del self.db[wlength]['y'][y]
        if wlength in self.db and not self.db[wlength]:
            del self.db[wlength]
        return val_to_return
    
    def k_mers(self,k):
        words = self.get_allWords()
        k_db = db()
        for word in words:
            subwords = self.convertor.constituents(word,k)
            if not subwords:
                continue
            for i in range(len(list(subwords))):
                k_db.add(subwords[i])
                k_db.append_info(subwords[i],word+[i])
        return k_db

    #### INFO
    def has(self,w,flg_complement=False):
        wlength,x,y = self.parse(w)
        try:
            count = self.db[wlength]['x'][x][y]
            return True
        except:
            if flg_complement:
                wlength,x,y = self.convertor.revcomplement([wlength,x,y])
                return self.has(wlength,x,y)
            else:
                return False
        
    def count_entries(self):
        count = 0
        for wlength in self.db:
            for x in self.db[wlength]['x']:
                count += len(self.db[wlength]['x'][x])
        return count

    def count_entries_perX(self,wlength,x):
        return len(self.db[wlength]['x'][x])

    def count_entries_perY(self,wlength,y):
        return len(self.db[wlength]['y'][y])
    
    def size(self):
        count = 0
        for wlength in self.db:
            for x in self.db[wlength]['x']:
                for y in self.db[wlength]['x'][x]:
                    count += 1
        return count

    def get_topX(self,wlength):
        count = 0
        topX = 0
        for x in self.db[wlength]['x']:
            val = self.count_entries_perX(wlength,x)
            if val > count:
                count = val
                topX = x
        return topX,count

    def get_topY(self,wlength):
        count = 0
        topY = 0
        for y in self.db[wlength]['y']:
            val = self.count_entries_perY(wlength,y)
            if val > count:
                count = val
                topY = y
        return topY,count

    def get_bottomX(self,wlength):
        count = "0"
        bottomX = 0
        for x in self.db[wlength]['x']:
            val = self.count_entries_perX(wlength,x)
            if val < count:
                count = val
                bottomX = x
        return bottomX,count

    def get_bottomY(self,wlength):
        count = "0"
        bottomY = 0
        for y in self.db[wlength]['y']:
            val = self.count_entries_perY(wlength,y)
            if val < count:
                count = val
                bottomY = y
        return bottomY,count
    
    def get_mode(self,wl=0,flg_revcomp=False):
        C = 0
        X = Y = 0
        if wl and wl not in self.db:
            return None
        for wlength in self.db:
            if wl and wlength != wl:
                continue
            for x in self.db[wlength]['x']:
                for y in self.db[wlength]['x'][x]:
                    count = self.get([wlength,x,y])
                    if flg_revcomp:
                        rw = self.convertor.revcomplement([wlength,x,y])
                        if rw != [wlength,x,y]:
                            count = self.get(rw)
                    if count > C:
                        c = count
                        X = x
                        Y = y
        return [wlength,x,y],count
    
    def get_stat(self,wl=0):
        output = []
        word_range = list(self.db.keys())
        word_range.sort()
        for wlength in word_range:
            if wl and wlength != wl:
                continue
            count = 0
            for x in self.db[wlength]['x']:
                for y in self.db[wlength]['x'][x]:
                    count += self.get([wlength,x,y])
            output.append(str(wlength)+"-mers\t"+str(count))
        return "\n".join(output)
###############################################################################
def format_number(num,dig,zoom=0):
    try:
        num = float(num)
    except:
        return num
    return str(int((10**(dig+zoom))*num)/float(10**dig))

def abbreviate(name,names,n=10):
    if len(list(name)) <= n and name not in names:
        return name
    if len(list(name)) > n:
        key = name[:n]
    i = 1
    while key in names:
        key = key[:n-len(list(str(i)))-1]+"_"+str(i)
        i += 1
    return key

def remove_duplicates(table):
    if len(table) < 2:
        return table
    table.sort()
    for i in range(len(table)-1,0,-1):
        if table[i] == table[i-1]:
            del table[i]
    return table

# pattern_type = n1_4mer
def parse_pattern_type(pattern_type):
    p_type = pattern_type[0]
    norm = int(pattern_type[1:pattern_type.find("_")])
    wl = int(pattern_type[pattern_type.find("_")+1:pattern_type.find("mer")])
    return p_type,norm,wl

def get_EuclidianDistance(first_coordinates,second_coordinates):
    if len(first_coordinates) != len(second_coordinates):
        return None
    sum = 0
    for j in range(len(first_coordinates)):
        sum += (first_coordinates[j]-second_coordinates[j])**2
    return math.sqrt(sum)

def get_residualShift(first_confidence=1.0,second_confidence=1.0,D12=1.0):
    if not first_confidence and not second_confidence:
        return D12
    return 0.2866*math.log(math.sqrt(2/(first_confidence**2+second_confidence**2)))*D12
    
def get_distanceByCoordinates(first_coordinates,second_coordinates,first_confidence=1.0,second_confidence=1.0,D12=1.0):
    distance = get_EuclidianDistance(first_coordinates,second_coordinates)
    if distance==None:
        return None
    # correction for the confidence leads that database entries with low confidence are never identified 
    # even in self identification test
    return D12*(distance + get_residualShift(first_confidence,second_confidence))
        
# Get list of values and return [mean value,standard deviation]
def getMeanAndStDev(DataList,ind=2,binome=None):
    if not DataList:
        return [0,0,0]
    if len(DataList) == 1:
        return [DataList[0],0,0]
    sum = 0.0
    sum_sq = 0.0
    sum_cube = 0.0
    # if there is a number of lists, recalculate stdDev for all lists and return the average and the skew for the last list
    if type(DataList[0]) == type([1,2]):
        dm = len(DataList)
        N = len(DataList[0])
    else:
        dm = 1
        N = len(DataList)
        DataList = [DataList,]
    StdDev = 0
    for k in range(dm):
        for i in range(N):
            val = float(DataList[k][i])
            sum = sum + val
            if ind > 1:
                sum_sq = sum_sq + val*val
            if ind > 2:
                sum_cube = sum_cube + val*val*val
        N = float(N)
        mean = sum/N
        if ind > 1:
            if binome:
                if mean > 1:
                    bmean = mean/100.0
                    stdev = math.sqrt(bmean*(1-bmean))*100.0
                else:
                    stdev = math.sqrt(mean*(1-mean))
            else:
                var = sum_sq/N - mean*mean
                if var>0:
                    stdev = math.sqrt(var)
                else:
                    stdev = 0
        else:
            stdev = None
        if ind > 2:
            if binome:
                skew = 0
            else:
                if stdev:
                    skew = (sum_cube - 3.0*sum*sum_sq/N + 2.0*sum*sum*sum/N/N)/stdev/stdev/stdev/N
                else:
                    skew = 0
        else:
            skew = None
        if stdev:
            StdDev = math.sqrt(stdev*stdev + StdDev*StdDev)
        else:
            StdDev = 0
    return [mean,StdDev,skew]
