usage = "\n"
usage += "-------------------------------------------------------------------------------------\n"
usage += "scattering.sage\n"
usage += "A sage code for computing scattering diagrams and broken lines\n"
usage += "by Tim Gräfnitz (timgraefnitz.com, tim.graefnitz@gmx.de)\n"
usage += "version: 07.03.2023\n"
usage += "-------------------------------------------------------------------------------------\n"
usage += "\n"
usage += "For interactive usage type 'scattering()'\n"
usage += "\n"
usage += "There are two classes: Diagram and Ray\n"
usage += "Initialize function ring QQ(x,y)[t_1,...t_r] by 'initialize(r)'\n"
usage += "Initialize ray by Ray(base,direction,function), e.g. 'Ray((-1,0),(1,0),1+t_0*x)'\n"
usage += "Initialize diagram by Diagram(ray_list), e.g. 'D = Diagram([((-1,0),(1,0)1+t_0*x),((0,-1),(0,1),1+t_1*y)])'\n"
usage += "Initialize predefined diagrams e.g. by 'D = Diagram(case='P2',order=3)' or 'D = Diagram('P2',3)'\n"
usage += "\n"
usage += "Compute scattering of D up to given order by D.scattering(order), e.g. 'D1 = D.scattering(3)'\n"
usage += "For predefined diagrams the computation may be optimized, so type D.scattering(order,case), e.g. 'D1 = D.scattering(3,'P2')\n"
usage += "You can print, draw and give tex code: 'print(D)', 'D.draw()' or 'D.tex()'\n"
usage += "print() and draw() take optional arguments: special:list(Ray), colors:boolean, function:boolean, directions:list(tuple)\n"
usage += "Example input: 'Diagram('P2',3).scattering(3,'P2',True).draw(colors=True,directions=[(1,0)])' draws the diagram of (P2,E) for degree 3.\n"
usage += "\n"
usage += "A broken line is encoded as a list of rays in reversed order; "
usage += "the base of the first ray is the endpoint, the direction of the last ray is the negative of the incoming direction\n"
usage += "Compute broken lines to given order ending in a point pt in a Diagram D by D.brokenlines(pt,order), e.g. 'D.brokenlines((10,0.123),2)'\n"
usage += "Again, for predefined diagrams this may be optimized, so type D.brokenlines(pt,order,case) instead, e.g. 'D.broken_lines((10,0.132),2,'P2')'\n"
usage += "You can give the optinal arguments end:list(tuple) and in:list(tuple) to specify incoming and ending directions.\n"
usage += "Example input: 'Diagram('P2',2).scattering(2,'P2').broken_lines((10,0.132),2,'P2',True,end=[(-1,0)],in=[(5,0)])' prints the broken lines of (P2,E)\n"
usage += "to compute the 2-marked log Gromov-Witten invariant R_1,5 of conics meeting an elliptic curve E in a fixed point with order 1 and another point with order 5.\n"
#print(usage)

import re
import numpy as np
import itertools
import time
import os.path
from sage.plot.line import Line
warnings.simplefilter(action='ignore', category=DeprecationWarning)

### INTERACTIVE USAGE METHODS ###

def scattering():
    """ Function for interactive calculation of scattering diagrams."""
    pre_own = input("predefined or own?  ")
    if pre_own == "predefined":
        pre_case = input("choose: del_pezzo, standard, primitive  ")
        if pre_case == "del_pezzo":
            case = input("choose: P2, P1xP1, dP{k}  ")
        order = int(input("calulcation to order? "))
        print_bool = input("print result? [y/n]  ")
        draw_bool = input("draw result? [y/n]  ")
        tex_bool = input("print tex code? [y/n]  ")
        D = Diagram(case,order)
        D1 = D.scattering(order,case)
        if print_bool:
            print(D1)
        if draw_bool:
            D1.draw()
        if tex_bool:
            D1.tex()
    elif pre_own == "own":
        num_vars = int(input("number of t-variables?  "))
        initialize(num_vars)
        ray_data = input("input diagram of the form [((-1,0),(1,0)1+t_0*x),((0,-1),(0,1),1+t_1*y)]  ")
        order = int(input("calulcation to order?  "))
        print_bool = input("print result? [y/n]  ")
        draw_bool = input("draw result? [y/n]  ")
        tex_bool = input("print tex code? [y/n]  ")
        D = Diagram(ray_data)
        D1 = D.scattering(order)
        if print_bool == "y":
            print(D1)
        if draw_bool == "y":
            D1.draw()
        if tex_bool == "y":
            D1.tex()

""" Initializes the ring for wall functions: Laurent series in x,y over the polynomial ring in t_1,...,t_r """
def initialize(number_of_t_vars:int):
    global x,y,h,t,Q,R
    x,y,h = var('x,y,h')
    t = [var('t_%d' % i) for i in range(number_of_t_vars)]
    P = PowerSeriesRing(QQ,'h',order='negdeglex')
    Q = LaurentPolynomialRing(P,'x,y',order='negdeglex')
    R = PolynomialRing(Q,t,order='negdeglex')
    
### DIAGRAM CLASS ###

"""Class of wall structures/scattering diagrams"""

class Diagram():
    def __init__ (self, rays_or_case = [], monodromy_rays_or_order = [], path = "calculations.txt", max_trafos = 'infty'):
        # if input is list of rays
        if isinstance(rays_or_case,list) or isinstance(rays_or_case,set):
            self.rays = [ray if isinstance(ray,Ray) else Ray(*ray) for ray in rays_or_case]
            self.monodromy_rays = [ray if isinstance(ray,Ray) else Ray(*ray) for ray in monodromy_rays_or_order]
            self.case = ''
            self.order = 'infty'
        # if input is case
        else:
            self.case = rays_or_case
            self.rays = []
            self.monodromy_rays = []
            self.order = 'infty'
        self.diagrams = [self]
        self.broken_lines = {}
        self.max_trafos = max_trafos
        # Initialize diagram if case is given
        if not self.case == '':
            if self.case == 'std':
                initialize(3)
                m,n = monodromy_rays_or_order
                self.rays = [Ray((-1,0),(1,0),(1+t_1*x)^m),Ray((0,-1),(0,1),(1+t_2*y)^n)]
                self.diagrams = self.__get_diagrams(self.case,self.order,path)
                self.case = ''
            elif self.case == 'exp':
                initialize(3)
                m,n = monodromy_rays_or_order
                self.rays = [Ray((-1,0),(1,0),1+t_1*x^m),Ray((0,-1),(0,1),1+t_2*y^n)]
                self.diagrams = self.__get_diagrams(self.case,self.order,path)
                self.case = ''
            elif self.case == 'det':
                initialize(3)
                m = monodromy_rays_or_order
                self.rays = [Ray((-1,0),(1,0),1+t_1*x),Ray((1,-m),(-1,m),1+t_2*x^(-1)*y^m)]
                self.diagrams = self.__get_diagrams(self.case,self.order,path)
                self.case = ''
            else: # del Pezzos
                self.order = monodromy_rays_or_order
                initialize(4*len(self.__kinks(self.case))*self.order)
                self.rays = self.__get_rays(self.case,self.order)
                self.monodromy_rays = self.__get_monodromy_rays(self.case,self.order)
                self.diagrams = self.__get_diagrams(self.case,self.order,path)
                self.order = self.__order(self.case)*self.order
                self.classes = self.__class_names(self.case)
        # transform all rays according to mondromoy rays
        for ray in self.rays:
            self.__transform(ray)
    
