import math
import numpy as np
from numpy import exp
from lmfit import Model

#Defining Verhulst Model (Logistic)
def verhulst(x, mut, k):
    return (((2*k)/(1+exp(-mut*x))) - k)
def verhulst2(x, mut, k, s, a1, b1):
    return ((((2*k)/(1+exp(-mut*x+s))) - k) + a1*(b1 + ((2*k)/(1+exp(-mut*x+s))) - k))
def verhulst3(x, mut, k, s, a1, b1, a2, b2):
    return ((((2*k)/(1+exp(-mut*x+s))) - k) + a1*(b1 + ((2*k)/(1+exp(-mut*x+s))) - k) + a2*(b2 + ((2*k)/(1+exp(-mut*x+s))) - k))

class PlotGraph(object):
    def __init__(self,width=675,height=450,border=40,legend=125):
        self.datalist = []
        self.width = width
        self.height = height
        self.top_border = border
        self.left_border = border
        self.right_border = width-legend
        self.bottom_border = height-border
        self.font_size = 14

    def svg(self,vmodel):
        return "\n".join(["<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % 
                    (self.width,self.height,self.width,self.height),self._draw_plot(vmodel),"</svg>"])
                
    def _draw_plot(self,vmodel):
        # add axes
        maxX = np.around(max(vmodel.x)+0.05,decimals=-1)
        maxY = np.around(max(vmodel.y)+0.05,decimals=1)
        svg = ["<rect width=\"1\" height=\"1\" style=\"fill:wight;stroke-width:0;stroke-opacity:0.0\" />"]
        svg.append(self._svg_cartesian(self.left_border,self.top_border,self.right_border,
                    self.bottom_border,maxX,maxY,self.left_border,self.top_border))         
        x_origin = 2*self.left_border
        y_origin = self.top_border+10
        width = self.right_border
        height = self.bottom_border
        x_axis = width - (x_origin+10)
        y_axis = self.bottom_border - (y_origin+10)
        for i in range(len(self.datalist)):
            x = self.x[i]
            y = self.y[i]
            color = "black"
            title = "%s | %s" % (self.datalist[i][0],self.datalist[i][1])
            svg.append(self._svg_spot(x,x_origin,y,self.bottom_border,maxX,maxY,x_axis,y_axis,color,title))
        color = ["red","blue","green"]
        output = []
        if vmodel.minchisqr == 0: 
            output.append(round(vmodel.result1.chisqr,3))
            if output[0] < 0.01:
                msg = "VG"
            elif 0.01 <= output[0] < 0.1:
                msg = "Good"
            elif 0.05 <= output[0] < 0.5:
                msg = "Mod"
            elif 0.5 <= output[0] < 1:
                msg = "Bad"
            else: 
                msg = "VB"
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,210,"ChiSq: "+ msg)) 
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,225,"G: "+ str(round(vmodel.result1.best_values['mut'],2))))         
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,240,"K: "+ str(round(vmodel.result1.best_values['k'],2))))      
            svg.append(self._verhulst_graph(vmodel.result1,maxX,maxY,color[0],width,y_axis,x_origin,self.bottom_border))
        if vmodel.minchisqr == 1:
            output.append(round(vmodel.result21.chisqr,3))
            output.append(round(vmodel.result22.chisqr,3))
            if output[0] < 0.01:
                msg1 = "VG"
            elif 0.01 <= output[0] < 0.1:
                msg1 = "Good"
            elif 0.05 <= output[0] < 0.5:
                msg1 = "Mod"
            elif 0.5 <= output[0] < 1:
                msg1 = "Bad"
            else: 
                msg1 = "VB"
            if output[1] < 0.01:
                msg2 = "VG"
            elif 0.01 <= output[1] < 0.1:
                msg2 = "Good"
            elif 0.05 <= output[1] < 0.5:
                msg2 = "Mod"
            elif 0.5 <= output[1] < 1:
                msg2 = "Bad"
            else: 
                msg2 = "VB"
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,195,"ChiSq(R): "+ msg1))
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,215,"ChiSq(B): "+ msg2))
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,235,"G: "+ str(round(vmodel.result22.best_values['mut'],2))))         
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,255,"K: "+ str(round(vmodel.result22.best_values['k'],2))))  
            svg.append(self._verhulst_graph(vmodel.result21,maxX,maxY,color[0],width,y_axis,x_origin,self.bottom_border))
            svg.append(self._verhulst_graph(vmodel.result22,maxX,maxY,color[1],width,y_axis,x_origin,self.bottom_border))
        if vmodel.minchisqr == 2:
            output.append(round(vmodel.result31.chisqr,3))
            output.append(round(vmodel.result32.chisqr,3))
            output.append(round(vmodel.result33.chisqr,3))
            if output[0] < 0.01:
                msg1 = "VG"
            elif 0.01 <= output[0] < 0.1:
                msg1 = "Good"
            elif 0.05 <= output[0] < 0.5:
                msg1 = "Mod"
            elif 0.5 <= output[0] < 1:
                msg1 = "Bad"
            else: 
                msg1 = "VB"
            if output[1] < 0.01:
                msg2 = "VG"
            elif 0.01 <= output[1] < 0.1:
                msg2 = "Good"
            elif 0.05 <= output[1] < 0.5:
                msg2 = "Mod"
            elif 0.5 <= output[1] < 1:
                msg2 = "Bad"
            else: 
                msg2 = "VB"
            if output[2] < 0.01:
                msg3 = "VG"
            elif 0.01 <= output[2] < 0.1:
                msg3 = "Good"
            elif 0.05 <= output[2] < 0.5:
                msg3 = "Mod"
            elif 0.5 <= output[2] < 1:
                msg3 = "Bad"
            else: 
                msg3 = "VB"
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,185,"ChiSq(R): "+ msg1))
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,205,"ChiSq(B): "+ msg2))
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,225,"ChiSq(G): "+ msg3))
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,245,"G: "+ str(round(vmodel.result33.best_values['mut'],2))))         
            svg.append("<text x=\"%f\" y=\"%f\">%s</text>" % (550,265,"K: "+ str(round(vmodel.result33.best_values['k'],2))))  
            svg.append(self._verhulst_graph(vmodel.result31,maxX,maxY,color[0],width,y_axis,x_origin,self.bottom_border))
            svg.append(self._verhulst_graph(vmodel.result32,maxX,maxY,color[1],width,y_axis,x_origin,self.bottom_border))
            svg.append(self._verhulst_graph(vmodel.result33,maxX,maxY,color[2],width,y_axis,x_origin,self.bottom_border))       
        return "/n".join(svg)
        
    def _svg_spot(self,X,x_origin,Y,y_origin,maxX,maxY,x_axis,y_axis,color,title="",opacity=1.0,r=3):
        svg = ""
        if title:
            svg = "<a title=\"%s\">\n" % title
        proportion_x = X/maxX
        proportion_y = Y/maxY
        svg += ("<circle style=\"stroke:%s;fill:%s;stroke-width:%f;opacity:%f\" cx=\"%f\" cy=\"%f\" r=\"%f\" />" % 
                (color,color,1.0,opacity,proportion_x*x_axis+x_origin,y_origin-proportion_y*y_axis,r))
        if title:
            svg += "\n</a>"
        return svg
    
    def _verhulst_graph(self,result,maxX,maxY,color,width,y_axis,x_origin,y_origin,length=50,stroke_width=2):
        def calculate(result,x):
            k = float(result.best_values['k'])
            g = float(result.best_values['mut'])
            if len(result.best_values) == 2:
                return (((2*k)/(1+exp(-g*x))) - k)
            if len(result.best_values) == 5:
                s = float(result.best_values['s'])
                a1 = float(result.best_values['a1'])
                b1 = float(result.best_values['b1'])
                return ((((2*k)/(1+exp(-g*x+s))) - k) + a1*(b1 + ((2*k)/(1+exp(-g*x+s))) - k))
            if len(result.best_values) == 7:
                s = float(result.best_values['s'])
                a1 = float(result.best_values['a1'])
                b1 = float(result.best_values['b1'])
                a2 = float(result.best_values['a2'])
                b2 = float(result.best_values['b2'])
                return ((((2*k)/(1+exp(-g*x+s))) - k) + a1*(b1 + ((2*k)/(1+exp(-g*x+s))) - k) + a2*(b2 + ((2*k)/(1+exp(-g*x+s))) - k))
        svg = "<path style=\"stroke:%s;stroke-width:%f;fill:%s\" d=\"M%f,%f " % (color,stroke_width,"none",x_origin,y_origin)
        r = float(maxX)/length
        s = float(width-x_origin)/length
        for i in range(1,length):
            x = x_origin+i*s
            proportion = calculate(result,i*r)/maxY
            y = y_origin-proportion*y_axis
            if  proportion > 1:
                continue
            svg += "L%f,%f" % (x,y)
        svg += "\" />"
        return svg
          
    def _svg_cartesian(self,X,Y,width,height,maxX,maxY,x_indent=30,y_indent=30,flg_finish=False):
        svg = []
        # vertical axis
        svg.append(("<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" fill=\"none\" stroke=\"black\" stroke-width=\"%f\" />" %
             (X+x_indent,Y+10,X+x_indent,height,3.0)))
        # vertical arrow
        svg.append(self._svg_arrow(X+x_indent,Y+10,math.pi/2.0,2.0))
        # horizontal axis
        svg.append(("<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" fill=\"none\" stroke=\"black\" stroke-width=\"%f\" />" %
            (X+x_indent,height,width,height,3.0)))
        # horizontal arrow
        svg.append(self._svg_arrow(width,height,0,2.0))
        for i in range(1,11):
            xrange = (width - (X+x_indent+10))/10
            yrange = (height - (Y+20))/10
            svg.append(("<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" fill=\"none\" stroke=\"black\" stroke-width=\"%f\" />" %
             (X+x_indent+i*xrange,height,X+x_indent+i*xrange,height+8,3.0)))
            svg.append(("<text x=\"%f\" y=\"%f\" dx=\"%f\" dy=\"%f\" >%s</text>") % (X+x_indent+i*xrange,height,-20,+25,str(round(i*maxX/10,1))))
            svg.append(("<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" fill=\"none\" stroke=\"black\" stroke-width=\"%f\" />" %
             (X+x_indent-8,height-i*yrange,X+x_indent,height-i*yrange,3.0)))
            svg.append(("<text x=\"%f\" y=\"%f\" dx=\"%f\" dy=\"%f\" >%s</text>") % (X+x_indent,height-i*yrange,-45,+5,str(round(i*maxY/10,2)))) 	
        if flg_finish:
            svg.insert(0,"<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewbox=\"0 0 %d %d\">" % (width,height,width,height))
            svg.append("</svg>")
        return "\n".join(svg)


    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:black;stroke-width:%f;fill:black\" d=\"M%f,%fL%f,%fL%f,%fL%f,%fZ\" />" %
                    (1.0,ax,ay,bx,by,cx,cy,dx,dy))
    
        
