Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

cython-lint for calculus/ folder #39636

Open
wants to merge 2 commits into
base: develop
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 6 additions & 5 deletions src/sage/calculus/integration.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ cdef double c_f(double t, void *params) noexcept:

def numerical_integral(func, a, b=None,
algorithm='qag',
max_points=87, params=[], eps_abs=1e-6,
max_points=87, params=None, eps_abs=1e-6,
eps_rel=1e-6, rule=6):
r"""
Return the numerical integral of the function on the interval
Expand Down Expand Up @@ -267,8 +267,6 @@ def numerical_integral(func, a, b=None,
"""
cdef double abs_err # step size
cdef double result
cdef int i
cdef int j
cdef double _a, _b
cdef PyFunctionWrapper wrapper # struct to pass information into GSL C function

Expand All @@ -288,6 +286,9 @@ def numerical_integral(func, a, b=None,
cdef gsl_integration_workspace* W
W = NULL

if params is None:
params = []

if True:
from sage.rings.infinity import Infinity
try:
Expand Down Expand Up @@ -619,7 +620,7 @@ def monte_carlo_integral(func, xl, xu, size_t calls, algorithm='plain',
if len(vars) < target_dim:
raise ValueError(("The function to be integrated depends on "
"{} variables {}, and so cannot be "
"integrated in {} dimensions. Please fix "
"integrated in {} dimensions. Please fix "
"additional variables with the 'params' "
"argument").format(len(vars), tuple(vars),
target_dim))
Expand All @@ -628,7 +629,7 @@ def monte_carlo_integral(func, xl, xu, size_t calls, algorithm='plain',
"{} variables {}, and so cannot be "
"integrated in {} dimensions. Please add "
"more items in upper and lower limits"
).format(len(vars), tuple(vars), target_dim))
).format(len(vars), tuple(vars), target_dim))

from sage.structure.element import Expression
if isinstance(func, Expression):
Expand Down
4 changes: 2 additions & 2 deletions src/sage/calculus/ode.pxd
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
cdef class ode_system:
cdef int c_j(self, double , double *, double *, double *) noexcept
cdef int c_j(self, double , double *, double *, double *) noexcept