### PUBLIC METHODS ###

    def class_rays(self, case, order):
        result = []
        kinks = self.__kinks(case)
        classes = self.__classes(case)
        lengths = self.__lengths(case)
        base_top = (0,1-lengths[0])
        base_bottom = (0,0+lengths[-1])
        slope_top = 0
        slope_bottom = 0
        for i in range(order):
            for j in range(len(kinks)):
                slope_top += kinks[j]
                slope_bottom += kinks[-1-j]
                base_top = np.add(base_top,(-slope_top,lengths[j]))
                base_bottom = np.add(base_bottom,(-slope_bottom,-lengths[-1-j]))
                result.append(Ray(base_top,(1,0),classes[j]))
                result.append(Ray(base_bottom,(1,0),classes[-1-j]))
        return result

    """Calculates the diagrams consistent up to the given order and returns the highest one"""
    def scattering(self, order:int, case:str = ''):
        start = time.time()
        ind_start = 1
        if len(self.diagrams) > 1:
            print("Already calculated to order", len(self.diagrams)-1)
            ind_start = len(self.diagrams)
        if ind_start <= order:
            print("scattering:")
            if len(self.__intersections()) == 1:
                for i in range(ind_start,order+1):
                    diag = self.diagrams[i-1].__scatter_step_one_vertex(i,order,case)
                    self.diagrams.append(diag)
            else:
                for i in range(ind_start,order+1):
                    diag = self.diagrams[i-1].__scatter_step(i,order,case)
                    self.diagrams.append(diag)
            print("time elapsed: " + str(time.time() - start) + "s")
        for ray in self.diagrams[order].rays:
            self.diagrams[order].__calc_trop(ray)
        if self.case != '':
            self.diagrams[order].case = self.case
            self.diagrams[order].__calc_classes(order)
        return self.diagrams[order]

    """Scattering for primitive diagrams"""
    def scattering_prim(self, order = float('inf')):
        start = time.time()
        result = self.rays
        new_rays = []
        for ray1,ray2 in set(itertools.combinations(self.rays,2)):
            new_ray = self.__scatter_prim(ray1,ray2)
            if isinstance(new_ray,Ray) and new_ray.t_order <= order:
                new_rays.append(new_ray)
        while len(new_rays) > 0:
            new_rays2 = []
            for ray1 in result:
                for ray2 in new_rays:
                    new_ray = self.__scatter_prim(ray1,ray2)
                    if isinstance(new_ray,Ray) and new_ray.t_order <= order:
                        new_rays2.append(new_ray)
            for ray1,ray2 in set(itertools.combinations(new_rays,2)):
                new_ray = self.__scatter_prim(ray1,ray2)
                if isinstance(new_ray,Ray) and new_ray.t_order <= order:
                    new_rays2.append(new_ray)
            result += new_rays
            new_rays = new_rays2 
        print("time elapsed: " + str(time.time() - start) + "s")
        return Diagram(result)

    """Gives the broken lines up to a given order with endpoint pt"""
    def brokenlines(self, pt:tuple, order:int, case:str = '', exp_end = 'all', exp_in = 'all'):
        start = time.time()
        self.broken_lines[pt] = []
        print("broken lines:")
        x_order = order * self.__order(case) if case != '' else order
        if exp_end == 'all':
            exp_list = [(a,b) for a in range(-x_order,x_order+1) for b in range(-x_order,x_order+1) if (a,b) != (0,0)]
        elif exp_end == 'infty':
            exp_list = [(a,0) for a in range(-x_order+1,0)]
        else:
            exp_list = exp_end
        for exponent in exp_list:
            if exponent == (0,0):
                continue
            print("ending exponent:",exponent)
            direction = (exponent[0]/gcd(exponent[0],exponent[1]),exponent[1]/gcd(exponent[0],exponent[1]))
            ray_incoming = Ray(pt,direction,x^exponent[0]*y^exponent[1])
            incoming_order = -exponent[0] if case != '' else 0
            self.__transform(ray_incoming)
            result = []
            self.__broken_iter(ray_incoming,[ray_incoming],result,order,case,exp_in)
            if self.case != '':
                for broken_line in result:
                    self.__calc_class_broken(broken_line)
            self.broken_lines[pt] += result
        print("\ntime elapsed: " + str(time.time() - start) + "s")
        return self

    """ Contributed by Eric Zaslow """
    def wallcrossing(self, ray, order):
        int_dict = self.__int_dict(ray)
        int_pts = sorted(int_dict.keys(),key=lambda pt: np.linalg.norm(np.subtract(pt,ray.base)))
        for int_pt in int_pts:
            for r in int_dict[int_pt]:
                W = ray.function
                f = r.function
                n = r.direction
                ray.function = expand(self.__cross(W,f,n,order))
                print(int_pt,r)
                print(self.__coeffs(ray.function))
        print(self.__coeffs(ray.function))
    
    def print_list(self):
        output = ""
        output += "rays:\n" + str([(ray.base,ray.direction,ray.function) for ray in self.rays])
        output += "\nmonodromy rays:\n" + str([(ray.base,ray.direction,ray.function) for ray in self.monodromy_rays])
        output += "\nbroken lines:\n" + str([(ray.base,ray.direction,ray.function) for pt, line_list in self.broken_lines.items() for line in line_list for ray in line])
        return output

    def save(self, path = "calculations.txt"):
        """Writes the output of print(self) to the given path (path to a .txt file)"""
        with open(path, 'w') as f:
            f.write(str(self.print_list()))
    
    def draw(self, special = [], colors = False, functions = False, directions = []):
        """ Displays a diagram as png image. """
        color_list = ['red','black','blue','green','brown','gray']
        special = special if isinstance(special,Diagram) else Diagram(special)
        G = Graphics()
        G.axes(show=None)
        for ray in self.__diag_conv(special,colors,functions,directions):
            if ray.special:
                color = color_list[ray.order % 6]
            else:
                color = Color(255,floor(255*(1-1/1.2^ray.order)),0)
            if ray.arrow:
                G += arrow(ray.base,ray.endpoint,thickness=1 if ray.special else 0.5,rgbcolor=color)
            else:
                G += line([ray.base,ray.endpoint],thickness=1 if ray.special else 0.5,rgbcolor=color)
        for ray in Diagram(self.monodromy_rays).__diag_conv():
            G += line([ray.base,ray.endpoint],linestyle='dashed')
        for pt, broken_line_list in self.broken_lines.items():
            G += circle((pt[0],pt[1]),0.005,fill=True,rgbcolor='green')
            for broken_line in broken_line_list:
                for i in range(len(broken_line)-1):
                    for part in broken_line[i].transforms:
                        G += line([(part.base[0],part.base[1]),(part.endpoint[0],part.endpoint[1])],thickness=2,color='green')
                for part in broken_line[-1].transforms:
                    G += line([(part.base[0],part.base[1]),(part.base[0]+part.direction[0],part.base[1]+part.direction[1])],thickness=2,color='green')
        G.show(figsize=[40,20])

    """ Gives LaTex code for a tikz picture of the diagram """
    def tex(self, special = [], colors = False, functions = False, directions = []):
        color_list = ['red','black','blue','green','brown','gray']
        special = special if isinstance(special,Diagram) else Diagram(special)
        output = "\\begin{figure}[h!]\n"
        output += "\\centering\n"
        output += "\\begin{tikzpicture}[scale=1"
        if colors:
            output += ",define rgb/.code={\\definecolor{mycolor}{RGB}{#1}}, rgb color/.style={define rgb={#1},mycolor}]"
        output += "]\n"
        for ray in Diagram(self.monodromy_rays).__diag_conv():
            output += "\\draw[dashed] " + str(tuple(map(lambda x: round(x,3),ray.base))) + " -- " + str(tuple(map(lambda x: round(x,3),ray.endpoint))) + ";\n"
        for ray in self.__diag_conv(special,colors,functions,directions):
            arrow = "->" if ray.arrow else "-"
            thick = ",thick" if ray.special else ""
            color = ""
            if colors:
                if ray.direction in directions:
                    color = ",color=" + color_list[ray.order % 6]
                else:
                    color = ",rgb color={255,"+str(floor(255*(1-1/1.2^ray.order)))+",0}"
            output += "\\draw["+arrow+thick+color+"] " + str(tuple(map(lambda x: round(x,3),ray.base)))
            output += " -- " + str(tuple(map(lambda x: round(x,3),ray.endpoint)))
            if functions and (directions == [] or ray.direction in directions):
                output += " node[" + ray.node_position + "]{$" + ray.function + "$}"
            output += ";\n"
        for pt, broken_line_list in self.broken_lines.items():
            for broken_line in broken_line_list:
                output += "\\fill[black!40!green] " + str(tuple(map(lambda x: round(x,3),pt))) + " circle (2pt);\n"
                for i in range(0,len(broken_line)-1):
                    for part in broken_line[i].transforms:
                        output += "\\draw[black!40!green] " + str(tuple(map(lambda x: round(x,3),part.base))) + " -- " + str(tuple(map(lambda x: round(x,3),part.endpoint))) + ";\n"
                for part in broken_line[-1].transforms:
                    output += "\\draw[black!40!green] " + str(tuple(map(lambda x: round(x,3),part.base))) + " -- " + str(tuple(map(lambda x: round(x,3),np.add(part.base,part.direction)))) + ";\n"
        output += "\\end{tikzpicture}\n"
        output += "\\caption{A scattering diagram.}\n"
        output += "\\end{figure}"
        print(output)

    ### STANDARD METHODS ###
        
    def __repr__ (self):
        self.rays = sorted(sorted(self.rays, key=lambda ray: ray.direction[1]), key=lambda ray: ray.direction[0])
        """Representation function for diagrams, used e.g. in print()"""
        output = "Diagram:" + '\n'.join(["Monodromy " + str(r) for r in self.monodromy_rays]) + "\n"
        output += '\n'.join([str(r) + self.__class_string(r) for r in self.rays]) + "\n"
        output += '\n'.join(["Broken Line: " + self.__class_string(broken_line[-1]) + "\n" + '\n'.join([str(part) for part in broken_line]) for pt, line_list in self.broken_lines.items()for broken_line in line_list])
        return output
        
    def __eq__ (self, other):
        """Equality for diagrams means equality of their rays"""
        return sorted(self.rays) == sorted(other.rays)
        
    def __add__ (self, other):
        """Add a ray to the diagram by putting it in the ray list"""
        self.rays = self.rays + other.rays
        return self
    
    def __contains__(self, ray):
        """Check whether a ray is contained in the ray list"""
        return ray if isinstance(ray,Ray) else Ray(*ray) in self.rays

    def __len__(self):
        """The length of the diagram is the length of its ray list"""
        return len(self.rays)

    def __iter__(self):
        """Iteration formula inherited from ray list"""
        return iter(self.rays)
    
    def add(self, ray):
        self.rays.append(ray if isinstance(ray,Ray) else Ray(*ray))
        return self
    
    def remove(self, ray):
        self.rays.remove(ray if isinstance(ray,Ray) else Ray(*ray))
        return self
    
    def update(self, other):
        return self + other
    
    ### PRIVATE METHODS - SCATTERING ###
    
    def __scatter_step(self, order, order_max, case):
        """Calculates scattering at a particular order"""
        print("order",order)
        result = deepcopy(self)
        new_rays = self.rays
        while True:
            local_diagrams = result.__local_diagrams(order_max,case,'all')
            # progress bar
            bar_width = len(local_diagrams.keys())
            sys.stdout.write("[%s]" % (" " * bar_width))
            sys.stdout.flush()
            sys.stdout.write("\b" * (bar_width+1))
            new_rays = []
            for pt, local_diagram in local_diagrams.items():
                scat = local_diagram.__scatter(order,order_max,case)
                local_diagrams[pt] = Diagram(local_diagram.rays + scat.rays)
                # for symmetric cases we only computed rays in a particular area, now we duplicate these rays
                if case != '':
                    scat = scat.__translate(order_max,case)
                new_rays = new_rays + scat.rays
                # update progess bar
                sys.stdout.write("-")
                sys.stdout.flush()
            sys.stdout.write("\n")          
            if len(new_rays) == 0:
                break
            result.rays += new_rays
            if case != '':
                break
            result.__unify(order_max)
        return result
    
    def __unify(self,order):
        result = Diagram()
        temp = deepcopy(self.rays)
        while True:
            if len(temp) == 0:
                break
            ray = temp[0]
            # combine functions
            function = ray.function
            remove = []
            for ray2 in temp:
                if all(np.isclose(ray.base,ray2.base)) and all(np.isclose(ray.direction,ray2.direction)):
                    function *= ray2.function
                    remove.append(ray2)
            for ray2 in remove:
                temp.remove(ray2)
            for ti in t:
                function = function.taylor(ti,0,order)
            new_ray = Ray(ray.base,ray.direction,function)
            self.__transform(new_ray)
            result.add(new_ray)
        return result

    def __scatter_step_one_vertex(self, order, order_max, case):
        """Calculates scattering at a particular order for diagrams with only one vertex"""
        result = deepcopy(self)
        print("order",order)
        new_rays = self.rays
        local_diagram = deepcopy(result)
        local_diagram.rays = [local_ray for ray in local_diagram.rays for local_ray in self.__local_rays(ray,(0,0))]
        result.rays += local_diagram.__scatter(order,order_max,case).rays
        return result
        
    def __scatter(self, order, order_max, case):
        """The heart of the code. Computes the scattering for a local diagram (i.e. all base points are equal)"""
        result = Diagram()
        # calculate thetaD up to order k
        X,Y = x,y
        for ray in sorted(self.rays,key=lambda ray: arctan2(-ray.direction[1],-ray.direction[0])):
            X = X.substitute(x=x*ray.function^(-ray.direction[1]),y=y*ray.function^ray.direction[0])
            Y = Y.substitute(x=x*ray.function^(-ray.direction[1]),y=y*ray.function^ray.direction[0])
            for ti in t:
                X = X.taylor(ti,0,order)
                Y = Y.taylor(ti,0,order)
            # for del Pezzos directly exclude most of the terms to speed up
            if case != "":
                X = taylor(X,x,0,self.__order(case)*order+1)
                Y = taylor(Y,x,0,self.__order(case)*order)
        X,Y = X/x, Y/y
        # add ray for each term
        exp_t = sorted(set(R(X).exponents() + R(Y).exponents()))
        for i in range(1,len(exp_t)):
            t_factor = prod([t[j]^exp_t[i][j] for j in range(len(t))])
            coeff_x = SR(X).coefficient(t_factor)
            coeff_y = SR(Y).coefficient(t_factor)
            for ti in t:
                coeff_x = coeff_x.substitute(ti==0)
                coeff_y = coeff_y.substitute(ti==0)
            exp_xy = sorted(set(Q(coeff_x).exponents() + Q(coeff_y).exponents()))
            for j in range(len(exp_xy)):
                exp_x = exp_xy[j][0]
                exp_y = exp_xy[j][1]
                gcd_xy = gcd(exp_x,exp_y) if gcd(exp_x,exp_y) !=0 else 1
                if exp_x == 0 and exp_y == 0:
                    coeff2_x = coeff_x
                    coeff2_y = coeff_y
                else:
                    coeff2_x = SR(coeff_x).coefficient(x^exp_x*y^exp_y)
                    coeff2_y = SR(coeff_y).coefficient(x^exp_x*y^exp_y)
                coeff = gcd(coeff2_x,coeff2_y)
                sign = sgn(coeff2_x)*sgn(exp_y) if coeff2_x != 0 else -sgn(coeff2_y)*sgn(exp_x)
                function = 1+sign*coeff*t_factor*x^exp_x*y^exp_y
                base = next(iter(self.rays)).base
                ray = Ray(base,(exp_x/gcd_xy,exp_y/gcd_xy),function)
                if function == 1:
                    continue
                if case != "":
                    ord = exp_x/self.__order(case)
                    if ord > order:
                        continue
                    ray = Ray(next(iter(self.rays)).base,(exp_x/gcd_xy,exp_y/gcd_xy),function,order=ord)
                self.__transform(ray)
                result.add(ray)
        return result

    def __scatter_prim(self, ray1, ray2):
        for int_pt in ray1.intersection(ray2):
            if isinstance(int_pt,Ray):
                continue
            f1 = R(ray1.function)
            f2 = R(ray2.function)
            if len(f1.exponents()) == 2 and len(f2.exponents()) == 2 and max(np.add(f1.exponents()[1],f2.exponents()[1])) == 1:
                term = (ray1.function-1)*(ray2.function-1)*abs(np.linalg.det([ray1.direction,ray2.direction]))
                exp_t = R(term).exponents()
                if len(exp_t) == 1 and max(exp_t[0]) == 1:
                    t_factor = prod([t[j]^exp_t[0][j] for j in range(len(t))])
                    coeff = SR(term).coefficient(t_factor)
                    exp_xy = Q(coeff).exponents()
                    if len(exp_xy) == 1:
                        exp_x = exp_xy[0][0]
                        exp_y = exp_xy[0][1]
                        w_out = gcd(exp_x,exp_y)
                        base = int_pt
                        direction = (exp_x/w_out,exp_y/w_out)
                        function = 1 + w_out*term
                        return Ray(base,direction,function)

    ### PRIVATE METHODS - BROKEN LINES ###

    def __broken_iter(self, ray_incoming, ray_list, result, order, case, exponents_in):
        x_order = order * self.__order(case) if case != '' else order
        exp_t_in = R(ray_incoming.transforms[-1].function).exponents()[0] # TODO shorter
        t_factor_in = prod([t[j]^exp_t_in[j] for j in range(len(t))])
        coeff_in = SR(ray_incoming.transforms[-1].function).coefficient(t_factor_in)
        incoming_orders = [exp_xy[0] for exp_xy in Q(coeff_in).exponents()]
        incoming_order = max(incoming_orders) if len(incoming_orders) > 0 else 0
        if exponents_in == 'all' or ray_incoming.transforms[-1].exponent in exponents_in:
            result.append(ray_list)
        int_dict = self.__int_dict(ray_incoming)
        int_pts = sorted(int_dict.keys(),key=lambda pt: np.linalg.norm(np.subtract(pt,ray_incoming.base)))
        for int_pt in int_pts:
            for ray in int_dict[int_pt]:
                ray_loc = ray.local_rep(int_pt)
                ray_in = ray_incoming.local_rep(int_pt)
                if all(np.isclose(int_pt,ray_in.base)): # TODO
                    continue
                # check if breaking would be of too high order
                ray_temp = Ray(int_pt,np.add(ray_loc.exponent,ray_in.exponent),1)
                min_order = ray_temp.direction[0] if case == '' else ray_temp.direction[0] + ray_temp.direction[1] * self.__speed_up(ray_temp,case,order)
                if min_order > x_order: # TODO weighted direction
                    continue
                exp_t_ray_in = R(ray_in.function).exponents()[0] # TODO shorter
                t_factor_ray_in = prod([t[j]^exp_t_ray_in[j] for j in range(len(t))])
                coeff_ray_in = SR(ray_in.function).coefficient(t_factor_ray_in)
                p_exp = Q(ray_in.function).exponents()[0] if sum(exp_t_ray_in) == 0 else Q(coeff_ray_in).exponents()[0]
                term = ray_loc.function^abs(ray_loc.direction[1]*p_exp[0]-ray_loc.direction[0]*p_exp[1])*ray_in.function
                for ti in t:
                    term = term.taylor(ti,0,order)
                if case != '':
                    term = taylor(term,x,0,x_order)
                # add ray for each term
                exp_t = R(term).exponents()
                for i in range(1,len(exp_t)):
                    t_factor = prod([t[j]^exp_t[i][j] for j in range(len(t))])
                    coeff = SR(term).coefficient(t_factor)
                    exp_xy = Q(coeff).exponents()
                    for j in range(len(exp_xy)):
                        exp_x = exp_xy[j][0]
                        exp_y = exp_xy[j][1]
                        gcd_xy = gcd(exp_x,exp_y)
                        coeff2 = SR(coeff).coefficient(x^exp_x*y^exp_y)
                        incoming_new = Ray(int_pt,(exp_x/gcd_xy,exp_y/gcd_xy),coeff2*t_factor*x^exp_x*y^exp_y)
                        self.__transform(incoming_new)
                        exp_t_last = R(incoming_new.transforms[-1].function).exponents()[0] # TODO shorter
                        t_factor_last = prod([t[j]^exp_t_last[j] for j in range(len(t))])
                        coeff_last = SR(incoming_new.transforms[-1].function).coefficient(t_factor_last)
                        if case == '':
                            new_orders = [exp_xy[0] for exp_xy in Q(coeff_last).exponents()]
                        else:
                            new_orders = [exp_xy[0]+self.__speed_up(incoming_new,case,order)*abs(exp_xy[1]) for exp_xy in Q(coeff_last).exponents()]
                        new_order = max(new_orders) if len(new_orders) > 0 else 0
                        if new_order <= x_order: # add
                            list_new = [deepcopy(ray) for ray in ray_list]
                            # set endpoint and cut transforms-list
                            ind = list_new[-1].transforms.index(list_new[-1].local_rep(int_pt))
                            list_new[-1].transforms[ind].endpoint = int_pt
                            list_new[-1].transforms = list_new[-1].transforms[:ind+1]
                            # add new line part and do next recursive step
                            list_new.append(incoming_new)
                            self.__broken_iter(incoming_new,list_new,result,order,case,exponents_in)
                            sys.stdout.write("-")
                            sys.stdout.flush()
    
    ### PRIVATE METHODS - WALLCROSSING ###
    
    def __cross(self, W, f, n, order):
        Wtemp = 0
        Cs = self.__coeffs(W)    
        for C in Cs:
            a,b,c = C
            p = n[0]*a + n[1]*b
            Wtemp += c*x^a*y^b*pow(f,p)
        TempCs = self.__coeffs(Wtemp)
        Wnew = 0
        for C in TempCs:
            a,b,c = C
            if abs(a)<3*(order+1):
                Wnew += c*x^a*y^b
        return expand(Wnew)

    def __coeffs(self, W):
    # return list with entries [a,b,c] for each summand c*x^a*y^b of W
        Cs = []
        A = W.coefficients(x)
        for p in A:
            a = p[1]
            B = p[0].coefficients(y)
            for q in B:
                b = q[1]
                Cs.append([a,b,q[0]])
        return Cs

    ### PRIVATE METHODS - CLASSES

    def __calc_classes(self, order):
        self.class_rays = self.__class_rays(self.case,order)
        for ray in self.rays:
            if ray.tropical_curve == []:
                self.__calc_trop(ray)
            if ray.direction == (1,0): # TODO
                self.__calc_class(ray)

    def __calc_class_broken(self, broken_line):
        trop_curve = []
        for ind, part in enumerate(broken_line):
            if ind > 0: # to get right multiplicity of tropical curve # TODO
                part.function = part.function / broken_line[ind-1].function
            self.__calc_trop(part)
            trop_curve += part.tropical_curve
            part.tropical_curve = []
            if ind > 0:
                part.function = part.function * broken_line[ind-1].function
        #save tropical curve and class of the broken line in its last line segment
        broken_line[-1].tropical_curve = trop_curve
        self.__calc_class(broken_line[-1])

    def __calc_class(self, ray):
        ray.curve_class = tuple([0 for i in self.class_rays[0].function])
        for class_ray in self.class_rays:
            # calculate intersection points
            int_pts = []
            for edge in ray.tropical_curve:
                intersections = class_ray.intersection(edge)  
                if isinstance(intersections,Ray) or len(intersections) == 0:
                    continue
                for int_pt in intersections:
                    if isinstance(int_pt,Ray):
                        continue
                    is_new = True
                    for pt in int_pts:
                        if all(np.isclose(int_pt,pt)):
                            is_new = False
                    if is_new:
                        int_pts.append(int_pt)
            for int_pt in int_pts:
                adjacent = []
                for edge in ray.tropical_curve:
                    if edge.contains(int_pt):
                        adjacent.append(edge)
                weights = [edge.function for edge in adjacent]
                mult = max(weights)
                ray.curve_class = tuple(np.add(ray.curve_class,np.multiply(mult,class_ray.function)))

    def __calc_trop(self, ray):
        outgoing_edge = deepcopy(ray)
        outgoing_edge.function = 1 # has to be 1 to get right mult
        self.__transform(outgoing_edge)
        ray.tropical_curve = [outgoing_edge]
        # t_exp: R(ray.function).exponents()[0] for broken lines, R(ray.function).exponents()[1] for walls
        t_exp = R(ray.function).exponents()[-1]
        for ray2 in self.rays: # we use that self.rays is ordered such that parents come before their childs
            if ray2.contains(ray.base) and not all(np.isclose(ray.base,ray2.base)) and not ray.direction == ray2.direction:
                t_exp2 = R(ray2.function-1).exponents()[0]
                weight = 0
                while min(t_exp) > -1:
                    t_exp = np.subtract(t_exp,t_exp2)
                    weight += 1
                t_exp = np.add(t_exp,t_exp2) # TODO ugly
                weight -= 1
                if weight > 0:
                    edge = deepcopy(ray2)
                    edge.endpoint = ray.base
                    edge.function = weight
                    self.__transform(edge)
                    ray.tropical_curve.append(edge)
                    for edge2 in ray2.tropical_curve:
                        if edge2.endpoint != (): # means if not the same
                            ray.tropical_curve.append(edge2)
                if sum(t_exp) == 0:
                    break
    
    ### PRIVATE METHODS - OUTPUT ###

    """Combines rays with the same base and direction"""
    def __combine_rays(self, condition = 'same'):
        result = copy(self)
        ray_list = self.rays
        i = 0
        while i < len(ray_list):
            ray1 = ray_list[i]
            function = ray1.function
            j = i+1
            while j < len(ray_list):
                ray2 = Dn[j]
                combine = False
                if all(np.isclose(ray1.direction,ray2.direction)):
                    if condition == 'same' and all(np.isclose(ray1.base,ray2.base)) and all(np.isclose(ray1.endpoint,ray2.endpoint)):
                        combine = True
                    if condition == 'overlapping':
                        if ray1.contains(ray2.base):
                            combine = True
                        if ray2.contains(ray1.base):
                            combine = True
                            ray1,ray2 = ray2,ray1
                    function *= ray2.function
                    del ray_list[j]
                else:
                    j = j+1
            # create new ray and add to result
            ray = deepcopy(ray1)
            ray.function = function
            result.add(ray)
            i = i+1
        return result

    """Converts rays in a diagram to a printable form"""
    def __diag_conv(self, special = [], colors = False, functions = False, directions = []):
        if len(self.rays) == 0:
            return Diagram()
        special = special if isinstance(special,Diagram) else Diagram(special)
        ray_parts = [part for ray in self.rays for part in ray.transforms]
        mon_parts = list(self.monodromy_rays)
        broken_parts = [part for broken_lines in list(self.broken_lines.values()) for broken_line in broken_lines for ray in broken_line for part in ray.transforms]
        special_parts = [part for ray in special for part in ray.transforms]
        # Determine bounds of the displayed diagram
        bound_l = min(min([(ray.base[0],ray.base[0]+ray.direction[0]) for ray in ray_parts + mon_parts + broken_parts + special_parts]))
        bound_r = max(max([(ray.base[0],ray.base[0]+ray.direction[0]) for ray in ray_parts + mon_parts + broken_parts + special_parts]))
        bound_u = max(max([(ray.base[1],ray.base[1]+ray.direction[1]) for ray in ray_parts + mon_parts + broken_parts + special_parts]))
        bound_d = min(min([(ray.base[1],ray.base[1]+ray.direction[1]) for ray in ray_parts + mon_parts + broken_parts + special_parts]))
        # add rays
        result = Diagram()
        for ray in ray_parts:
            length = self.__length(ray,bound_l,bound_r,bound_u,bound_d)
            print_fct = ray.direction in directions
            new_ray = self.__ray_conv(ray,length,print_fct)
            new_ray.special = False
            result.add(new_ray)
        for ray in special_parts:
            length = self.__length(ray,bound_l,bound_r,bound_u,bound_d)
            print_fct = ray.direction in directions
            new_ray = self.__ray_conv(ray,length,print_fct)
            new_ray.special = True
            result.add(new_ray)
        return result

    def __length(self, ray, bound_l, bound_r, bound_u, bound_d):
        # Determine length factor of displayed lines.
        lengths = []
        if ray.direction[0] != 0:
            lengths.append((bound_l-ray.base[0])/ray.direction[0] if ray.direction[0] < 0 else (bound_r-ray.base[0])/ray.direction[0])
        if ray.direction[1] != 0:
            lengths.append((bound_d-ray.base[1])/ray.direction[1] if ray.direction[1] < 0 else (bound_u-ray.base[1])/ray.direction[1])
        return min(lengths) if len(lengths) > 0 else 0

    def __ray_conv(self, ray, length, print_fct):
        # Determine function to be displayed
        function = ""
        if print_fct and len(ray.function) > 0:
            term = str(R(ray.function-1)).replace("(","").replace(")","").replace("*","")
            term = re.sub('t_[0-9]+(\\^[0-9]+)*','',term)
            term = re.sub('([\\^_])([0-9][0-9]+)',r'\1{\2}',term)
            term = re.sub('^([^-])(.*)',r'+\1\2',term)
            function = function + "(1" + term + ")"
            function = str(function)
        result = Ray(ray.base,ray.direction,function,order=ray.order)
        # Determine main direction for node positioning.
        if abs(ray.direction[0]) > abs(ray.direction[1]):
            node_position = "right" if ray.direction[0] > 0 else "left"
        else:
            node_position = "above" if ray.direction[1] > 0 else "below"
        result.node_position = node_position
        result.arrow = ray.endpoint == ()
        # Determine endpoint
        endpoint = tuple(np.add(ray.base,np.multiply(length,ray.direction)))
        if len(ray.endpoint) > 0 and np.linalg.norm(np.subtract(ray.base,ray.endpoint)) < np.linalg.norm(np.subtract(ray.base,endpoint)):
            result.endpoint = ray.endpoint
        else:
            result.endpoint = endpoint
        return result

    def __class_string(self, ray):
        class_names = self.__class_names(self.case)
        if ray.curve_class == ():
            return ""
        if self.case == "":
            return " curve class: " + str(ray.curve_class)
        class_str = ""
        if ray.curve_class[0] == 1:
            class_str += class_names[0]
        elif not ray.curve_class[0] == 0:
            class_str += str(ray.curve_class[0]) + class_names[0]
        for i in range(1,len(ray.curve_class)):
            if not ray.curve_class[i-1] == 0 and ray.curve_class[i] > 0:
                class_str += "+"
            if ray.curve_class[i] == 1:
                class_str += class_names[i]
            elif ray.curve_class[i] == -1:
                class_str += "-" + class_names[i]
            elif not ray.curve_class[i] == 0:
                class_str += str(ray.curve_class[i]) + class_names[i]
        return " curve class: " + class_str

    ### PRIVATE METHODS - UTILITY

    def __intersections(self):
        result = []
        for ray1,ray2 in set(itertools.combinations(self.rays,2)):
            intersections = ray1.intersection(ray2)
            if isinstance(intersections,Ray) or len(intersections) == 0:
                continue
            for int_pt in intersections:
                result.append(int_pt)
        return result

    def __int_dict(self, ray):
        result = {}
        for r in self.rays:
            intersections = r.intersection(ray)
            if isinstance(intersections,Ray) or len(intersections) == 0:
                continue
            for int_pt in intersections:
                is_new = true
                for key in result.keys():
                    if all(np.isclose(key,int_pt)):
                        is_new = False
                        result[key].append(r)
                if is_new:
                    result[int_pt] = [r]
        return result

    def __local_diagrams(self, order, case, new_rays):
        """Returns a map from intersection points to a list of corresponding rays"""
        int_dict = {}
        for ray1,ray2 in set(itertools.combinations(self.rays,2)):
            if new_rays != 'all' and ray1 not in new_rays and ray2 not in new_rays:
                continue
            intersections = ray1.intersection(ray2)
            if isinstance(intersections,Ray) or len(intersections) == 0:
                continue
            for int_pt in intersections:
                # for symmetric cases only look at singularities between 0 and 0.5, then use symmetry
                if case in ["(9)","9","P2","(8')","8'","P1xP1","(3)","3","(3a)","3a","Cubic","cubic"]:
                    if (int_pt[1] < 0 and not np.isclose(int_pt[1],0)) or (int_pt[1] > 0.5 and not np.isclose(int_pt[1],0.5)) or (int_pt[0] < 0 and not np.isclose(int_pt[0],0)):
                            continue
                # check if new
                is_new = True
                for key in int_dict.keys():
                    if all(np.isclose(key,int_pt)):
                        is_new = False
                        int_dict[key] += self.__local_rays(ray1,key) + self.__local_rays(ray2,key)
                if is_new:
                    int_dict[int_pt] = self.__local_rays(ray1,int_pt) + self.__local_rays(ray2,int_pt)
        result = {}
        for key in int_dict.keys():
            # throw away multiples
            int_dict[key] = set(int_dict[key])
            # if case is given, check if scattering at key produces only high order terms
            if case != '':
                high_order = True
                for ray1,ray2 in set(itertools.combinations(int_dict[key],2)):
                    direction_sum = np.add(ray1.direction,ray2.direction)
                    ord = direction_sum[0]
                    if ord <= self.__order(case) * order:
                        high_order = False
                if high_order:
                    continue
            result[key] = Diagram(int_dict[key]) # TODO return iterator
        return result
    
    def __local_rays(self, ray, pt):
        # if pt is not the base, split ray
        r = ray.local_rep(pt)
        result = []
        if all(np.isclose(ray.base,pt)):
            result.append(r)
        else:
            result.append(Ray(pt,r.direction,r.function))
            result.append(Ray(pt,(-r.direction[0],-r.direction[1]),r.function))
        return result

    def __transform(self, ray):
        ind = 0
        while self.max_trafos == 'infty' or ind < self.max_trafos:
            ind = ind + 1
            last = ray.transforms[-1]
            int_dict = dict([(mon,min([int_pt for int_pt in last.intersection(mon) if not all(np.isclose(int_pt,last.base))], key = lambda int_pt: np.linalg.norm(np.subtract(int_pt,last.base)))) for mon in self.monodromy_rays if len([int_pt for int_pt in last.intersection(mon) if not all(np.isclose(int_pt,last.base))]) > 0])
            if len(int_dict) == 0:
                break
            mon = min(int_dict.keys(), key = lambda mon: np.linalg.norm(np.subtract(int_dict[mon],last.base)))
            trafo = mon.transform(last)
            last.endpoint = int_dict[mon]
            ray.transforms.append(trafo)
    
    ### PRIVATE METHODS - DEL PEZZOS ###
    
    def __kinks(self, case):
        if case in ["P2","(9)","9"]: return [3]
        if case in ["P1xP1","(8')","8'"]: return [2,2]
        if case in ["(8'a)","8'a"]: return [2,2,2]
        if case in ["dP1","F1","H1","(8)","8"]: return [1,2,3,2]
        if case in ["dP2","(7)","7"]: return [1,2,2,1,1]
        if case in ["6c"]: return [1,2,3]
        if case in ["Cubic","cubic","(3a)","3a","(3)","3"]: return [1]
    
    def __order(self, case):
        #return max(self.__kinks(case))
        return max(sum(self.__kinks(case)[:ceil(len(self.__kinks(case))/2)]),sum(self.__kinks(case)[floor(len(self.__kinks(case))/2):]))
    
    def __classes(self, case):
        if case in ["P2","(9)","9"]: return [(1,)]
        if case in ["P1xP1","(8')","8'"]: return [(1,0),(0,1)]
        if case in ["(8'a)","8'a"]: return []
        if case in ["dP1","F1","H1","(8)","8"]: return [(0,1),(1,-1),(1,0),(1,-1)]
        if case in ["dP2","(7)","7"]: return [(0,1,0),(1,-1,0),(1,0,-1),(0,0,1),(1,-1,-1)]
        if case in ["6c"]: return [(0,1),(1,-1),(1,0)]
        if case in ["Cubic","cubic","(3a)","3a","(3)","3"]: return []
      
    def __lengths(self, case):
        if case in ["P2","(9)","9"]: return [1]
        if case in ["P1xP1","(8')","8'"]: return [1]
        if case in ["(8'a)","8'a"]: return [1,1,2]
        if case in ["dP1","F1","H1","(8)","8"]: return [1,1,1,1]
        if case in ["dP2","(7)","7"]: return [1,1,1,1,1]
        if case in ["6c"]: return [3,1,2]
        if case in ["Cubic","cubic","(3a)","3a","(3)","3"]: return [3]
    
    def __class_names(self, case):
        if case in ["P2","(9)","9"]: return ["L"]
        if case in ["P1xP1","(8')","8'"]: return ["L1","L2"]
        if case in ["(8'a)","8'a"]: return []
        if case in ["dP1","F1","H1","(8)","8"]: return ["L","E"]
        if case in ["dP2","(7)","7"]: return ["L","E1","E2"]
        if case in ["6c"]: return ["L","E"]
        if case in ["Cubic","cubic","(3a)","3a","(3)","3"]: return []
        
    def __class_rays(self, case, order):
        eps = (-order,-0.01) #TODO
        result = []
        kinks = self.__kinks(case)
        classes = self.__classes(case)
        lengths = self.__lengths(case)
        base_top = (0,1-lengths[0])
        base_bottom = (0,0+lengths[-1])
        slope_top = 0
        slope_bottom = 0
        for i in range(order):
            for j in range(len(kinks)):
                slope_top += kinks[j]
                slope_bottom += kinks[-1-j]
                base_top = np.add(base_top,(-slope_top,lengths[j]))
                base_bottom = np.add(base_bottom,(-slope_bottom,-lengths[-1-j]))
                result.append(Ray(np.add(base_top,eps),(1,0),classes[j]))
                result.append(Ray(np.add(base_bottom,eps),(1,0),classes[-1-j]))
        return result

    def __get_rays(self, case, order):
        kinks = self.__kinks(case)
        lengths = self.__lengths(case)
        ray_list = []
        base_top = (0,lengths[0]/2)
        base_bottom = (0,lengths[0]/2)
        dir_top = (0,-1)
        dir_bottom = (0,1)
        ind = 0
        while True:
            j = ind%len(kinks)
            if ind >= ceil(len(kinks)/2) and dir_top[0]+kinks[j] > self.__order(case)*order and dir_bottom[0]+kinks[-1-j] > self.__order(case)*order:
                break
            # rays
            ray_list.append(Ray(base_top, (-dir_top[0],-dir_top[1]), (1+t[4*ind]*x^(-dir_top[0])*y^(-dir_top[1]))^lengths[j]))
            ray_list.append(Ray(base_bottom, (-dir_bottom[0],-dir_bottom[1]), (1+t[4*ind+1]*x^(-dir_bottom[0])*y^(-dir_bottom[1]))^lengths[j]))
            # calculate new bases and directions
            base_top = np.add(base_top,(-dir_top[0]-kinks[j]/2,lengths[j]/2+lengths[(j+1)%len(lengths)]/2))
            base_bottom = np.add(base_bottom,(-dir_bottom[0]-kinks[-1-j]/2,-lengths[-j]/2-lengths[(-j-1)%len(lengths)]/2))
            dir_top = (dir_top[0]+kinks[j],-1)
            dir_bottom = (dir_bottom[0]+kinks[-1-j],1)
            # rays
            ray_list.append(Ray(base_top,dir_top,(1+t[4*ind+3]*x^dir_top[0]*y^dir_top[1])^lengths[j]))
            ray_list.append(Ray(base_bottom,dir_bottom,(1+t[4*ind+2]*x^dir_bottom[0]*y^dir_bottom[1])^lengths[j]))
            # increment index
            ind += 1
        return ray_list

    def __get_monodromy_rays(self, case, order):
        kinks = self.__kinks(case)
        lengths = self.__lengths(case)
        monodromy_ray_list = []
        base_top = (0,lengths[0]/2)
        base_bottom = (0,lengths[0]/2)
        dir_top = (0,-1)
        dir_bottom = (0,1)
        ind = 0
        while True:
            j = ind%len(kinks)
            if ind >= ceil(len(kinks)/2) and dir_top[0]+kinks[j] > self.__order(case)*order and dir_bottom[0]+kinks[-1-j] > self.__order(case)*order:
                break
            # monodromy rays
            dir_mon_top = (-2-dir_top[0]*lengths[j],lengths[j]) if lengths[j] % 2 == 1 else (-1-dir_top[0]*lengths[j]/2,lengths[j]/2)
            dir_mon_bottom = (-2-dir_bottom[0]*lengths[-j],-lengths[-j]) if lengths[-j] % 2 == 1 else (-1-dir_bottom[0]*lengths[-j]/2,-lengths[-j]/2)
            matrix_top = [[1-lengths[j]*dir_top[0],-lengths[j]*dir_top[0]^2],[lengths[j],1+lengths[j]*dir_top[0]]]
            matrix_bottom = [[1+lengths[-j]*dir_top[0],lengths[-j]*dir_top[0]^2],[-lengths[-j],1-lengths[-j]*dir_top[0]]]
            end_top = np.add(base_top,np.multiply(1/2,dir_mon_top))
            end_bottom = np.add(base_bottom,np.multiply(1/2,dir_mon_bottom))
            monodromy_ray_list.append(Ray(base_top, dir_mon_top, matrix_top, endpoint = end_top))
            monodromy_ray_list.append(Ray(base_bottom, dir_mon_bottom, matrix_bottom, endpoint = end_bottom))
            # calculate new bases and directions
            base_top = np.add(base_top,(-dir_top[0]-kinks[j]/2,lengths[j]/2+lengths[(j+1)%len(lengths)]/2))
            base_bottom = np.add(base_bottom,(-dir_bottom[0]-kinks[-1-j]/2,-lengths[-j]/2-lengths[(-j-1)%len(lengths)]/2))
            dir_top = (dir_top[0]+kinks[j],-1)
            dir_bottom = (dir_bottom[0]+kinks[-1-j],1)
            # monodromy rays
            j_1 = (j+1)%len(lengths)
            dir_mon_top = (-2+dir_top[0]*lengths[j_1],-lengths[j_1]) if lengths[j_1] % 2 == 1 else (-1+dir_top[0]*lengths[j_1]/2,-lengths[j_1]/2)
            dir_mon_bottom = (-2+dir_bottom[0]*lengths[-j-1],lengths[-j-1]) if lengths[-j-1] % 2 == 1 else (-1+dir_bottom[0]*lengths[-j-1]/2,lengths[-j-1]/2)
            matrix_top = [[1+lengths[j_1]*dir_top[0],lengths[j_1]*dir_top[0]^2],[-lengths[j_1],1-lengths[j_1]*dir_top[0]]]
            matrix_bottom = [[1+lengths[-j-1]*dir_top[0],-lengths[-j-1]*dir_top[0]^2],[lengths[-j-1],1-lengths[-j-1]*dir_top[0]]]
            end_top = np.add(base_top,np.multiply(1/2,dir_mon_top))
            end_bottom = np.add(base_bottom,np.multiply(1/2,dir_mon_bottom))
            monodromy_ray_list.append(Ray(base_top, dir_mon_top, matrix_top, endpoint = end_top))
            monodromy_ray_list.append(Ray(base_bottom, dir_mon_bottom, matrix_bottom, endpoint = end_bottom))
            # increment index
            ind += 1
        return monodromy_ray_list
    
    def __get_diagrams(self, case, order, path):
        result = [self]
        if not os.path.isfile(path):
            return result
        text_file = open(path,'r')
        lines = text_file.readlines()
        ind_first = 0
        ind_last = len(lines) - 1
        # find starting index
        for i in range(len(lines)):
            if lines[i] == case + "\n":
                ind_first = i+1
            if lines[i] == "\n" and ind_first > 0:
                ind_last = i-1
        if ind_first == 0:
            #print("Did not find data.")
            text_file.close()
            return result
        if not (ind_last - ind_first + 1) % 5 == 0:
            print("Wrong format.")
            text_file.close()
            return result
        for i in range((ind_last - ind_first + 1) / 5):
            if i+1 > order:
                break
            if not lines[ind_first+5*i] == "order " + str(i+1) + "\n" or not lines[ind_first+5*i+1] == "rays\n" or not lines[ind_first+5*i+3] == "broken lines\n":
                print("Wrong format.")
                text_file.close()
                return result
                break
            diag_string = lines[ind_first+5*i+2]
            if diag_string == "[]\n":
                break
            ray_list = [(ray[0],ray[1],ray[2]) for ray in eval(lines[ind_first+5*i+2].replace('^','**'))]
            diagram = Diagram(ray_list,self.__get_monodromy_rays(case,i+1))
            # add broken lines
            for br_lines in eval(lines[ind_first+5*i+4].replace('^','**')):
                broken_lines = []
                for br_line in br_lines:
                    broken_line = []
                    for ray in br_line:
                        broken_line.append(Ray(np.array(ray[0], dtype=float),np.array(ray[1], dtype=float),ray[2]))
                    broken_lines = append(broken_line)
                diagram.broken_lines.append(broken_lines) # TODO dict
            result.append(diagram)
        text_file.close()
        return result

    # y factor for cases
    def __speed_up(self, ray, case, order):
        if ray.direction[1] == 0:
            return 0
        kinks = self.__kinks(case)
        lengths = self.__lengths(case)
        base_top = (0,lengths[0]/2)
        base_bottom = (0,lengths[0]/2)
        dir_top = (0,-1)
        dir_bottom = (0,1)
        for ind in range(order):
            ray_top = Ray(base_top,(-dir_top[0],-dir_top[1]),1)
            ray_bottom = Ray(base_bottom,(-dir_bottom[0],-dir_bottom[1]),1)
            if ray.direction[1] > 0 and len(ray_bottom.intersection(ray)) > 0:
                return ray_bottom.direction[0]
            if ray.direction[1] < 0 and len(ray_top.intersection(ray)) > 0:
                return ray_top.direction[0]
            j = ind%len(kinks)
            base_top = (base_top[0]-dir_top[0]-kinks[j]/2,base_top[1]+lengths[j]/2+lengths[(j+1)%len(lengths)]/2)
            base_bottom = (base_bottom[0]-dir_bottom[0]-kinks[-1-j]/2,base_bottom[1]-lengths[-j]/2-lengths[(-j-1)%len(lengths)]/2)
            dir_top = (dir_top[0]+kinks[j],-1)
            dir_bottom = (dir_bottom[0]+kinks[-1-j],1)
            ray_top = Ray(base_top,dir_top,1)
            ray_bottom = Ray(base_bottom,dir_bottom,1)
            if ray.direction[1] > 0 and len(ray_top.intersection(ray)) > 0:
                return ray_top.direction[0]
            if ray.direction[1] < 0 and len(ray_bottom.intersection(ray)) > 0:
                return ray_bottom.direction[0]
        return order*self.__order(case)

    # used in __translate
    def __t_shift(self, i, j): # TODO slow
        if i >= 0:
            for k in range(i+1):
                j = self.__shift_left(j)
        if i < 0:
            for k in range(abs(i)):
                j = self.__shift_right(j)
        return j
    
    def __shift_left(self, j):
        if j == 1:
            return 3
        if j == 2:
            return 0
        if j%4 == 0:
            return j+4
        if j%4 == 1:
            return j-4
        if j%4 == 2:
            return j-4
        if j%4 == 3:
            return j+4
    
    def __shift_right(self, j):
        if j == 0:
            return 2
        if j == 3:
            return 1
        if j%4 == 0:
            return j-4
        if j%4 == 1:
            return j+4
        if j%4 == 2:
            return j+4
        if j%4 == 3:
            return j-4

    # for (P2,E) do calculations in a given domain, then add other rays by affine transformation via this method
    def __translate(self, k, case):
        result = []
        # reflect
        D2 = self.rays
        for ray in self.rays:
            if (0 < ray.base[1] and ray.base[1] < 0.5 and not np.isclose(ray.base[1],0.5)) or np.isclose(ray.base[1],0):
                f = ray.function
                f = f.substitute(y == y^-1)
                f = f.substitute([t[i] == t[i + (-1)^i] for i in range(len(t))])
                ray2 = Ray((ray.base[0],1-ray.base[1]),(ray.direction[0],-ray.direction[1]),f,order=ray.order)
                D2.append(ray2)
        # affine trafo
        for ray in D2:
            result.append(ray)
            i_min = -k+1 #-ray.order
            i_max = k-1 #ray.order
            for i in range(i_min,i_max):
                if i >= 0 and np.isclose(ray.base[1],0):
                    continue
                if i < 0 and np.isclose(ray.base[1],1):
                    continue
                e = 1 if i >= 0 else 0
                base = (ray.base[0]-self.__order(case)*(i+e)*ray.base[1]-self.__order(case)/2*i^2-self.__order(case)/2*abs(i),ray.base[1]+i+e)
                direction = (ray.direction[0]-self.__order(case)*(i+e)*ray.direction[1],ray.direction[1])
                f = ray.function
                f = f.substitute([t[j] == t[self.__t_shift(i,j)%len(t)] for j in range(len(t))])
                if ray.direction[0] != 0:
                    f = f.substitute(x==x^(direction[0]/ray.direction[0]))
                if ray.direction[1] != 0:
                    f = f.substitute(y==y^(direction[1]/ray.direction[1]))
                ray2 = Ray(base,direction,f,order=ray.order+abs(i))
                result.append(ray2)
        return Diagram(set(result))

