# This licensed under the CCA license detailed here https://creativecommons.org/licenses/by/3.0/ from collections import Counterfrom multiprocessing import Poolimport itertools if __name__ == "__main__": output = open ('Output.txt', "r+") H,E1,E2,E3,E4,E5,d = var('H E1 E2 E3 E4 E5 d') N = 8 PROCESSES = 7 eps=1.15 """This implements all the variables we will require in our calculations. Note that the variable d is a placeholder which records the intersection with the boundary, it is the term about which we may Taylor expand. N is the maximal order of d to which we will attempt to expand. eps is a small error term to avoid collisions""" BaseRays = [vector([1,0]), vector([0,1]), vector([-1,0]), vector([-1,-1])] BaseRayMonomials = {(1,0):H-E1-E2, (0,1):H-E3-E4, (-1,0):E4, (-1,-1):H - E4 - E5} """ This is a list of lines in the base of the affine manifold together with what the corresponding toric divisor is. We use it here to calculate what direction a Ray object should end up pointing in, i.e. to define the function r""" Monomials = [H,E1,E2,E3,E4,E5] DotSigns = {H:1, E1:-1, E2:-1, E3:-1,E4:-1,E5:-1} """ A list of all the monomials we will want, together with the intersection form on them. Note that really we are conflating the monoid we are working with and the degree one part of the associated monoid ring""" def dot (v1,v2): product = (v1*v2).expand() return sum ([product.coefficient(m^2)*DotSigns[m] for m in Monomials]) """This gives a dot product method for elements of the monoid we are working over. Note that this is different from the dot product of two vectors as provided by Sage""" def project (ringelt): monoidelt = sum ([ringelt.degree(monomial)*monomial for monomial in Monomials]) return sum([ray * dot(BaseRayMonomials[(ray[0],ray[1])], monoidelt) for ray in BaseRays]) """This is the function r we talked about earlier, it provides a way to project from a curve class, explicitly telling us which components of the boundary it is meeting""" class Ray: def __init__ (self, series, direction, n = N, variable = d, preservedfactor = 1): self.series = series self.direction = direction self.normal = vector ([self.direction[1],-self.direction[0]]) self.angle = float(atan2 (direction[1],direction[0])) self.seriespowers = {} self.preservedfactor = preservedfactor def __lt__(self,other): if self.angle < other.angle: return 1 else: return 0 def __le__(self,other): if self.angle <= other.angle: return 1 else: return 0 def __eq__(self,other): if self.angle == other.angle: return 1 else: return 0 def __mul__ (self,other): if not self == other: raise Exception("These two rays are not equal") else: excess = (self.series*other.series)/(self.preservedfactor*other.preservedfactor) print (excess*(self.preservedfactor*other.preservedfactor)) return Ray (excess*(self.preservedfactor*other.preservedfactor), self.direction) #.taylor(d,0,N).expand() def __add__ (self,other): if not self == other: raise Exception("These two rays are not equal") else: return Ray ((self.series+other.series).expand()-1, self.direction) def __hash__ (self): return hash(self.series) def __repr__ (self): return str(self.direction) + " , " + str (self.series) def reset (self): self.seriespowers = {} def scatter (self, function, n): expdt = {} arguments = zip(itertools.repeat (self.normal), list((function+eps).expand().operands()[-2::-1])) sorteddts = pool.map (sortbyintersection, arguments) for dt in sorteddts: for key in dt.keys(): if key in expdt: expdt[key] += dt[key] else: expdt[key] = dt[key] maxexponent = max ([abs(key) for key in expdt.keys()]) for count in range (0, int(log(maxexponent+1, 2))+2): if count == 0: self.seriespowers[0] = 1 + d - d elif count == 1: self.seriespowers[1] = self.series else: self.seriespowers[count] = (self.seriespowers[count-1] *self.seriespowers[count-1]).taylor (d,0,n) self.seriespowers[-count] = (1/self.seriespowers[count]).taylor(d,0,n) arguments = zip(itertools.repeat(self), [exp for exp in expdt], [expdt[exp] for exp in expdt], itertools.repeat(n)) newfunction = pool.map (monomialscatter, arguments) return sum (newfunction) + d - d def store (self): return 'Ray (' + str(self.series) + ', vector ([' + str(self.direction[0])+ ',' + str(self.direction[1]) + '])), ' def sortbyintersection (args): normal, monomial = args[0], args[1] expdt = {} if not normal*project(monomial) in expdt: expdt[normal*project(monomial)] = monomial else: expdt[normal*project(monomial)] += monomial return expdtdef monomialscatter (package): ray, exp, monomial, n = package[0], package[1], package[2], package[3] if monomial.degree (d) > n: return 0 else: if exp < 0: sign = -1 else: sign = 1 base2list = ZZ(exp*sign).digits(2) counter = 1 runningtotal = 1 for exponent in base2list: if exponent == 1: runningtotal = runningtotal*ray.seriespowers[sign*counter] counter += 1 newfunction = ((monomial*runningtotal)).taylor(d,0,n) return newfunction """This function has to be separated out to ensure that it interacts well with pool. This is an annoying optimisation issue, but well worth the efficiency savings""" pool = Pool (PROCESSES) class ScatteringDiagram: def __init__(self, rays, bends): self.rays = rays self.bends = bends for ray in self.rays: ray.series = self.twist(ray.series) ray.preservedfactor = ray.series self.bends.sort() self.order() def __repr__ (self): name = "" for ray in self.rays: name += str(ray.direction) + " , " +str(self.twistinv(ray)) name += "\n" return name def twist (self, series): newseries = 0 for monomial in list(series.operands()[-1::-1]): newseries += monomial * (self.phi (project (monomial))) return newseries #This is the backend of the twisting by phi, we automatically twist the rays when they are instantiated into the scattering diagram def twistinv (self, ray): print (ray, "#") if 1 in ray.series.operands(): bracketsum = 1 for monomial in ray.series.operands()[-2::-1]: monomialdirection = project(monomial) correction = self.phi(ray.direction) exp = monomialdirection/ray.direction bracketsum += monomial / (correction^exp) newseries = bracketsum else: newseries = 1 for bracket in list(ray.series.operands()): bracketsum = 1 for monomial in bracket.operands()[-2::-1]: monomialdirection = project(monomial) correction = self.phi(ray.direction) exp = monomialdirection/ray.direction bracketsum += monomial / (correction^exp) newseries = newseries*bracketsum return newseries #This un-twists the rays when we come to print the rays of the scattering diagram. It takes each term and undoes the application of phi based on the direction the ray was heading in. def order (self): newrays = {} for ray in self.rays: if ray.angle in newrays.keys(): newrays[ray.angle] = newrays[ray.angle]*ray else: newrays[ray.angle] = ray self.rays = [item[1] for item in newrays.items()] self.rays.sort() #This orders the rays based on which angle they make clockwise measured from the negative x-axis. If there are two rays with the same angle we multiply the associated series. def scatter (self, functions, n): for ray in self.rays: print (ray.direction) functions = [ray.scatter(function,n).taylor(d,0,n) for function in functions] ray.reset() return functions #For every ray in the diagram in order we scatter the function of it def phi (self, basevector): angle = float(atan2 (basevector[1],basevector[0])) phi = 1 for ray in self.bends: if ray.angle < angle: phi = phi * (ray.series)^(ray.normal*basevector) return phi #This defines the function phi based on the angle each line make with the negative x-axis. def reset (self): for ray in self.rays: ray.reset() def store (self): string = '' for ray in self.rays: string+= ray.store() return stringif __name__ == "__main__": defects = [Ray (1+E1*d, vector ([1,0])), Ray (1+E2*d, vector ([1,0])), Ray (1+E3*d, vector ([0,1])), Ray (1+E5*d, vector ([-1,-1]))] bends = [Ray(H, vector([0,1])), Ray(H, vector([1,0])), Ray(H, vector([-1,0])), Ray(H, vector([-1,-1]))] print (defects) scatteringdiagram = ScatteringDiagram (defects, bends) """You would think that we have to twist all these by phi to produce the correct answer. This is all managed behind the scenes for us, so you just need to enter 1) the defect rays and 2) the bends of the toric model. All the defects carry a single weight of d to represent the fact that these curves all meet the boundary to order one. """ X = E1 Y = E5 vX = project(X) vY = project(Y) print (scatteringdiagram.rays) print ("<><><><>") print (scatteringdiagram) for n in range (0,N): print (n) print (scatteringdiagram.rays) output.write(str(n) + "\n") scatteredX, scatteredY = scatteringdiagram.scatter([X,Y], n) dx = ((scatteredX - X)/X).expand() dy = ((scatteredY - Y)/(Y)).expand() errorx = dx.coefficient (d,n).expand() errory = dy.coefficient (d,n).expand() print ("\n") corrections = [] for monomial in (errorx+eps).operands()[-2::-1]: basevector = project (monomial) normalisedbasevector = basevector / gcd (basevector[0], basevector[1]) normalvector = vector([normalisedbasevector[1], -normalisedbasevector[0]]) corrections.append (Ray (1+(monomial/(normalvector*vX))*d^n, -normalisedbasevector)) for monomial in (errory+eps).operands()[-2::-1]: basevector = project (monomial) normalisedbasevector = basevector / gcd (basevector[0], basevector[1]) normalvector = vector([normalisedbasevector[1], -normalisedbasevector[0]]) if vX*normalvector == 0: corrections.append (Ray (1+(monomial/(normalvector*vY))*d^n, -normalisedbasevector)) #This second if clause is needed to select only the rays which are missing, i.e. those which should be purely in the y direction. print (corrections) scatteringdiagram.rays += corrections scatteringdiagram.order() scatteringdiagram.reset() output.write (str(scatteringdiagram) + "\n")pool.close()pool.join()