cdef int c_f(self, double t, double* , double* ) noexcept
cdef int c_f(self, double t, double* , double*) noexcept
10 changes: 6 additions & 4 deletions src/sage/calculus/ode.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -342,7 +342,9 @@ class ode_solver():
sage: with NamedTemporaryFile(suffix='.png') as f:
....: T.plot_solution(i=0, filename=f.name)
"""
def __init__(self, function=None, jacobian=None, h=1e-2, error_abs=1e-10, error_rel=1e-10, a=False, a_dydt=False, scale_abs=False, algorithm='rkf45', y_0=None, t_span=None, params=[]):
def __init__(self, function=None, jacobian=None, h=1e-2, error_abs=1e-10,
error_rel=1e-10, a=False, a_dydt=False, scale_abs=False,
algorithm='rkf45', y_0=None, t_span=None, params=None):
self.function = function
self.jacobian = jacobian
self.h = h
Expand All @@ -354,7 +356,7 @@ class ode_solver():
self.algorithm = algorithm
self.y_0 = y_0
self.t_span = t_span
self.params = params
self.params = [] if params is None else params
self.solution = []

def __setattr__(self, name, value):
Expand Down Expand Up @@ -407,15 +409,15 @@ class ode_solver():
else:
G.save(filename=filename)

def ode_solve(self, t_span=False, y_0=False, num_points=False, params=[]):
def ode_solve(self, t_span=False, y_0=False, num_points=False, params=None):
cdef double h # step size
h = self.h
cdef int i
cdef int j
cdef int type
cdef int dim
cdef PyFunctionWrapper wrapper # struct to pass information into GSL C function
self.params = params
self.params = [] if params is None else params

if t_span:
self.t_span = t_span
Expand Down
116 changes: 68 additions & 48 deletions src/sage/calculus/riemann.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ from math import pi
from math import sin
from math import cos

from math import log # used for complex plot lightness
from math import log # used for complex plot lightness
from math import atan

from cmath import exp
Expand Down Expand Up @@ -201,7 +201,7 @@ cdef class Riemann_Map:
cdef exterior

def __init__(self, fs, fprimes, COMPLEX_T a, int N=500, int ncorners=4,
opp=False, exterior = False):
opp=False, exterior=False):
"""
Initialize the ``Riemann_Map`` class. See the class :class:`Riemann_Map`
for full documentation on the input of this initialization method.
Expand All @@ -227,7 +227,7 @@ cdef class Riemann_Map:
self.f = fs[0]
self.a = a
self.ncorners = ncorners
self.N = N # Number of collocation pts
self.N = N # Number of collocation pts
self.opp = opp
self.exterior = exterior
self.tk = np.array(np.arange(N) * TWOPI / N + 0.001 / N,
Expand All @@ -236,14 +236,14 @@ cdef class Riemann_Map:
for i in range(N):
self.tk2[i] = self.tk[i]
self.tk2[N] = TWOPI
self.B = len(fs) # number of boundaries of the figure
self.B = len(fs) # number of boundaries of the figure
if self.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)
dtype=COMPLEX)
cdef np.ndarray[COMPLEX_T,ndim=2] dps = np.zeros([self.B, N],
dtype=COMPLEX)
dtype=COMPLEX)
# Find the points on the boundaries and their derivatives.
if self.exterior:
for k in range(self.B):
Expand Down Expand Up @@ -321,22 +321,21 @@ cdef class Riemann_Map:
hconj = h.conjugate()
g = -sadp * hconj
normalized_dp=dp/adp
C = I / N * sadp # equivalent to -TWOPI / N * 1 / (TWOPI * I) * sadp
errinvalid = np.geterr()['invalid'] # checks the current error handling for invalid
errdivide = np.geterr()['divide'] # checks the current error handling for divide
C = I / N * sadp # equivalent to -TWOPI / N * 1 / (TWOPI * I) * sadp
errinvalid = np.geterr()['invalid'] # checks the current error handling for invalid
errdivide = np.geterr()['divide'] # checks the current error handling for divide
np.seterr(divide='ignore',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=np.complex128)
np.seterr(divide=errdivide,invalid=errinvalid) # resets the error handling
(normalized_dp[t]/(cp-cp[t])).conjugate())
for t in np.arange(NB)], dtype=np.complex128)
np.seterr(divide=errdivide,invalid=errinvalid) # resets the error handling
for i in range(NB):
K[i, i] = 1
# Nystrom Method for solving 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])
start = 0
# Finding the theta correspondence using phase. Misbehaves for some
# regions.
if B != 1:
Expand Down Expand Up @@ -431,7 +430,6 @@ cdef class Riemann_Map:
sage: sz0 = m.get_szego(boundary=0)
sage: sz1 = m.get_szego(boundary=1)
"""
cdef int k, B
if boundary < 0:
temptk = self.tk
for i in range(self.B - 1):
Expand Down Expand Up @@ -706,10 +704,11 @@ cdef class Riemann_Map:
else:
return mapped

def plot_boundaries(self, plotjoined=True, rgbcolor=[0,0,0], thickness=1):
def plot_boundaries(self, plotjoined=True, rgbcolor=None, thickness=1):
"""
Plots the boundaries of the region for the Riemann map. Note that
this method DOES work for multiply connected domains.
Plot the boundaries of the region for the Riemann map.

Note that this method DOES work for multiply connected domains.

INPUT:

Expand Down Expand Up @@ -747,6 +746,9 @@ cdef class Riemann_Map:
"""
from sage.plot.all import list_plot

if rgbcolor is None:
rgbcolor = [0, 0, 0]

plots = list(range(self.B))
for k in range(self.B):
# This conditional should be eliminated when the thickness/pointsize
Expand Down Expand Up @@ -827,8 +829,9 @@ cdef class Riemann_Map:

@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, min_mag=0.001, **options):
rgbcolor=None, thickness=1,
plotjoined=True, withcolor=False,
plot_points=200, min_mag=0.001, **options):
"""
Generate a traditional "spiderweb plot" of the Riemann map.

Expand Down Expand Up @@ -944,12 +947,16 @@ cdef class Riemann_Map:
cdef int k, i
if self.exterior:
raise ValueError(
"Spiderwebs for exterior maps are not currently supported")
"Spiderwebs for exterior maps are not currently supported")

if rgbcolor is None:
rgbcolor = [0, 0, 0]