#Verhulst Class object
class Verhulst(PlotGraph):
    
    #Initializing OUP, WGS and Specie name
    def __init__(self,datalist,cluster):
        PlotGraph.__init__(self)
        try:
            self.x = [float(item[2]) for item in datalist]
            self.y = [float(item[3]) for item in datalist]
        except:
            try:
                datalist = datalist[1:]
                self.x = [float(item[2]) for item in datalist]
                self.y = [float(item[3]) for item in datalist]
            except:
                print()
                print("Datalist is in wrong format!")
                return
        self.x21 = []
        self.y21 = []
        self.x22 = []
        self.y22 = []
        self.x31 = []
        self.y31 = []
        self.x32 = []
        self.y32 = []
        self.x33 = []
        self.y33 = []
        self.datalist = datalist
        self.datalist21 = []
        self.datalist22 = []
        self.datalist31 = []
        self.datalist32 = []
        self.datalist33 = []
        self.header = [' ']
        self.verhulst_model = []
        self.data_matrix = []
        self.result1 = ''
        self.result21 = ''
        self.result22 = ''
        self.result31 = ''
        self.result32 = ''
        self.result33 = ''
        self.minchisqr = 0
        self.cluster = cluster
        self.g_param = 0
        self.k_param = 0

    #Creating Specie Matrix of Zeros
    def set_matrix(self):
        self.header.append(self.datalist[0][0])
        for row in self.datalist:
            counter = 0
            for specie in self.header:
                if specie == row[1]:
                    counter += 1
            if counter == 0:
                self.header.append(row[1])

        self.data_matrix = [self.header]

        for column in range(1,len(self.header)):
            insertrow = [self.header[column]] + [0]*(len(self.header)-1)
            self.data_matrix.append(insertrow)

    #Assigning default parameters:
    def assign_default(self,g,k):
        self.g_param = g
        self.k_param = k

    #Calculating Verhulst Model Fitting based on Shape Method:
    #Using Verhulst Model Shape to divide data into clusters used for Model Fitting	
    def verhulst_fitting(self):
        median = sum(self.y) / float(len(self.y))
        
        gmod = Model(verhulst)
        self.result1 = gmod.fit(self.y, x=self.x, mut=1, k=median)
        
        init_mut = float(self.result1.best_values['mut'])	
        init_k = float(self.result1.best_values['k'])
        
        for combinations in self.datalist:
            emp_y = (((2*init_k)/(1+math.exp(-init_mut*float(combinations[2])))) - init_k)
            if float(combinations[3]) <= emp_y:
                self.datalist21.append(combinations)		
            else:
                self.datalist22.append(combinations)
                        
        for specie in self.datalist21:
            self.x21.append(float(specie[2]))
            self.y21.append(float(specie[3]))

        median1 = sum(self.y21) / float(len(self.y21))
        gmod1 = Model(verhulst)
        self.result21 = gmod1.fit(self.y21, x=self.x21, mut=1, k=median1)

        init_mut2 = float(self.result21.best_values['mut'])	
        init_k2 = float(self.result21.best_values['k'])

        for y in self.datalist21:
            emp_y = (((2*init_k2)/(1+exp(-init_mut2*float(y[2])))) - init_k2)
            abs_dif = abs(emp_y-float(y[3]))
            y.append(abs_dif)			

        counter = 0 
        self.datalist21.sort(key = lambda row: (row[3]))
        for distance in self.datalist21:
            counter += 1
            distance.append(counter)

        counter = 0 
        self.datalist21.sort(key = lambda row: (row[3]))
        for position in self.datalist21:
            counter += 1
            position.append(counter)

        for x in self.datalist21:
            average_sort = (x[-2] + x[-1])/2
            x.append(average_sort)
	
        self.datalist21.sort(key = lambda row: row[-1])
        self.datalist22.append(self.datalist21[0])

        for specie in self.datalist22:
            self.x22.append(float(specie[2]))
            self.y22.append(float(specie[3]))
	
        self.datalist22.sort(key = lambda row: float(row[2]))
        min_x = float(self.datalist22[0][2])

        self.datalist22.sort(key = lambda row: float(row[2]))
        min_x = float(self.datalist22[0][2])
        min_y = float(self.datalist22[0][3])

        for y in self.datalist21:
            if float(y[2]) <= min_x and float(y[3]) <= min_y:
                self.datalist22.append(y)

        for specie in self.datalist22:
            self.x22.append(float(specie[2]))
            self.y22.append(float(specie[3]))

        median2 = sum(self.y22) / float(len(self.y22))
        gmod2 = Model(verhulst2)

        pars = gmod2.make_params(mut=init_mut2, k=init_k2, s=1, a1=0.001, b1=100)
        pars['mut'].vary = False
        pars['k'].vary = False
        pars['a1'].min = 0
        pars['s'].min = 1

        self.result22 = gmod2.fit(self.y22, x=self.x22, params=pars)

        init_mut22 = float(self.result22.best_values['mut'])	
        init_k22 = float(self.result22.best_values['k'])
        init_a1 = float(self.result22.best_values['a1'])	
        init_b1 = float(self.result22.best_values['b1'])
        init_s = float(self.result22.best_values['s'])

        for combinations in self.datalist:
            emp_y1 = (((2*init_k22)/(1+exp(-init_mut22*float(combinations[2])+init_s))) - init_k22) + init_a1*(init_b1 + ((2*init_k22)/(1+exp(-init_mut22*float(combinations[2])+init_s))) - init_k22)
            emp_y2 = (((2*init_k22)/(1+exp(-init_mut22*float(combinations[2])))) - init_k22)
            if float(combinations[3]) <= emp_y2:	
                self.datalist31.append(combinations)
            elif emp_y2 <= float(combinations[3]) <= emp_y1:
                self.datalist32.append(combinations)
            else:
                self.datalist33.append(combinations)

        for specie in self.datalist31:
            self.x31.append(float(specie[2]))
            self.y31.append(float(specie[3]))
        
        #Checking for overlapping Cluster points
        self.datalist32.sort(key = lambda row: float(row[2]))
        min_x = float(self.datalist32[0][2])
        min_y = float(self.datalist32[0][3])	

        for y in self.datalist31:
            if float(y[2]) <= min_x and float(y[3]) <= min_y:
                self.datalist32.append(y)	

        for combinations in self.datalist31:
            emp_y = (((2*init_k22)/(1+exp(-init_mut22*float(combinations[2])))) - init_k22)
            abs_dif = abs(emp_y-float(combinations[3]))
            combinations.append(abs_dif)

        counter = 0 
        self.datalist31.sort(key = lambda row: (row[3]))
        for distance in self.datalist31:
            counter += 1
            distance.append(counter)

        counter = 0 
        self.datalist31.sort(key = lambda row: (row[3]))
        for position in self.datalist31:
            counter += 1
            position.append(counter)

        for x in self.datalist31:
            average_sort = (x[-2] + x[-1])/2
            x.append(average_sort)

        self.datalist31.sort(key = lambda row: row[-1])
        self.datalist32.append(self.datalist21[0])

        for specie in self.datalist32:
            self.x32.append(float(specie[2]))
            self.y32.append(float(specie[3]))

        for combinations in self.datalist32:
            emp_y = (((2*init_k22)/(1+exp(-init_mut22*float(combinations[2])))) - init_k22) + init_a1*(init_b1 + ((2*init_k22)/(1+exp(-init_mut22*float(combinations[2])+init_s))) - init_k22)
            abs_dif = abs(emp_y-float(combinations[3]))
            combinations.append(abs_dif)

        counter = 0 
        self.datalist32.sort(key = lambda row: (row[3]))
        for distance in self.datalist32:
            counter += 1
            distance.append(counter)

        counter = 0 
        self.datalist32.sort(key = lambda row: (row[3]))
        for position in self.datalist32:
            counter += 1
            position.append(counter)

        for x in self.datalist32:
            average_sort = (x[-2] + x[-1])/2
            x.append(average_sort)
	
        self.datalist32.sort(key = lambda row: row[-1])
        self.datalist33.append(self.datalist21[0])

        self.datalist33.sort(key = lambda row: float(row[2]))
        min_x = float(self.datalist33[0][2])	
        min_y = float(self.datalist33[0][3])

        for y in self.datalist32:
            if float(y[2]) <= min_x and float(y[3]) <= min_y:
                self.datalist33.append(y)

        for specie in self.datalist33:
            self.x33.append(float(specie[2]))
            self.y33.append(float(specie[3]))

        median1 = sum(self.y31) / float(len(self.y31))
        median2 = sum(self.y32) / float(len(self.y32))
        median3 = sum(self.y33) / float(len(self.y33))

        gmod1 = Model(verhulst)
        gmod2 = Model(verhulst2)
        gmod3 = Model(verhulst3)

        self.result31 = gmod1.fit(self.y31, x=self.x31, mut=1, k=median1)
        init_mut3 = float(self.result31.best_values['mut'])	
        init_k3 = float(self.result31.best_values['k'])
        
        pars = gmod2.make_params(mut=init_mut3, k=init_k3, s=1, a1=0.01, b1=1)
        pars['mut'].vary = False
        pars['k'].vary = False
        pars['a1'].min = 0
        pars['s'].min = 1


        pars2 = gmod3.make_params(mut=init_mut3, k=init_k3, s=1, a1=0.001, b1=10, a2=0.002, b2=20)
        pars2['mut'].vary = False
        pars2['k'].vary = False
        pars2['a1'].min = 0
        pars2['a2'].min = 0
        pars2['s'].min = 1

        #Model Fitting
        self.result32 = gmod2.fit(self.y32, x=self.x32, params=pars)
        self.result33 = gmod3.fit(self.y33, x=self.x33, params=pars2)	
        
        #self.result32 = result32.fit_report()

    #Determining best Verhulst Model in terms of Chi-Square criterion and output Textfile
    def chisqr_criterion(self): 		
        model2chisqr = self.result21.chisqr + self.result22.chisqr
        model3chisqr = self.result31.chisqr + self.result32.chisqr + self.result33.chisqr
        chisqr = [self.result1.chisqr,model2chisqr,model3chisqr]
        chisqr_dif = [model2chisqr-self.result1.chisqr,model3chisqr-model2chisqr]
        #if chisqr_dif[1] < 0.05:
        #    self.minchisqr = 0
        if chisqr_dif[1] < 0.025:
            self.minchisqr = 1
        else:
            self.minchisqr = chisqr.index(min(chisqr))
        if self.cluster != 0:
            self.minchisqr = int(self.cluster)-1	
        if self.minchisqr == 0:
            self.verhulst_model.append(self.result1)
        elif self.minchisqr == 1:
            self.verhulst_model.append(self.result21)
            self.verhulst_model.append(self.result22)
        elif self.minchisqr == 2:
            self.verhulst_model.append(self.result31)
            self.verhulst_model.append(self.result32)
            self.verhulst_model.append(self.result33)
            
    #Converting Distance between species into clusters number based on Verhulst Model and level of Jump into a Matrix
    def fill_matrix(self):
        for combination in self.datalist21:
            for combination2 in self.datalist22:
                if combination[0] == combination2[0] and combination[1] == combination2[1]:
                    self.datalist22.remove(combination2)
            
        for combination in self.datalist32:
            for combination2 in self.datalist33:
                if combination[0] == combination2[0] and combination[1] == combination2[1]:
                    self.datalist33.remove(combination2)
            
        for combination in self.datalist31:
            for combination2 in self.datalist32:
                if combination[0] == combination2[0] and combination[1] == combination2[1]:
                    self.datalist32.remove(combination2)
        
        if self.minchisqr == 1:
                for combination in self.datalist22:
                    rowpos = self.header.index(combination[1])
                    columnpos = self.header.index(combination[0])
                    self.data_matrix[columnpos][rowpos] = 1
                    self.data_matrix[rowpos][columnpos] = 1

                for combination in self.datalist21:
                    rowpos = self.header.index(combination[1])
                    columnpos = self.header.index(combination[0])
                    self.data_matrix[columnpos][rowpos] = 0
                    self.data_matrix[rowpos][columnpos] = 0
        
        if self.minchisqr == 2:
                for combination in self.datalist33:
                    rowpos = self.header.index(combination[1])
                    columnpos = self.header.index(combination[0])
                    self.data_matrix[columnpos][rowpos] = 2
                    self.data_matrix[rowpos][columnpos] = 2

                for combination in self.datalist32:
                    rowpos = self.header.index(combination[1])
                    columnpos = self.header.index(combination[0])
                    self.data_matrix[columnpos][rowpos] = 1
                    self.data_matrix[rowpos][columnpos] = 1

                for combination in self.datalist31:
                    rowpos = self.header.index(combination[1])
                    columnpos = self.header.index(combination[0])
                    self.data_matrix[columnpos][rowpos] = 0
                    self.data_matrix[rowpos][columnpos] = 0


if __name__ == "__main__":
        #Example		
        input_fname = "example.txt"
        data = open(input_fname)
        x = []
        y = []
        datalist = []

        for line in data:
                tmp = line.replace("\n","").split("\t")
                datalist.append(tmp)
                x.append(float(tmp[2]))
                y.append(float(tmp[3]))

        enterobacteria = Verhulst(datalist,y,x)
        enterobacteria.set_matrix()
        enterobacteria.verhulst_fitting()
        enterobacteria.chisqr_criterion()
        enterobacteria.fill_matrix()
        print((enterobacteria.data_matrix))

        #LMFIT Module functions
        #.chisqr() - Get Chi-Square criterion value
        #.ci_out() - Get Confidence interval
        #.params() - Parameter of best fit
        #.resudual() - Sum of residual value 
        #.fit_report() - Verhulst Model fitting Summary
        print((enterobacteria.verhulst_model[0].fit_report()))

