Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/master' into pr-64/aboucaud/remo…
Browse files Browse the repository at this point in the history 
…ve-python2-future-imports
  • Loading branch information
@EiffL
EiffL committed Oct 17, 2020
2 parents 60e51f2 + e0c7b37 commit ec580b0
Show file tree
Hide file tree
Showing 4 changed files with 101 additions and 93 deletions.
170 changes: 85 additions & 85 deletions jax_cosmo/scipy/integrate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,33 +21,33 @@

def _difftrap1(function, interval):
"""
Perform part of the trapezoidal rule to integrate a function.
Assume that we had called difftrap with all lower powers-of-2
starting with 1. Calling difftrap only returns the summation
of the new ordinates. It does _not_ multiply by the width
of the trapezoids. This must be performed by the caller.
'function' is the function to evaluate (must accept vector arguments).
'interval' is a sequence with lower and upper limits
of integration.
'numtraps' is the number of trapezoids to use (must be a
power-of-2).
"""
Perform part of the trapezoidal rule to integrate a function.
Assume that we had called difftrap with all lower powers-of-2
starting with 1. Calling difftrap only returns the summation
of the new ordinates. It does _not_ multiply by the width
of the trapezoids. This must be performed by the caller.
'function' is the function to evaluate (must accept vector arguments).
'interval' is a sequence with lower and upper limits
of integration.
'numtraps' is the number of trapezoids to use (must be a
power-of-2).
"""
return 0.5 * (function(interval[0]) + function(interval[1]))


def _difftrapn(function, interval, numtraps):
"""
Perform part of the trapezoidal rule to integrate a function.
Assume that we had called difftrap with all lower powers-of-2
starting with 1. Calling difftrap only returns the summation
of the new ordinates. It does _not_ multiply by the width
of the trapezoids. This must be performed by the caller.
'function' is the function to evaluate (must accept vector arguments).
'interval' is a sequence with lower and upper limits
of integration.
'numtraps' is the number of trapezoids to use (must be a
power-of-2).
"""
Perform part of the trapezoidal rule to integrate a function.
Assume that we had called difftrap with all lower powers-of-2
starting with 1. Calling difftrap only returns the summation
of the new ordinates. It does _not_ multiply by the width
of the trapezoids. This must be performed by the caller.
'function' is the function to evaluate (must accept vector arguments).
'interval' is a sequence with lower and upper limits
of integration.
'numtraps' is the number of trapezoids to use (must be a
power-of-2).
"""
numtosum = numtraps // 2
h = (1.0 * interval[1] - 1.0 * interval[0]) / numtosum
lox = interval[0] + 0.5 * h
Expand All @@ -58,75 +58,75 @@ def _difftrapn(function, interval, numtraps):

