# HG changeset patch
diff -r d305a40fdd8e -r 1b058b7dfc19 sage/calculus/riemann.pyx
 a AUTHORS: - Ethan Van Andel (2009): initial version - Ethan Van Andel (2009-2011): initial version and upgrades - Robert Bradshaw (2009): his "complex_plot" was adapted for plot_colored """ #***************************************************************************** #       Copyright (C) 2009 Ethan Van Andel , #       Copyright (C) 2011 Ethan Van Andel , # #  Distributed under the terms of the GNU General Public License (GPL) # include "../ext/stdsage.pxi" include "../ext/interrupt.pxi" from sage.plot.primitive import GraphicPrimitive from sage.misc.decorators import options from sage.plot.plot import list_plot, Graphics, setup_for_eval_on_grid from sage.plot.plot import list_plot, Graphics from sage.ext.fast_eval import fast_callable from sage.plot.complex_plot import ComplexPlot from sage.gsl.integration import numerical_integral import numpy as np cimport numpy as np from math import pi from math import sin from math import cos from math import sqrt from math import log # used for complex plot lightness from math import atan from cmath import exp from cmath import phase cdef double PI = pi cdef double TWOPI = 2*PI cdef I = complex(0,1) from random import uniform # used for accuracy tests FLOAT = np.float64 ctypedef np.float64_t FLOAT_T COMPLEX = np.complex128 ctypedef np.complex128_t COMPLEX_T cdef FLOAT_T PI = pi cdef FLOAT_T TWOPI = 2*PI cdef COMPLEX_T I = complex(0,1) cdef class Riemann_Map: """ The Riemann_Map class computes a Riemann or Ahlfors map from supplied data. It also has various methods to provide information about the map. A Riemann map conformally maps a simply connected region in The Riemann_Map class computes an interior or exterior Riemann map, or an Ahlfors map of a region given by the supplied boundary curve(s) and center ' point. The class also provides various methods to evaluate, visualize, or extract data from the map. A Riemann map conformally maps a simply connected region in the complex plane to the unit disc. The Ahlfors map does the same thing for multiply connected regions. Note that all the methods are numeric rather than analytic, for unusual regions or insufficient collocation points may give very inaccurate results. Note that all the methods are numerical. As a result all answers have some imprecision. Moreover, maps computed with small number of collocation points, or for unusually shaped regions may be very inaccurate. Error computations for the ellipse can be found in the documentation for analytic_boundary() and analytic_interior(). [BSV] provides an overview of the Riemann map and discusses the research that lead to the creation of this module. INPUT: Must be in the same order as fs. - a -- Complex, the center of the Riemann map. Will be mapped to the origin of the unit disc. the origin of the unit disc. Note that a MUST be within the region in order for the results to be mathematically valid. The following inputs must all be passed in as named parameters: - N -- integer (default: 500), the number of collocation points used to compute the map. More points will give more accurate results, especially near the boundaries, but will take longer to compute. - exterior -- boolean (default: False), if set to True, the exterior map will be computed, mapping the exterior of the region to the exterior of the unit circle. The following inputs may be passed as named parameters in unusual circumstances: - ncorners -- integer (default: 4), if mapping a figure with (equally t-spaced) corners, better results may be obtained by (equally t-spaced) corners--corners that make a significant change in the direction of the boundary, better results may be sometimes obtained by accurately giving this parameter. Used to add the proper constant to the theta correspondance function. the theta correspondence function. - opp -- boolean (default: False), set to True in very rare cases where the theta correspondance function is off by pi, that cases where the theta correspondence function is off by pi, that is, if red is mapped left of the origin in the color plot. EXAMPLES: .. [KT] N. Kerzman and M. R. Trummer. "Numerical Conformal Mapping via the Szego kernel". Journal of Computational and Applied Mathematics, 14(1-2): 111--123, 1986. .. [BSV] M. Bolt, S. Snoeyink, E. Van Andel. "Visual representation of the Riemann map and Ahlfors map via the Kerzman-Stein equation". Involve 3-4 (2010), 405-420. """ cdef int N, B, ncorners cdef f cdef opp cdef double complex a cdef np.ndarray tk, tk2, cps, dps, szego, p_vector, pre_q_vector cdef COMPLEX_T a cdef np.ndarray tk, tk2 cdef np.ndarray cps, dps, szego, p_vector, pre_q_vector cdef np.ndarray p_vector_inverse, sinalpha, cosalpha, theta_array cdef x_range, y_range cdef exterior def __init__(self, fs, fprimes, a, int N=500, int ncorners=4, opp=False): def __init__(self, fs, fprimes, COMPLEX_T a, int N=500, int ncorners=4, opp=False, exterior = False): """ Initializes the Riemann_Map class. See the class Riemann_Map for full documentation on the input of this initialization method. sage: f(t) = e^(I*t) - 0.5*e^(-I*t) sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t) sage: m = Riemann_Map([f], [fprime], 0) sage: m = Riemann_Map([f], [fprime], 0, N = 10) """ a = np.complex128(a) cdef double xmax, xmin, ymax, ymin, space cdef int i, k if N <= 2: raise ValueError( "The number of collocation points must be > 2.") if exterior and a!=0: raise ValueError("The exterior map requires a=0") try: fs = [fast_callable(f, domain=CDF) for f in fs] fprimes = [fast_callable(f, domain=CDF) for f in fprimes] self.ncorners = ncorners self.N = N  # Number of collocation pts self.opp = opp cdef int i, k self.exterior = exterior self.tk = np.array(np.arange(N) * TWOPI / N + 0.001 / N, dtype=np.float64) self.tk2 = np.zeros(N + 1, dtype=np.float64) dtype=FLOAT) self.tk2 = np.zeros(N + 1, dtype=FLOAT) for i in xrange(N): self.tk2[i] = self.tk[i] self.tk2[N] = TWOPI self.B = len(fs) # number of boundaries of the figure self.cps = np.zeros([self.B, N], dtype=np.complex128) self.dps = np.zeros([self.B, N], dtype=np.complex128) if exterior and (self.B > 1): raise ValueError( "The exterior map is undefined for multiply connected domains") cdef np.ndarray[COMPLEX_T,ndim=2] cps = np.zeros([self.B, N], dtype=COMPLEX) cdef np.ndarray[COMPLEX_T,ndim=2] dps = np.zeros([self.B, N], dtype=COMPLEX) # Find the points on the boundaries and their derivatives. for k in xrange(self.B): for i in xrange(N): self.cps[k, i] = np.complex(fs[k](self.tk[i])) self.dps[k, i] = np.complex(fprimes[k](self.tk[i])) cdef double xmax = self.cps.real.max() cdef double xmin = self.cps.real.min() cdef double ymax = self.cps.imag.max() cdef double ymin = self.cps.imag.min() cdef double space = 0.1 * max(xmax - xmin, ymax - ymin) #The default plotting window, mainly for color plot. if self.exterior: for k in xrange(self.B): for i in xrange(N): fk = fs[k](self.tk[N-i-1]) cps[k, i] = np.complex(1/fk) dps[k, i] = np.complex(1/fk**2*fprimes[k](self.tk[N-i-1])) else: for k in xrange(self.B): for i in xrange(N): cps[k, i] = np.complex(fs[k](self.tk[i])) dps[k, i] = np.complex(fprimes[k](self.tk[i])) if exterior: xmax = (1/cps).real.max() xmin = (1/cps).real.min() ymax = (1/cps).imag.max() ymin = (1/cps).imag.min() else: xmax = cps.real.max() xmin = cps.real.min() ymax = cps.imag.max() ymin = cps.imag.min() space = 0.1 * max(xmax - xmin, ymax - ymin) #The default plotting window for this map. self.cps = cps self.dps = dps self.x_range = (xmin - space, xmax + space) self.y_range = (ymin - space, ymax + space) # Now we do some more computation self._generate_theta_array() self._generate_interior_mapper() self._generate_inverse_mapper() def _repr_(self): """ Return a string representation of this Riemann_Map object. sage: f(t) = e^(I*t) - 0.5*e^(-I*t) sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t) sage: isinstance(Riemann_Map([f], [fprime], 0)._repr_(), str)  # long time sage: isinstance(Riemann_Map([f], [fprime], 0, N = 10)._repr_(), str)  # long time True """ return "A Riemann mapping of a figure to the unit circle." return "A Riemann or Ahlfors mapping of a figure to the unit circle." cdef _generate_theta_array(self): """ Generates the essential data for the Riemann map, primarily the Szego kernel and boundary correspondance. Szego kernel and boundary correspondence. See [KT] for the algorithm. TESTS:: sage: f(t) = e^(I*t) - 0.5*e^(-I*t) sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t) sage: m = Riemann_Map([f], [fprime], 0) sage: m = Riemann_Map([f], [fprime], 0, N = 10) """ cp = self.cps.flatten() dp = self.dps.flatten() cdef np.ndarray[COMPLEX_T,ndim =1] cp = self.cps.flatten() cdef np.ndarray[COMPLEX_T,ndim =1] dp = self.dps.flatten() cdef int N = self.N cdef int NB = N * self.B cdef int B = self.B cdef int i, k cdef FLOAT_T saa, t0 cdef np.ndarray[FLOAT_T, ndim=1] adp cdef np.ndarray[FLOAT_T, ndim=1] sadp cdef np.ndarray h cdef np.ndarray hconj cdef np.ndarray g cdef np.ndarray K cdef np.ndarray phi cdef np.ndarray[FLOAT_T, ndim=1] adp, sadp cdef np.ndarray[COMPLEX_T,ndim =1] h, hconj, g, normalized_dp,C, phi cdef np.ndarray[COMPLEX_T,ndim =2] K cdef np.ndarray[FLOAT_T, ndim=2] theta_array # Seting things up to use the Nystrom method adp = abs(dp) sadp = np.sqrt(adp) h = 1 / (TWOPI * I) * ((dp / adp) / (self.a - cp)) hconj = np.array(map(np.complex.conjugate, h), dtype=np.complex128) hconj = np.array(map(np.complex.conjugate, h), dtype=COMPLEX) g = -sadp * hconj normalized_dp=dp/adp C = I / N * sadp # equivalent to -TWOPI / N * 1 / (TWOPI * I) * sadp olderr = np.geterr()['invalid'] # checks the current error handling np.seterr(invalid='ignore') np.seterr(invalid='ignore') #there are some 0 divides that get over written later K = np.array( [C * sadp[t] * (normalized_dp/(cp-cp[t]) - (normalized_dp[t]/(cp-cp[t])).conjugate()) for t in np.arange(NB)], dtype=np.complex128) for t in np.arange(NB)], dtype=COMPLEX) np.seterr(invalid=olderr) # resets the error handling for i in xrange(NB): K[i, i] = 1 phi = np.linalg.solve(K, g) / NB * TWOPI  # Nystrom phi = np.linalg.solve(K, g) / NB * TWOPI  # Nystrom Method for solvin 2nd kind integrals # the all-important Szego kernel szego = np.array(phi.flatten() / np.sqrt(dp), dtype=np.complex128) szego = np.array(phi.flatten() / np.sqrt(dp), dtype=COMPLEX) self.szego = szego.reshape([B, N]) start = 0 # Finding the theta correspondance using phase. Misbehaves for some # Finding the theta correspondence using phase. Misbehaves for some # regions. if B != 1: theta_array = np.zeros([1, NB]) [theta_array.reshape([B, N]), np.zeros([B, 1])], axis=1) for k in xrange(B): self.theta_array[k, N] = self.theta_array[k, 0] + TWOPI # Finding the theta correspondance using ab. Well behaved, but # Finding the theta correspondence using abs. Well behaved, but # doesn't work on multiply connected domains. else: phi2 = phi.reshape([self.B, N]) sage: s = spline(points) sage: s(3*pi / 4) 0.00121587378429... sage: plot(s,0,2*pi) # plot the kernel The unit circle with a small hole:: """ Returns an array of points of the form [t value, theta in e^(I*theta)], that is, a discretized version of the theta/boundary correspondence function. For multiply connected domains, get_theta_points will list the points for each boundary in the order that they were supplied. of the theta/boundary correspondence function. In other words, a point in this array [t1, t2] represents that the boundary point given by f(t1) is mapped to a point on the boundary of the unit circle given by e^(I*t2). For multiply connected domains, get_theta_points will list the points for each boundary in the order that they were supplied. INPUT: The following input must all be passed in as named parameters: - boundary -- integer (default: -1) if < 0, - boundary -- integer (default: -1) if < 0, get_theta_points() will return the points for all boundaries. If >= 0, get_theta_points() will return only the points for the boundary specified. Extending the points by a spline:: sage: s = spline(points) sage: s(3*pi / 4) sage: s(3*pi / 4) 1.62766037996... The unit circle with a small hole:: sage: hfprime(t) = 0.5*-I*e^(-I*t) sage: m = Riemann_Map([f, hf], [hf, hfprime], 0.5 + 0.5*I) Getting the szego for a specifc boundary:: Getting the boundary correspondence for a specifc boundary:: sage: tp0 = m.get_theta_points(boundary=0) sage: tp1 = m.get_theta_points(boundary=1) cdef _generate_interior_mapper(self): """ Generates the data necessary to use the reimann_map() function. As much setup as possible is done here to minimize what has to be done in riemann_map(). Generates the data necessary to use the riemann_map function. As much setup as possible is done here to minimize the computation that must be done in riemann_map. TESTS:: sage: f(t) = e^(I*t) - 0.5*e^(-I*t) sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t) sage: m = Riemann_Map([f], [fprime], 0) sage: m = Riemann_Map([f], [fprime], 0, N = 10) """ cdef int N = self.N cdef double complex coeff = 3*I / (8*N) cdef np.ndarray[double complex, ndim=1] pq = self.cps[:,list(range(N))+[0]].flatten() self.pre_q_vector = pq cpdef riemann_map(self, pt): cpdef riemann_map(self, COMPLEX_T pt): """ Returns the Riemann mapping of a point. That is, given pt on the interior of the mapped region, riemann_map() will return the interior of the mapped region, riemann_map will return the point on the unit disk that pt maps to. Note that this method only works for interior points; it breaks down very close method only works for interior points; accuracy breaks down very close to the boundary. To get boundary corrospondance, use get_theta_points(). sage: m.riemann_map(np.complex(-3, 0.0001)) (1.405757...e-05+8.06...e-10j) """ pt1 = np.complex(pt) cdef np.ndarray[double complex, ndim=1] q_vector = 1 / ( self.pre_q_vector - pt1) return -np.dot(self.p_vector, q_vector) cdef COMPLEX_T pt1 cdef np.ndarray[COMPLEX_T, ndim=1] q_vector if self.exterior: pt1 = 1/pt q_vector = 1 / ( self.pre_q_vector - pt1) return -1/np.dot(self.p_vector, q_vector) else: pt1 = pt q_vector = 1 / ( self.pre_q_vector - pt1) return -np.dot(self.p_vector, q_vector) cdef _generate_inverse_mapper(self): """ Generates the data necessary to use the inverse_reimann_map() function. As much setup as possible is done here to minimize what has to be done in inverse_riemann_map() function. As much setup as possible is done here to minimize the computation that must be done in inverse_riemann_map(). TESTS:: sage: f(t) = e^(I*t) - 0.5*e^(-I*t) sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t) sage: m = Riemann_Map([f], [fprime], 0) sage: m = Riemann_Map([f], [fprime], 0, N = 10) """ cdef int N = self.N cdef int B = self.B for i in xrange(N): self.cosalpha[k, i] = cos(-theta_array[k, i]) def inverse_riemann_map(self, pt): cpdef inverse_riemann_map(self, COMPLEX_T pt): """ Returns the inverse Riemann mapping of a point. That is, given pt on the interior of the unit disc, inverse_reimann_map() will on the interior of the unit disc, inverse_riemann_map() will return the point on the original region that would be Riemann mapped to pt. .. NOTE:: This method does not work for multiply connected domains. mapped to pt. Note that this method does not work for multiply connected domains. INPUT: sage: m.inverse_riemann_map(np.complex(-0.2, 0.5)) (-0.156280570579...+0.321819151891...j) """ pt = np.complex128(pt) if self.exterior: pt = 1/pt r = abs(pt) if r == 0: stheta = 0 stheta = pt.imag / r ctheta = pt.real / r k = 0 return (1 - r**2) / TWOPI * np.dot( mapped = (1 - r**2) / TWOPI * np.dot( self.p_vector_inverse[k] * self.cps[k], 1 / (1 + r**2 - 2*r * (ctheta * self.cosalpha[k] - stheta * self.sinalpha[k]))) if self.exterior: return 1/mapped else: return mapped def plot_boundaries(self, plotjoined=True, rgbcolor=[0,0,0], thickness=1): """ """ plots = range(self.B) for k in xrange(self.B): # This should be eliminated when the thickness/pointsize issue # This conditional should be eliminated when the thickness/pointsize issue # is resolved later. Same for the others in plot_spiderweb(). if plotjoined: plots[k] = list_plot( comp_pt(self.cps[k], 1), rgbcolor=rgbcolor, pointsize=thickness) return sum(plots) cpdef compute_on_grid(self, plot_range, int x_points): """ Computes the riemann map on a grid of points. Note that these points are complex of the form z = x + y*i. INPUT: - plot_range -- a tuple of the form [xmin, xmax, ymin, ymax]. If the value is [], the default plotting window of the map will be used. - x_points -- int, the size of the grid in the x direction The number of points in the y_direction is scaled accordingly OUTPUT: - a tuple containing [z_values, xmin, xmax, ymin, ymax] where z_values is the evaluation of the map on the specified grid. EXAMPLES: General usage:: sage: f(t) = e^(I*t) - 0.5*e^(-I*t) sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t) sage: m = Riemann_Map([f], [fprime], 0) sage: data = m.compute_on_grid([],5) sage: print data[0][8,1] (-0.0879...+0.9709...j) """ cdef FLOAT_T xmin, xmax, xstep, ymin, ymax, ystep cdef int y_points if plot_range == []: xmin = self.x_range[0] xmax = self.x_range[1] ymin = self.y_range[0] ymax = self.y_range[1] else: xmin = plot_range[0] xmax = plot_range[1] ymin = plot_range[2] ymax = plot_range[3] xstep = (xmax - xmin) / x_points ystep = xstep y_points = int((ymax-ymin)/ ystep) cdef Py_ssize_t i, j cdef COMPLEX_T pt cdef np.ndarray[COMPLEX_T] pre_q_vector = self.pre_q_vector cdef np.ndarray[COMPLEX_T] p_vector = self.p_vector cdef np.ndarray[COMPLEX_T, ndim=2] z_values = np.empty( [y_points, x_points], dtype=np.complex128) if self.exterior: for i in xrange(x_points): for j in xrange(y_points): pt = 1/(xmin + 0.5*xstep + i*xstep + I*(ymin + 0.5*ystep + j*ystep)) z_values[j, i] = 1/(-np.dot(p_vector,1/(pre_q_vector - pt))) else: for i in xrange(x_points): for j in xrange(y_points): pt = xmin + 0.5*xstep + i*xstep + I*(ymin + 0.5*ystep + j*ystep) z_values[j, i] = -np.dot(p_vector,1/(pre_q_vector - pt)) return z_values, xmin, xmax, ymin, ymax @options(interpolation='catrom') def plot_spiderweb(self, spokes=16, circles=4, pts=32, linescale=0.99, rgbcolor=[0,0,0], thickness=1, plotjoined=True): rgbcolor=[0,0,0], thickness=1, plotjoined=True, withcolor = False,plot_points = 200, **options): """ Generates a traditional "spiderweb plot" of the Riemann map. Shows what concentric circles and radial lines map to. Note that this method DOES NOT work for multiply connected domains. Generates a traditional "spiderweb plot" of the Riemann map. Shows what concentric circles and radial lines map to. The radial lines may exhibit erratic behavior near the boundary, if this occurs, decreasing linescale may mitigate the problem. Note that this method requires significantly more computation for multiply connected domains. INPUT: plot. Each radial line is made by 1*pts points, each circle has 2*pts points. Note that high values may cause erratic behavior of the radial lines near the boundaries. - only for simply connected domains - linescale -- float between 0 and 1. Shrinks the radial lines away from the boundary to reduce erratic behavior. - only for simply connected domains - rgbcolor -- float array (default: [0,0,0]) the red-green-blue color of the spiderweb. - plotjoined -- boolean (default: True) If False, discrete points will be drawn; otherwise they will be connected by lines. - only for simply connected domains - withcolor -- boolean (default: False) If True, The spiderweb will be overlaid on the basic color plot. - plot_points -- integer (default: 200) the size of the grid in the x direction The number of points in the y_direction is scaled accordingly. Note that very large values can cause this function to run slowly. - only for multiply connected domains EXAMPLES: General usage:: sage: m = Riemann_Map([f], [fprime], 0, 1000) sage: m.plot_spiderweb() """ edge = self.plot_boundaries(plotjoined=plotjoined, rgbcolor=rgbcolor, thickness=thickness) circle_list = range(circles) theta_array = self.theta_array[0] s = spline(np.column_stack([self.theta_array[0], self.tk2]).tolist()) tmax = self.theta_array[0, self.N] tmin = self.theta_array[0, 0] cdef int k, i for k in xrange(circles): temp = range(pts*2) for i in xrange(2*pts): temp[i] = self.inverse_riemann_map( (k + 1) / (circles + 1.0) * exp(I*i * TWOPI / (2*pts))) if plotjoined: circle_list[k] = list_plot( comp_pt(temp, 1), rgbcolor=rgbcolor, thickness=thickness, plotjoined=True) if self.B == 1: #The efficient simply connected edge = self.plot_boundaries(plotjoined=plotjoined, rgbcolor=rgbcolor, thickness=thickness) circle_list = range(circles) theta_array = self.theta_array[0] s = spline(np.column_stack([self.theta_array[0], self.tk2]).tolist()) tmax = self.theta_array[0, self.N] tmin = self.theta_array[0, 0] for k in xrange(circles): temp = range(pts*2) for i in xrange(2*pts): temp[i] = self.inverse_riemann_map( (k + 1) / (circles + 1.0) * exp(I*i * TWOPI / (2*pts))) if plotjoined: circle_list[k] = list_plot( comp_pt(temp, 1), rgbcolor=rgbcolor, thickness=thickness, plotjoined=True) else: circle_list[k] = list_plot( comp_pt(temp, 1), rgbcolor=rgbcolor, pointsize=thickness) line_list = range(spokes) for k in xrange(spokes): temp = range(pts) angle = (k*1.0) / spokes * TWOPI if angle >= tmax: angle -= TWOPI elif angle <= tmin: angle += TWOPI for i in xrange(pts - 1): temp[i] = self.inverse_riemann_map( (i * 1.0) / (pts * 1.0) * exp(I * angle) * linescale) temp[pts - 1] = np.complex( self.f(s(angle)) if angle <= tmax else self.f(s(angle-TWOPI))) if plotjoined: line_list[k] = list_plot( comp_pt(temp, 0), rgbcolor=rgbcolor, thickness=thickness, plotjoined=True) else: line_list[k] = list_plot( comp_pt(temp, 0), rgbcolor=rgbcolor, pointsize=thickness) if withcolor: return edge + sum(circle_list) + sum(line_list) + self.plot_colored(plot_points=plot_points) else: circle_list[k] = list_plot( comp_pt(temp, 1), rgbcolor=rgbcolor, pointsize=thickness) line_list = range(spokes) for k in xrange(spokes): temp = range(pts) angle = (k*1.0) / spokes * TWOPI if angle >= tmax: angle -= TWOPI elif angle <= tmin: angle += TWOPI for i in xrange(pts - 1): temp[i] = self.inverse_riemann_map( (i * 1.0) / (pts * 1.0) * exp(I * angle) * linescale) temp[pts - 1] = np.complex( self.f(s(angle)) if angle <= tmax else self.f(s(angle-TWOPI))) if plotjoined: line_list[k] = list_plot( comp_pt(temp, 0), rgbcolor=rgbcolor, thickness=thickness, plotjoined=True) else: line_list[k] = list_plot( comp_pt(temp, 0), rgbcolor=rgbcolor, pointsize=thickness) return edge + sum(circle_list) + sum(line_list) return edge + sum(circle_list) + sum(line_list) else: # The more difficult multiply connected z_values, xmin, xmax, ymin, ymax = self.compute_on_grid([], plot_points) xstep = (xmax-xmin)/plot_points ystep = (ymax-ymin)/plot_points dr, dtheta= get_derivatives(z_values, xstep, ystep) # clean later g = Graphics() g.add_primitive(ComplexPlot(complex_to_spiderweb(z_values,dr,dtheta, spokes, circles, rgbcolor,thickness, withcolor), (xmin, xmax), (ymin, ymax),options)) return g + self.plot_boundaries(thickness = thickness) @options(interpolation='catrom') def plot_colored(self, plot_range=[], int plot_points=100, **options): """ Draws a colored plot of the Riemann map. A red point on the colored plot corresponds to a red point on the unit disc. Note that this method DOES work for multiply connected domains. Generates a colored plot of the Riemann map. A red point on the colored plot corresponds to a red point on the unit disc. INPUT: (xmin, xmax, ymin, ymax). Declare if you do not want the plot to use the default range for the figure. - plot_points -- integer (default: 100), number of points to plot in each direction of the grid. Note that very large values can cause this function to run slowly. - plot_points -- integer (default: 100), number of points to plot in the x direction. Points in the y direction are scaled accordingly. Note that very large values can cause this function to run slowly. EXAMPLES: sage: m = Riemann_Map([f], [fprime], 0, 1000) sage: m.plot_colored() """ cdef int i, j cdef double xmin cdef double xmax cdef double ymin cdef double ymax if plot_range == []: xmin = self.x_range[0] xmax = self.x_range[1] ymin = self.y_range[0] ymax = self.y_range[1] else: xmin = plot_range[0] xmax = plot_range[1] ymin = plot_range[2] ymax = plot_range[3] xstep = (xmax - xmin) / plot_points ystep = (ymax - ymin) / plot_points cdef np.ndarray z_values = np.empty( [plot_points, plot_points], dtype=np.complex128) for i in xrange(plot_points): for j in xrange(plot_points): z_values[j, i] = self.riemann_map( np.complex(xmin + 0.5*xstep + i*xstep, ymin + 0.5*ystep + j*ystep)) z_values, xmin, xmax, ymin, ymax = self.compute_on_grid(plot_range, plot_points) g = Graphics() g.add_primitive(ComplexPlot(complex_to_rgb(z_values), (xmin, xmax), (ymin, ymax),options)) return g cdef comp_pt(clist, loop=True): """ This function converts a list of complex numbers to the plottable Converts a list of complex numbers xderivs = get_derivatives(z_values, xstep, ystep)[0]to the plottable (x,y) form. If loop=True, then the first point will be added as the last, i.e. to plot a closed circle. if loop: list2[len(clist)] = list2[0] return list2 cpdef get_derivatives(np.ndarray[COMPLEX_T, ndim=2]z_values, FLOAT_T xstep, FLOAT_T ystep): """ Computes the r*e^(I*theta) form derivatives from the grid of points. The derivatives are computed using quick-and-dirty taylor expansion and assuming analyticy. As such get_derivatives is primarily intended to be used for comparisions in plot_spiderweb and not for applications that require great precision. INPUT: - z_values -- The values for a complex function evaluated on a grid in the complex plane, usually from compute_on_grid. - xstep -- float, the spacing of the grid points in the real direction OUTPUT: - A tuple of arrays, [dr, dtheta], each array is 2 less in both dimensions than z_values dr - the absolute value of the derivative of the function in the +r direction dtheta - the rate of accumulation of angle in the +theta direction EXAMPLES: Standard usage with compute_on_grid:: sage: from sage.calculus.riemann import get_derivatives sage: f(t) = e^(I*t) - 0.5*e^(-I*t) sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t) sage: m = Riemann_Map([f], [fprime], 0) sage: data = m.compute_on_grid([],19) sage: xstep = (data[2]-data[1])/19 sage: ystep = (data[4]-data[3])/19 sage: dr, dtheta = get_derivatives(data[0],xstep,ystep) sage: dr[8,8] 0.241... sage: dtheta[5,5] 5.907... cdef inline double mag_to_lightness(double r): """ Tweak this to adjust how the magnitude affects the color. Note this method is customized for riemann plots, the magnitude loops rather than fading to black. cdef np.ndarray[COMPLEX_T, ndim=2] xderiv cdef np.ndarray[FLOAT_T, ndim = 2] dr, dtheta, zabs imax = len(z_values)-2 jmax = len(z_values[0])-2 xderiv = (z_values[1:-1,2:]-z_values[1:-1,:-2])/(2*xstep) #(f(x+delta)-f(x-delta))/2delta dr = np.abs(xderiv) #b/c the function is analytic, we know it's abs(derivative) is equal in all directions zabs = np.abs(z_values[1:-1,1:-1]) # the abs(derivative) scaled by distance from origin dtheta = np.divide(dr,zabs) return dr, dtheta cpdef complex_to_spiderweb(np.ndarray[COMPLEX_T, ndim = 2]z_values, np.ndarray[FLOAT_T, ndim = 2] dr, np.ndarray[FLOAT_T, ndim = 2] dtheta, spokes, circles, rgbcolor, thickness, withcolor): """ Converts a grid of complex numbers into a matrix containing rgb data for the Riemann spiderweb plot. INPUT: - r -- a non-negative real number. - z_values -- A grid of complex numbers, as a list of lists. - dr -- grid of floats, the r derivative of z_values. Used to determine precision. - dtheta -- grid of floats, the theta derivative of z_values. Used to determine precision. - spokes -- integer the number of equallyspaced radial lines to plot. - circles -- integer the number of equally spaced circles about the center to plot. - rgbcolor -- float array the red-green-blue color of the lines of the spiderweb. - thickness -- positive float the thickness of the lines or points in the spiderweb. - withcolor -- boolean If True the spiderweb will be overlaid on the basic color plot. OUTPUT: A value between -1 (black) and +1 (white), inclusive. An N x M x 3 floating point Numpy array X, where X[i,j] is an (r,g,b) tuple. EXAMPLES: sage: from sage.calculus.riemann import complex_to_spiderweb sage: import numpy sage: zval = numpy.array([[0, 1, 1000],[.2+.3j,1,-.3j],[0,0,0]],dtype = numpy.complex128) sage: deriv = numpy.array([[.1]],dtype = numpy.float64) sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,False) array([[[ 1.,  1.,  1.], [ 1.,  1.,  1.], [ 1.,  1.,  1.]], [[ 1.,  1.,  1.], [ 0.,  0.,  0.], [ 1.,  1.,  1.]], [[ 1.,  1.,  1.], [ 1.,  1.,  1.], [ 1.,  1.,  1.]]]) EXAMPLES: sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,True) array([[[ 1.        ,  1.        ,  1.        ], [ 1.        ,  0.05558355,  0.05558355], [ 0.17301243,  0.        ,  0.        ]], [[ 1.        ,  0.96804683,  0.48044583], [ 0.        ,  0.        ,  0.        ], [ 0.77351965,  0.5470393 ,  1.        ]], [[ 1.        ,  1.        ,  1.        ], [ 1.        ,  1.        ,  1.        ], [ 1.        ,  1.        ,  1.        ]]]) This tests it implicitly:: sage: from sage.calculus.riemann import complex_to_rgb sage: import numpy sage: complex_to_rgb(numpy.array([[0, 1, 1000]],dtype = numpy.complex128)) array([[[   1.,    1.,    1.], [   1.,    0.,    0.], [-998.,    0.,    0.]]]) """ return 1 - r cpdef complex_to_rgb(np.ndarray z_values): r""" Convert from an array of complex numbers to its corresponding matrix of cdef Py_ssize_t i, j, imax, jmax cdef FLOAT_T x, y, mag, arg, width, target, precision, dmag, darg cdef COMPLEX_T z cdef FLOAT_T DMAX = 70 # change to adjust rate_of_change cutoff below cdef FLOAT_T MMIN = .0001 precision = thickness/150.0 imax = len(z_values) jmax = len(z_values[0]) cdef np.ndarray[FLOAT_T, ndim=3, mode="c"] rgb if withcolor: rgb = complex_to_rgb(z_values) else: rgb = np.zeros(dtype=FLOAT, shape=(imax, jmax, 3)) rgb += 1 if circles != 0: circ_radii = srange(0,1.0,1.0/circles) else: circ_radii = [] if spokes != 0: spoke_angles = srange(-PI,PI+TWOPI/spokes,TWOPI/spokes) # both -pi and pi are included else: spoke_angles = [] for i in xrange(imax-2): # the d arrays are 1 smaller on each side for j in xrange(jmax-2): z = z_values[i+1,j+1] mag = abs(z) arg = phase(z) dmag = dr[i,j] darg = dtheta[i,j] if darg < DMAX and mag > MMIN: #points that change too rapidly are presumed to be borders, points that are too small are presumed to be outside for target in circ_radii: if abs(mag - target)/dmag < precision: rgb[i+1,j+1] = rgbcolor break for target in spoke_angles: if abs(arg - target)/darg < precision: rgb[i+1,j+1] = rgbcolor break return rgb cpdef complex_to_rgb(np.ndarray[COMPLEX_T, ndim = 2] z_values): """ Converts an array of complex numbers to its corresponding matrix of RGB values. INPUT: OUTPUT: An N \times M \times 3 floating point Numpy array X, where An N x M x 3 floating point Numpy array X, where X[i,j] is an (r,g,b) tuple. EXAMPLES:: sage: from sage.calculus.riemann import complex_to_rgb sage: import numpy sage: complex_to_rgb(numpy.array([[0, 1, 1000]],dtype = numpy.complex128)) array([[[   1.,    1.,    1.], [   1.,    0.,    0.], [-998.,    0.,    0.]]]) array([[[ 1.        ,  1.        ,  1.        ], [ 1.        ,  0.05558355,  0.05558355], [ 0.17301243,  0.        ,  0.        ]]]) sage: complex_to_rgb(numpy.array([[0, 1j, 1000j]],dtype = numpy.complex128)) array([[[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00], [  5.00000000e-01,   1.00000000e+00,   0.00000000e+00], [ -4.99000000e+02,  -9.98000000e+02,   0.00000000e+00]]]) array([[[ 1.        ,  1.        ,  1.        ], [ 0.52779177,  1.        ,  0.05558355], [ 0.08650622,  0.17301243,  0.        ]]]) """ cdef unsigned int i, j, imax, jmax cdef double x, y, mag, arg cdef double lightness, hue, top, bot cdef double r, g, b cdef Py_ssize_t i, j, imax, jmax cdef FLOAT_T x, y, mag, arg cdef FLOAT_T lightness, hue, top, bot cdef FLOAT_T r, g, b cdef int ihue cdef z from cmath import phase cdef COMPLEX_T z imax = len(z_values) jmax = len(z_values[0]) cdef np.ndarray[np.float_t, ndim=3, mode="c"] rgb = np.empty( dtype=np.float, shape=(imax, jmax, 3)) cdef np.ndarray[FLOAT_T, ndim=3, mode="c"] rgb = np.empty( dtype=FLOAT, shape=(imax, jmax, 3)) sig_on() for i from 0 <= i < imax: for i from 0 <= i < imax: #replace with xrange? row = z_values[i] for j from 0 <= j < jmax: for j from 0 <= j < jmax: # replace with xrange? z = row[j] mag = abs(z) arg = phase(z) lightness = mag_to_lightness(mag) # tweak these levels to adjust how bright/dark the colors appear # output can range from -1 (black) to 1 (white) lightness = -(atan(log(mag*1.5+1)) * (4/PI) - 1) # in hsv, variable value, full saturation (s=1, v=1+lightness) if lightness < 0: bot = 0 rgb[i, j, 2] = b sig_off() return rgb cpdef analytic_boundary(FLOAT_T t, int n): """ Provides an exact (for n = infinity) riemann boundary correspondence for the ellipse with axes 1 + epsilon and 1 - epsilon. The boundary is therefore given by e^(I*t)+epsilon*e^(-I*t). It is primarily useful for testing the accuracy of the numerical Riemann Map. INPUT: - t -- The boundary parameter, from 0 to two*pi - n -- Integer, The number of terms to include. 10 is fairly accurate, 20 is very accurate. OUTPUT: A theta value from 0 to 2*pi, corresponding to the point on the circle e^(I*theta) TESTS: Checking the accuracy for different n values:: sage: from sage.calculus.riemann import analytic_boundary sage: t100 = analytic_boundary(pi/2,100) sage: abs(analytic_boundary(pi/2,10) - t100) < 10^-8 True sage: abs(analytic_boundary(pi/2,20) - t100) < 10^-15 True Using this to check the accuracy of the Riemann_Map boundary:: sage: f(t) = e^(I*t)+.3*e^(-I*t) sage: fp(t) = I*e^(I*t)-I*.3*e^(-I*t) sage: m = Riemann_Map([f], [fp],0,200) sage: s = spline(m.get_theta_points()) sage: test_pt = uniform(0,2*pi) sage: s(test_pt) - analytic_boundary(test_pt,20) < 10^-4 True """ cdef FLOAT_T i cdef FLOAT_T result = t for i from 1 <= i < n+1: result += (2*(-1)**i/i)*(.3**i/(1+.3**(2*i)))*sin(2*i*t) return result cpdef cauchy_kernel(t, args): """ Intermediate function for the integration in analytic_interior. INPUT: - t -- The boundary parameter, meant to be integrated over - args -- a tuple containing: - z -- Complex, the point to be mapped. - n -- Integer, The number of terms to include. 10 is fairly accurate, 20 is very accurate. - part -- will return the real ('r'), imaginary ('i') or complex ('c') value of the kernel TESTS: Primarily tested implicitly by analytic_interior, Simple test:: sage: from sage.calculus.riemann import cauchy_kernel sage: cauchy_kernel(.5,(.1+.2*I, 10,'c')) (-0.584136405997...+0.5948650858950...j) """ cdef COMPLEX_T result cdef COMPLEX_T z = args[0] cdef int n = args[1] part = args[2] result = exp(I*analytic_boundary(t,n))/(exp(I*t)+.3*exp(-I*t)-z)*(I*exp(I*t)-I*.3*exp(-I*t)) if part == 'c': return result elif part == 'r': return result.real elif part == 'i': return result.imag else: return None cpdef analytic_interior(COMPLEX_T z, int n): """ Provides a nearly exact compuation of the Riemann Map of an interior point of the ellipse with axes 1 + epsilon and 1 - epsilon. It is primarily useful for testing the accuracy of the numerical Riemann Map. INPUT: - z -- Complex, the point to be mapped. - n -- Integer, The number of terms to include. 10 is fairly accurate, 20 is very accurate. TESTS: Testing the accuracy of Riemann_Map:: sage: from sage.calculus.riemann import analytic_interior sage: f(t) = e^(I*t)+.3*e^(-I*t) sage: fp(t) = I*e^(I*t)-I*.3*e^(-I*t) sage: m = Riemann_Map([f],[fp],0,200) sage: abs(m.riemann_map(.5)-analytic_interior(.5,20)) < 10^-4 True sage: m = Riemann_Map([f],[fp],0,2000) sage: abs(m.riemann_map(.5)-analytic_interior(.5,20)) < 10^-6 True """ # evaluates the cauchy integral of the boundary, split into the real # and imaginary results because numerical_integral can't handle complex data. rp = 1/(TWOPI)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'i'])[0] ip = 1/(TWOPI*I)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'r'])[0] return rp + ip