# Ticket #11273: trac-11273-riemann_upgrade_part3.patch

File trac-11273-riemann_upgrade_part3.patch, 22.1 KB (added by evanandel, 10 years ago)

Formatting Fixes

• ## sage/calculus/riemann.pyx

```# HG changeset patch
# User Ethan Van Andel <evlutte@gmail.com>
# Date 1311270131 14400
# Node ID eb019dc2d8a6f2db9f1ae66377d38eb6025eb08a
# Parent  89043440f58380e7007e30d0e991351c3d2224c3
Trac 11273: additional patch, formatting fixes

diff -r 89043440f583 -r eb019dc2d8a6 sage/calculus/riemann.pyx```
 a cdef class Riemann_Map: """ 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. 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 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()``. 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: - ``fs`` -- A list of the boundaries of the region, given as complex-valued functions with domain ``0`` to ``2*pi``. Note that the complex-valued functions with domain `0` to `2*pi`. Note that the outer boundary must be parameterized counter clockwise (i.e. ``e^(I*t)``) while the inner boundaries must be clockwise (i.e. ``e^(-I*t)``). 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: 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--corners that make a significant change in cdef x_range, y_range cdef exterior def __init__(self, fs, fprimes, COMPLEX_T a, int N=500, int ncorners=4, opp=False, exterior = 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. 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) 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. if self.exterior: for k in xrange(self.B): 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') #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()) #there are some 0 divides that get over written later np.seterr(invalid='ignore') 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=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 Method for solvin 2nd kind integrals # Nystrom Method for solvin 2nd kind integrals phi = np.linalg.solve(K, g) / NB * TWOPI # the all-important Szego kernel szego = np.array(phi.flatten() / np.sqrt(dp), dtype=COMPLEX) self.szego = szego.reshape([B, N]) 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_riemann_map()`` will return the point on the original region that would be Riemann mapped to ``pt``. Note that this method does not work for multiply connected domains. 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 that this method does not work for multiply connected domains. INPUT: """ plots = range(self.B) for k in xrange(self.B): # This conditional should be eliminated when the thickness/pointsize issue # is resolved later. Same for the others in plot_spiderweb(). # 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), plotjoined=True, OUTPUT: - a tuple containing [z_values, xmin, xmax, ymin, ymax] where z_values is the evaluation of the map on the specified grid. - a tuple containing [z_values, xmin, xmax, ymin, ymax] where z_values is the evaluation of the map on the specified grid. EXAMPLES: @options(interpolation='catrom') def plot_spiderweb(self, spokes=16, circles=4, pts=32, linescale=0.99, rgbcolor=[0,0,0], thickness=1, plotjoined=True, withcolor = False,plot_points = 200, **options): 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. The radial lines """ cdef int k, i if self.B == 1: #The efficient simply connected edge = self.plot_boundaries(plotjoined=plotjoined, rgbcolor=rgbcolor, thickness=thickness) edge = self.plot_boundaries(plotjoined=plotjoined, rgbcolor=rgbcolor, thickness=thickness) circle_list = range(circles) theta_array = self.theta_array s = spline(np.column_stack([self.theta_array, self.tk2]).tolist()) 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) 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) 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) 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) return edge + sum(circle_list) + sum(line_list) + \ self.plot_colored(plot_points=plot_points) else: 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) 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)) 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) ``(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 the x direction. Points in the y direction are scaled accordingly. 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() """ z_values, xmin, xmax, ymin, ymax = self.compute_on_grid(plot_range, plot_points) 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)) g.add_primitive(ComplexPlot(complex_to_rgb(z_values), (xmin, xmax), (ymin, ymax),options)) return g cdef comp_pt(clist, loop=True): list2[len(clist)] = list2 return list2 cpdef get_derivatives(np.ndarray[COMPLEX_T, ndim=2]z_values, FLOAT_T xstep, FLOAT_T ystep): 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. 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``. - ``z_values`` -- The values for a complex function evaluated on a grid in the complexplane, 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 - A tuple of arrays, [``dr``, ``dtheta``], each array is 2 less in both dimensions than ``z_values`` ``dr`` - the abs 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) cdef np.ndarray[FLOAT_T, ndim = 2] dr, dtheta, zabs imax = len(z_values)-2 jmax = len(z_values)-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 #(f(x+delta)-f(x-delta))/2delta xderiv = (z_values[1:-1,2:]-z_values[1:-1,:-2])/(2*xstep) #b/c the function is analytic, we know it's abs(derivative) is equal #in all directions dr = np.abs(xderiv) # the abs(derivative) scaled by distance from origin zabs = np.abs(z_values[1:-1,1:-1]) 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): 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. - ``spokes`` -- integer the number of equallyspaced radial lines to plot. - ``circles`` -- integer the number of equally spaced circles about the center 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. - ``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. - ``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. - ``withcolor`` -- boolean If ``True`` the spiderweb will be overlaid on the basic color plot. OUTPUT: An `N x M x 3` floating point Numpy array ``X``, where ``X[i,j]`` is an (r,g,b) tuple. EXAMPLES: 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) else: circ_radii = [] if spokes != 0: spoke_angles = srange(-PI,PI+TWOPI/spokes,TWOPI/spokes) # both -pi and pi are included # both -pi and pi are included spoke_angles = srange(-PI,PI+TWOPI/spokes,TWOPI/spokes) else: spoke_angles = [] for i in xrange(imax-2): # the d arrays are 1 smaller on each side 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 #points that change too rapidly are presumed to be borders, #points that are too small are presumed to be outside if darg < DMAX and mag > MMIN: for target in circ_radii: if abs(mag - target)/dmag < precision: rgb[i+1,j+1] = rgbcolor ``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)) sage: complex_to_rgb(numpy.array([[0, 1, 1000]], dtype = numpy.complex128)) 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)) sage: complex_to_rgb(numpy.array([[0, 1j, 1000j]], dtype = numpy.complex128)) array([[[ 1.        ,  1.        ,  1.        ], [ 0.52779177,  1.        ,  0.05558355], [ 0.08650622,  0.17301243,  0.        ]]]) """ 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. 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: OUTPUT: A theta value from 0 to 2*pi, corresponding to the point on the circle e^(I*theta) 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 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) - ``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. - ``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 - ``part`` -- will return the real ('r'), imaginary ('i') or complex ('c') value of the kernel TESTS: Primarily tested implicitly by analytic_interior, 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 z = args cdef int n = args part = args result = exp(I*analytic_boundary(t,n))/(exp(I*t)+.3*exp(-I*t)-z)*(I*exp(I*t)-I*.3*exp(-I*t)) 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': 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. 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: 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) """ # 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']) ip = 1/(TWOPI*I)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'r']) rp = 1/(TWOPI)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'i']) ip = 1/(TWOPI*I)*numerical_integral(cauchy_kernel,0,2*pi, params = [z,n,'r']) return rp + ip
• ## sage/plot/complex_plot.pyx

`diff -r 89043440f583 -r eb019dc2d8a6 sage/plot/complex_plot.pyx`
 a """ self.xrange = xrange self.yrange = yrange #self.z_values = z_values self.x_count = len(rgb_data) self.y_count = len(rgb_data) self.rgb_data = rgb_data