Thanks for the help so far. I've come a fair way and implemented a good portion in cython. However, the speedup was only minor (< a factor of 2) and I suspect I'm missing something. If anyone can give me any tips on how to optimize the code, that would be most appreciated. Tips on the code itself or general good-practice would be great. At this point, from reading the cython user manual and numpy tutorial, I guess that I should do the following, but do not know how.
1: make the whole class into a struct or as C'ish as possible 2: ctypedef the np.complex128 type or something. 3: cdef my functions so that they run faster. Unfortunately I'm a bit stuck on how to do those things, so any clarification on those would be great. A brief description of what my code does: you create a Riemann_Map class by giving it a function that defines the boundary of a complex figure in the plane (like a circle or ellipse), the derivative of said function, a center point, and the number of points you want to use in the computation. Riemann_Map class maps your figure to the unit circle, giving you a boundary correspondence (theta_array) and function to map points on the figure's interior to the circle's interior (riemann_map(pt)). That's a poor, general description. These articles might give more info: An explanation of the theory can be found here: http://www.springerlink.com/content/r44467137l65h5x1/ I think this paper might have some partial examples. http://portal.acm.org/citation.cfm?id=901698 Thank you for any help you may be able to provide, Also, If this has left the realms of stuff that can reasonable be posted in sage support, I apologize and please direct me elsewhere. Here's my code: %cython %time import numpy as np cimport numpy as np import math #cimport math import cmath #cimport cmath from math import pi from cmath import exp cdef double PI = pi cdef double TWOPI = 2*PI cdef I = complex(0,1) class Riemann_Map(): def __init__(self,f,fprime,a,N, init_full = False): self.a = a self.N = N self.tk = np.array(np.arange(N)*TWOPI/N, dtype = np.float64) cp = np.zeros(N,dtype = np.complex128) dp = np.zeros(N,dtype = np.complex128) cdef int i for i in xrange(N): cp[i] = f(self.tk[i]) dp[i] = fprime(self.tk[i]) self.cp = cp self.dp = dp self.theta_array = None self.p_vector = None if init_full: self.generate_theta_array() self.generate_interior_mapper() pass def testf(self): print "N" print self.N print "tk" print self.tk print "cp" print self.cp print "dp" print self.dp print "theta_array" print self.theta_array def generate_theta_array(self): cdef np.ndarray cp = self.cp cdef np.ndarray dp = self.dp cdef unsigned N = self.N cdef np.ndarray adp cdef np.ndarray sadp cdef np.ndarray h cdef np.ndarray hconj cdef np.ndarray g cdef np.ndarray K cdef np.ndarray phi cdef np.ndarray temp adp = abs(dp) sadp = np.array(map(math.sqrt,adp),dtype = np.float64) h = 1/(TWOPI*I)*(dp/adp/(self.a-cp)) hconj = np.array(map(np.complex.conjugate,h),dtype = np.complex128) g = -sadp*hconj K = np.array([-TWOPI/N*sadp*sadp[t]*1/(TWOPI*I)*(dp/adp/(cp-cp [t])-(dp[t]/adp[t]/(cp-cp[t])).conjugate()) for t in np.arange(N)], dtype = np.complex128) for i in xrange(N): K[i,i] = 1 phi = np.linalg.solve(K,g) temp = -I*phi**2*dp self.theta_array = np.empty(N+1,dtype = np.float64) for i in range(N): self.theta_array[i] = math.atan2(temp[i].imag, temp [i].real) self.theta_array[N] = TWOPI return self.theta_array def get_theta_array(self): if self.theta_array is None: return "Theta array has not yet been generated. Call generate_theta_array()" else: return self.theta_array def get_theta_points(self): if self.theta_array == None: return "Theta array has not yet been generated. Call generate_theta_array()" else: return np.column_stack([np.concatenate([self.tk,np.array ([TWOPI])]),self.theta_array]) def generate_interior_mapper(self): #Generates the necessary code, call riemann_map to find value at a pt cdef unsigned N = self.N cdef coeff = 3*I/(8*N) dp = self.dp theta_array = self.theta_array p_vector = np.zeros(N+1,dtype = np.complex128) for k in range(N//3): p_vector[3*k] = 2*coeff*dp[3*k]*exp(I*theta_array[3*k]) p_vector[3*k+1] = 3*coeff*dp[3*k+1]*exp(I*theta_array[3*k +1]) p_vector[3*k+2] = 3*coeff*dp[3*k+2]*exp(I*theta_array[3*k +2]) p_vector[0] = 1*coeff*dp[0]*exp(I*theta_array[0]) if (N)%3 == 0: p_vector[N] = 1*coeff*dp[0]*exp(I*theta_array[0]) elif (N-2)%3 == 0: p_vector[N-2] = (coeff + I/(3*N))*dp[N-2]*exp(I*theta_array [N-2]) p_vector[N-1] = 4*I/(3*N)*dp[N-1]*exp(I*theta_array[N-1]) p_vector[N] = I/(3*N)*dp[0]*exp(I*theta_array[0]) else: p_vector[N-4] = (coeff + I/(3*N))*dp[N-4]*exp(I*theta_array [N-4]) p_vector[N-3] = 4*I/(3*N)*dp[N-3]*exp(I*theta_array[N-3]) p_vector[N-2] = 2*I/(3*N)*dp[N-2]*exp(I*theta_array[N-2]) p_vector[N-1] = 4*I/(3*N)*dp[N-1]*exp(I*theta_array[N-1]) p_vector[N] = I/(3*N)*dp[0]*exp(I*theta_array[0]) self.p_vector = p_vector def riemann_map(self, pt): if self.p_vector == None: print "please call generate_interior_mapper()" else: cp = self.cp q_vector = np.concatenate([1/(cp-pt),np.array([1/(cp[0]- pt)])]) return -np.dot(self.p_vector,q_vector) I call it in sage like so: alpha = .95 fsym(t) = (alpha^2*cos(2*t)+(1-alpha^4*sin(2*t)^2)^.5)^.5*e^(I*t) fpsym = fsym.derivative() def f(x): return np.complex(fsym(x)) def fprime(x): return np.complex(fpsym(x)) time m = Riemann_Map(f,fprime,-.75,512,init_full = True) --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support-unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~----------~----~----~----~------~----~------~--~---