class Ray:
    """Class of rays"""
    def __init__ (self, base:tuple = (), direction:tuple = (), function = 1, endpoint:tuple = (), order:int = 0):
        self.base = tuple(base)
        self.direction = tuple(direction)
        self.function = function
        self.endpoint = tuple(endpoint)
        self.order = order
        self.transforms = [self]
        self.tropical_curve = []
        self.curve_class = ()
        self.exponents = []
        self.exponent= ()
        self.t_order = 0
        if not isinstance(function,(list,tuple,str)):
            self.t_exps = R(function).exponents()
            t_exps = [max(t_exp) for t_exp in self.t_exps if not max(t_exp) == 0]
            if len(t_exps) > 0:
                self.t_order = min(t_exps)
                for ti in t:
                    function = function.substitute(ti==1)
            self.exponents = Q(function).exponents()
            self.exponent = tuple(self.exponents[-1])
        
    def __hash__ (self):
        return hash(str(self.base)+str(self.direction)+str(self.function))
    
    def __repr__ (self):
        output = "Ray: " + ' -- '.join([ray.__str_part() for ray in self.transforms])
        return output
        
    def __getitem__ (self, key):
        if key == 0:
            return self.base
        if key == 1:
            return self.direction
        if key == 2:
            return self.function
    
    def __len__ (self):
        return 2
        
    def __eq__ (self, other):
        return all(np.isclose(self.base,other.base)) and all(np.isclose(self.direction,other.direction)) and self.function == other.function
    
    def local_rep(self, pt):
        """Gives the local representative of the ray at the point pt."""
        for ray in self.transforms:
            if ray.__contains_point(pt):
                return ray
        return self #TODO
    
    def set_endpoint(self, pt):
        self.local_rep(pt).endpoint = pt
        ind = self.transforms.index(self.local_rep(pt))
        self.transforms = self.transforms[:ind+1]

    def intersection(self, other):
        result = []
        if self.contains(other):
            return other
        if other.contains(self):
            return self
        for part1 in self.transforms:
            for part2 in other.transforms:
                if part1.__intersection_part(part2) != ():
                    result.append(part1.__intersection_part(part2))
        return result

    def contains(self, other):
        """Checks whether the ray contains another ray or a point."""
        if not isinstance(other,Ray):
            for part in self.transforms:
                if part.__contains_point(other):
                    return True
            return False
        for part in self.transforms:
            if part.__contains_part(other): # If self contains other, it already contains the first part (no transform) of other
                return True
        return False

    def transform(self, other):
        intersections = [int_pt for int_pt in self.intersection(other) if not all(np.isclose(int_pt,other.base))]
        if len(intersections) == 0:
            return other
        int_pt = min(intersections, key = lambda int_pt: np.linalg.norm(np.subtract(int_pt,other.base)))
        if len(int_pt) == 0 or all(np.isclose(int_pt,other.base)):
            return other
        tf = self.function
        base = np.add(self.base,np.matmul(tf,np.subtract(int_pt,self.base)))
        direction = np.matmul(tf,other.direction)
        direction = (int(round(direction[0])),int(round(direction[1])))
        function = other.function.substitute([x==x^int(round(tf[0][0]))*y^int(round(tf[1][0])),y==x^int(round(tf[0][1]))*y^int(round(tf[1][1]))])
        return Ray(base,direction,function)
    
    def __intersection_part(self, other):
        """Returns the intersection of the first part of the rays ray1 and ray2: [], int(2) or one of the rays"""
        if self.__contains_part(other):
            return other
        if other.__contains_part(self):
            return self
        den = np.linalg.det([self.direction,other.direction])
        if den != 0:
            length1 = np.linalg.det([other.direction,np.subtract(self.base,other.base)])/den
            length2 = np.linalg.det([self.direction,np.subtract(self.base,other.base)])/den
            if (length1 > 0 or np.isclose(length1,0)) and (length2 > 0 or np.isclose(length2,0)):
                int_pt = (self.base[0]+length1*self.direction[0],self.base[1]+length1*self.direction[1])
                if len(self.endpoint) > 0 and np.linalg.norm(np.subtract(self.base,self.endpoint)) < np.linalg.norm(np.subtract(self.base,int_pt)):
                    return ()
                if len(other.endpoint) > 0 and np.linalg.norm(np.subtract(other.base,other.endpoint)) < np.linalg.norm(np.subtract(other.base,int_pt)):
                    return ()
                return int_pt
        return ()

    def __contains_part(self, other):
        """Checks whether the ray contains the other ray"""
        if all(np.isclose(self.direction,(0,0))):
            return all(np.isclose(self.base,other.base))
        if all(np.isclose(other.direction,(0,0))):
            return self.__contains_point(other.base)
        if not np.isclose(np.linalg.det([self.direction,other.direction]),0):
            return false
        return np.dot(self.direction,other.direction) > 0 and self.__contains_point(other.base)
    
    def __contains_point(self, pt):
        if len(self.endpoint) > 0 and np.linalg.norm(np.subtract(self.base,pt)) > np.linalg.norm(np.subtract(self.base,self.endpoint)): # new
            return false
        if np.isclose(self.direction[0],0):
            term = (pt[1]-self.base[1])/self.direction[1]
            return np.isclose(self.base[0],pt[0]) and (term > 0 or np.isclose(term,0))
        if np.isclose(self.direction[1],0):
            term = (pt[0]-self.base[0])/self.direction[0]
            return np.isclose(self.base[1],pt[1]) and (term > 0 or np.isclose(term,0))
        term1 = abs((pt[0]-self.base[0])/self.direction[0] - (pt[1]-self.base[1])/self.direction[1])
        term2 = (pt[0]-self.base[0])/self.direction[0]
        return (term1 < 0 or np.isclose(term1,0)) and (term2 > 0 or np.isclose(term2,0))
    
    def __str_part(self):
        result = ""
        result += "base : " + str((round(self.base[0],3),round(self.base[1],3))) if len(self.base) > 0 else "base: ()"
        result += " direction: " + str((int(round(self.direction[0])),int(round(self.direction[1])))) if len(self.direction) > 0 else "direction: ()"
        if isinstance(self.function,np.ndarray) or isinstance(self.function,list):
            result += " function: " + str(np.around(self.function).astype(int).tolist())
        else:
            try:
                result += " function: " + str(R(self.function))
            except:
                result += " function: " + str(self.function)
        if len(self.endpoint) > 0:
            result += " endpoint: " + str((round(self.endpoint[0],3),round(self.endpoint[1],3)))
        return result
       
case = 'P2'
order = 2
D = Diagram(case,order)
D1 = D.scattering(order,case)
pt = (3*order,0.432)
D2 = D1.brokenlines(pt,order,case,exp_end=[(-5,0)],exp_in=[(1,0)])
print(D2)
D2.draw(colors=True,directions=[(1,0)])
D2.tex(colors=True,directions=[(1,0)])

"""
initialize(3)
D = Diagram([((-1,0),(1,0),(1+t_1*x)^3),((0,-1),(0,1),(1+t_2*y)^3)])
D1 = D.scattering(4)
print(D1)
D1.draw()
"""