def _romberg_diff(b, c, k):
"""
Compute the differences for the Romberg quadrature corrections.
See Forman Acton's "Real Computing Made Real," p 143.
"""
Compute the differences for the Romberg quadrature corrections.
See Forman Acton's "Real Computing Made Real," p 143.
"""
tmp = 4.0 ** k
return (tmp * c - b) / (tmp - 1.0)


def romb(function, a, b, args=(), divmax=6, return_error=False):
"""
Romberg integration of a callable function or method.
Returns the integral of `function` (a function of one variable)
over the interval (`a`, `b`).
If `show` is 1, the triangular array of the intermediate results
will be printed. If `vec_func` is True (default is False), then
`function` is assumed to support vector arguments.
Parameters
----------
function : callable
Function to be integrated.
a : float
Lower limit of integration.
b : float
Upper limit of integration.
Returns
-------
results : float
Result of the integration.
Other Parameters
----------------
args : tuple, optional
Extra arguments to pass to function. Each element of `args` will
be passed as a single argument to `func`. Default is to pass no
extra arguments.
divmax : int, optional
Maximum order of extrapolation. Default is 10.
See Also
--------
fixed_quad : Fixed-order Gaussian quadrature.
quad : Adaptive quadrature using QUADPACK.
dblquad : Double integrals.
tplquad : Triple integrals.
romb : Integrators for sampled data.
simps : Integrators for sampled data.
cumtrapz : Cumulative integration for sampled data.
ode : ODE integrator.
odeint : ODE integrator.
References
----------
.. [1] 'Romberg's method' http://en.wikipedia.org/wiki/Romberg%27s_method
Examples
--------
Integrate a gaussian from 0 to 1 and compare to the error function.
>>> from scipy import integrate
>>> from scipy.special import erf
>>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
>>> result = integrate.romberg(gaussian, 0, 1, show=True)
Romberg integration of <function vfunc at ...> from [0, 1]
::
Steps StepSize Results
1 1.000000 0.385872
2 0.500000 0.412631 0.421551
4 0.250000 0.419184 0.421368 0.421356
8 0.125000 0.420810 0.421352 0.421350 0.421350
16 0.062500 0.421215 0.421350 0.421350 0.421350 0.421350
32 0.031250 0.421317 0.421350 0.421350 0.421350 0.421350 0.421350
The final result is 0.421350396475 after 33 function evaluations.
>>> print("%g %g" % (2*result, erf(1)))
0.842701 0.842701
"""
Romberg integration of a callable function or method.
Returns the integral of `function` (a function of one variable)
over the interval (`a`, `b`).
If `show` is 1, the triangular array of the intermediate results
will be printed. If `vec_func` is True (default is False), then
`function` is assumed to support vector arguments.
Parameters
----------
function : callable
Function to be integrated.
a : float
Lower limit of integration.
b : float
Upper limit of integration.
Returns
-------
results : float
Result of the integration.
Other Parameters
----------------
args : tuple, optional
Extra arguments to pass to function. Each element of `args` will
be passed as a single argument to `func`. Default is to pass no
extra arguments.
divmax : int, optional
Maximum order of extrapolation. Default is 10.
See Also
--------
fixed_quad : Fixed-order Gaussian quadrature.
quad : Adaptive quadrature using QUADPACK.
dblquad : Double integrals.
tplquad : Triple integrals.
romb : Integrators for sampled data.
simps : Integrators for sampled data.
cumtrapz : Cumulative integration for sampled data.
ode : ODE integrator.
odeint : ODE integrator.
References
----------
.. [1] 'Romberg's method' http://en.wikipedia.org/wiki/Romberg%27s_method
Examples
--------
Integrate a gaussian from 0 to 1 and compare to the error function.
>>> from scipy import integrate
>>> from scipy.special import erf
>>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
>>> result = integrate.romberg(gaussian, 0, 1, show=True)
Romberg integration of <function vfunc at ...> from [0, 1]
::
Steps StepSize Results
1 1.000000 0.385872
2 0.500000 0.412631 0.421551
4 0.250000 0.419184 0.421368 0.421356
8 0.125000 0.420810 0.421352 0.421350 0.421350
16 0.062500 0.421215 0.421350 0.421350 0.421350 0.421350
32 0.031250 0.421317 0.421350 0.421350 0.421350 0.421350 0.421350
The final result is 0.421350396475 after 33 function evaluations.
>>> print("%g %g" % (2*result, erf(1)))
0.842701 0.842701
"""
vfunc = jit(lambda x: function(x, *args))

n = 1
Expand Down
12 changes: 6 additions & 6 deletions jax_cosmo/scipy/interpolate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,15 @@
@functools.partial(vmap, in_axes=(0, None, None))
def interp(x, xp, fp):
"""
Simple equivalent of np.interp that compute a linear interpolation.
Simple equivalent of np.interp that compute a linear interpolation.
We are not doing any checks, so make sure your query points are lying
inside the array.
We are not doing any checks, so make sure your query points are lying
inside the array.
TODO: Implement proper interpolation!
TODO: Implement proper interpolation!
x, xp, fp need to be 1d arrays
"""
x, xp, fp need to be 1d arrays
"""
# First we find the nearest neighbour
ind = np.argmin((x - xp) ** 2)

Expand Down
9 changes: 8 additions & 1 deletion jax_cosmo/sparse.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -337,6 +337,10 @@ def slogdet(sparse):
Based on equation (2.2) of https://arxiv.org/abs/1112.4379
For a zero sparse matrix, the result of this computation is
currently undefined and will return nan.
TODO: support null matrix as input.
Parameters
----------
sparse : array
Expand All @@ -346,7 +350,7 @@ def slogdet(sparse):
-------
tuple
Tuple (sign, logdet) such that sign * exp(logdet) is the
determinant. If the determinant is zero, logdet = -inf.
determinant.
"""
sparse = check_sparse(sparse, square=True)
N, _, P = sparse.shape
Expand All @@ -368,6 +372,9 @@ def det(sparse):
Uses :func:`slogdet`.
For a zero sparse matrix, the result of this computation is
currently undefined and will return nan.
Parameters
----------
sparse : array
Expand Down
3 changes: 2 additions & 1 deletion tests/test_sparse.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,5 +81,6 @@ def test_det():
]
)
assert_array_equal(-det(-X), det(X))
assert_array_equal(det(0.0 * X), 0.0)
# TODO: Add proper support for 0 matrix
# assert_array_equal(det(0.0 * X), 0.0)
assert_allclose(det(X), np.linalg.det(to_dense(X)), rtol=1e-6)

0 comments on commit ec580b0

Please sign in to comment.