if self.B == 1: # The efficient simply connected
edge = self.plot_boundaries(plotjoined=plotjoined,
rgbcolor=rgbcolor, thickness=thickness)
rgbcolor=rgbcolor,
thickness=thickness)
circle_list = 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]
Expand All @@ -960,10 +967,13 @@ cdef class 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)
rgbcolor=rgbcolor,
thickness=thickness,
plotjoined=True)
else:
circle_list[k] = list_plot(comp_pt(temp, 1),
rgbcolor=rgbcolor, pointsize=thickness)
rgbcolor=rgbcolor,
pointsize=thickness)
line_list = list(range(spokes))
for k in range(spokes):
temp = list(range(pts))
Expand All @@ -989,24 +999,31 @@ cdef class Riemann_Map:
self.plot_colored(plot_points=plot_points)
else:
return edge + sum(circle_list) + sum(line_list)
else: # The more difficult multiply connected
else: # The more difficult multiply connected
z_values, xmin, xmax, ymin, ymax = self.compute_on_grid([],
plot_points)
plot_points)
xstep = (xmax-xmin)/plot_points
ystep = (ymax-ymin)/plot_points
dr, dtheta= get_derivatives(z_values, xstep, ystep) # clean later
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, min_mag),
(xmin, xmax), (ymin, ymax),options))
g.add_primitive(ComplexPlot(complex_to_spiderweb(z_values, dr,
dtheta, spokes,
circles,
rgbcolor,
thickness,
withcolor,
min_mag),
(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):
def plot_colored(self, plot_range=None, int plot_points=100, **options):
"""
Generates a colored plot of the Riemann map. A red point on the
colored plot corresponds to a red point on the unit disc.
Generate a colored plot of the Riemann map.

A red point on the colored plot corresponds to a red point on
the unit disc.

INPUT:

Expand Down Expand Up @@ -1053,11 +1070,14 @@ cdef class Riemann_Map:
from sage.plot.complex_plot import ComplexPlot
from sage.plot.all import Graphics

if plot_range is None:
plot_range = []

z_values, xmin, xmax, ymin, ymax = self.compute_on_grid(plot_range,
plot_points)
plot_points)
g = Graphics()
g.add_primitive(ComplexPlot(complex_to_rgb(z_values), (xmin, xmax),
(ymin, ymax),options))
(ymin, ymax), options))
return g

cdef comp_pt(clist, loop=True):
Expand Down Expand Up @@ -1089,8 +1109,8 @@ cdef comp_pt(clist, loop=True):
list2.append(list2[0])
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):
"""
Compute the r*e^(I*theta) form of derivatives from the grid of points. The
derivatives are computed using quick-and-dirty taylor expansion and
Expand Down Expand Up @@ -1132,21 +1152,21 @@ cpdef get_derivatives(np.ndarray[COMPLEX_T, ndim=2] z_values, FLOAT_T xstep,
"""
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
# (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 the magnitude of its
# 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)
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, min_mag):
np.ndarray[FLOAT_T, ndim = 2] dr,
np.ndarray[FLOAT_T, ndim = 2] dtheta,
spokes, circles, rgbcolor, thickness,
withcolor, min_mag):
"""
Convert a grid of complex numbers into a matrix containing rgb data
for the Riemann spiderweb plot.
Expand Down Expand Up @@ -1221,9 +1241,9 @@ cpdef complex_to_spiderweb(np.ndarray[COMPLEX_T, ndim = 2] z_values,
[1. , 1. , 1. ]]])
"""
cdef Py_ssize_t i, j, imax, jmax
cdef FLOAT_T x, y, mag, arg, width, target, precision, dmag, darg
cdef FLOAT_T mag, arg, target, precision, dmag, darg
cdef COMPLEX_T z
cdef FLOAT_T DMAX = 70 # change to adjust rate_of_change cutoff below
cdef FLOAT_T DMAX = 70 # change to adjust rate_of_change cutoff below
precision = thickness/150.0
imax = len(z_values)
jmax = len(z_values[0])
Expand All @@ -1234,7 +1254,7 @@ cpdef complex_to_spiderweb(np.ndarray[COMPLEX_T, ndim = 2] z_values,
rgb = np.zeros(dtype=FLOAT, shape=(imax, jmax, 3))
rgb += 1
if circles != 0:
circ_radii = srange(0,1.0,1.0/circles)
circ_radii = srange(0, 1.0, 1.0/circles)
else:
circ_radii = []
if spokes != 0:
Expand Down Expand Up @@ -1301,7 +1321,7 @@ cpdef complex_to_rgb(np.ndarray[COMPLEX_T, ndim=2] z_values):
TypeError: Argument 'z_values' has incorrect type (expected numpy.ndarray, got list)
"""
cdef Py_ssize_t i, j, imax, jmax
cdef FLOAT_T x, y, mag, arg
cdef FLOAT_T mag, arg
cdef FLOAT_T lightness, hue, top, bot
cdef FLOAT_T r, g, b
cdef int ihue
Expand Down Expand Up @@ -1494,7 +1514,7 @@ cpdef analytic_interior(COMPLEX_T z, int n, FLOAT_T epsilon):
# 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 = [epsilon,z,n,'i'])[0]
params = [epsilon,z,n,'i'])[0]
ip = 1/(TWOPI*I)*numerical_integral(cauchy_kernel,0,2*pi,
params = [epsilon,z,n,'r'])[0]
params = [epsilon,z,n,'r'])[0]
return rp + ip
Loading