From 18a50edb23691c62e4af1f89f9fd068f91144cce Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Thu, 4 Dec 2014 01:48:54 +0000 Subject: [PATCH] Simplified interpolations Made dea3 more robust --- pywafo/src/wafo/integrate.py | 2 +- pywafo/src/wafo/interpolate.py | 70 +++++++++++++++------ pywafo/src/wafo/misc.py | 2 +- pywafo/src/wafo/polynomial.py | 2 +- pywafo/src/wafo/stats/_continuous_distns.py | 4 +- 5 files changed, 57 insertions(+), 23 deletions(-) diff --git a/pywafo/src/wafo/integrate.py b/pywafo/src/wafo/integrate.py index 5df54a2..1e5b090 100644 --- a/pywafo/src/wafo/integrate.py +++ b/pywafo/src/wafo/integrate.py @@ -81,7 +81,7 @@ def dea3(v0, v1, v2): smalle2 = (abs(ss * E1) <= 1.0e-3).ravel() result = 1.0 * E2 - abserr = err1 + err2 + E2 * _EPS * 10.0 + abserr = err1 + err2 + abs(E2) * _EPS * 10.0 converged = (err1 <= tol1) & (err2 <= tol2).ravel() | smalle2 k4, = (1 - converged).nonzero() if k4.size > 0: diff --git a/pywafo/src/wafo/interpolate.py b/pywafo/src/wafo/interpolate.py index 030ae11..37947be 100644 --- a/pywafo/src/wafo/interpolate.py +++ b/pywafo/src/wafo/interpolate.py @@ -4,11 +4,8 @@ import numpy as np import scipy.signal # import scipy.sparse.linalg # @UnusedImport import scipy.sparse as sparse -from numpy.ma.core import ones, zeros, prod, sin -from numpy import diff, pi, inf # @UnresolvedImport -from numpy.lib.shape_base import vstack -from numpy.lib.function_base import linspace -from scipy.interpolate import PiecewisePolynomial +from numpy import ones, zeros, prod, sin, diff, pi, inf, vstack, linspace +from scipy.interpolate import PiecewisePolynomial, interp1d import polynomial as pl @@ -323,17 +320,11 @@ class PPform(object): indxs = indxs.clip(0, len(self.breaks)) pp = self.coeffs dx = xx - self.breaks.take(indxs) - if True: - v = pp[0, indxs] - for i in xrange(1, self.order): - v = dx * v + pp[i, indxs] - values = v - else: - V = np.vander(dx, N=self.order) - # values = np.diag(dot(V,pp[:,indxs])) - dot = np.dot - values = np.array([dot(V[k, :], pp[:, indxs[k]]) - for k in xrange(len(xx))]) + + v = pp[0, indxs] + for i in xrange(1, self.order): + v = dx * v + pp[i, indxs] + values = v res[mask] = values res.shape = saveshape @@ -1060,6 +1051,46 @@ class Pchip(PiecewisePolynomial): super(Pchip, self).__init__(x, zip(y, yp), orders=3) +def interp3(x, y, z, v, xi, yi, zi, method='cubic'): + """Interpolation on 3-D. x, y, xi, yi should be 1-D + and z.shape == (len(x), len(y), len(z))""" + q = (x, y, z) + qi = (xi, yi, zi) + for j in range(3): + pp = interp1d(q[j], v, axis=j, kind=method) + v = pp(qi[j]) + return v + + +def somefunc(x, y, z): + return x**2 + y**2 - z**2 + x*y*z + + +def test_interp3(): + # some input data + x = np.linspace(0, 1, 5) + y = np.linspace(0, 2, 6) + z = np.linspace(0, 3, 7) + v = somefunc(x[:, None, None], y[None, :, None], z[None, None, :]) + + # interpolate + xi = np.linspace(0, 1, 45) + yi = np.linspace(0, 2, 46) + zi = np.linspace(0, 3, 47) + vi = interp3(x, y, z, v, xi, yi, zi) + + import matplotlib.pyplot as plt + X, Y = np.meshgrid(xi, yi) + plt.figure(1) + plt.subplot(1, 2, 1) + plt.pcolor(X, Y, vi[:, :, 12].T) + plt.title('interpolated') + plt.subplot(1, 2, 2) + plt.pcolor(X, Y, somefunc(xi[:, None], yi[None, :], zi[12]).T) + plt.title('exact') + plt.show('hold') + + def test_smoothing_spline(): x = linspace(0, 2 * pi + pi / 4, 20) y = sin(x) # + np.random.randn(x.size) @@ -1074,7 +1105,7 @@ def test_smoothing_spline(): import matplotlib.pyplot as plb plb.plot(x, y, x1, y1, '.', x1, dy1, 'ro', x1, y01, 'r-') - plb.show() + plb.show('hold') pass # tck = interpolate.splrep(x, y, s=len(x)) @@ -1234,6 +1265,7 @@ def test_docstrings(): if __name__ == '__main__': # test_func() # test_doctstrings() - # test_smoothing_spline() + test_smoothing_spline() # compare_methods() - demo_monoticity() + # demo_monoticity() + # test_interp3() diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index 77ef33a..6cfc62b 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -2130,7 +2130,7 @@ def dea3(v0, v1, v2): smallE2 = (abs(ss * E1) <= 1.0e-3).ravel() result = 1.0 * E2 - abserr = err1 + err2 + E2 * _EPS * 10.0 + abserr = err1 + err2 + abs(E2) * _EPS * 10.0 converged = (err1 <= tol1) & (err2 <= tol2).ravel() | smallE2 k4, = (1 - converged).nonzero() if k4.size > 0: diff --git a/pywafo/src/wafo/polynomial.py b/pywafo/src/wafo/polynomial.py index 7294157..4f199c2 100644 --- a/pywafo/src/wafo/polynomial.py +++ b/pywafo/src/wafo/polynomial.py @@ -22,7 +22,7 @@ from plotbackend import plotbackend as plt import numpy as np from numpy.fft import fft, ifft from numpy import (zeros, ones, zeros_like, array, asarray, newaxis, arange, - logical_or, any, pi, cos, round, diff, all, exp, + logical_or, any, pi, cos, round, diff, all, exp, atleast_1d, where, extract, linalg, sign, concatenate, floor, isreal, conj, remainder, linspace, sum) from numpy.lib.polynomial import * # @UnusedWildImport diff --git a/pywafo/src/wafo/stats/_continuous_distns.py b/pywafo/src/wafo/stats/_continuous_distns.py index 38ca3f5..7c3bc79 100644 --- a/pywafo/src/wafo/stats/_continuous_distns.py +++ b/pywafo/src/wafo/stats/_continuous_distns.py @@ -11,7 +11,9 @@ from scipy.misc.doccer import inherit_docstring_from from scipy import special from scipy import optimize from scipy import integrate -from scipy.special import (gammaln as gamln, gamma as gam, boxcox, boxcox1p, log1p, expm1) +from scipy.special import (gammaln as gamln, gamma as gam, boxcox, boxcox1p, + log1p, expm1) + #inv_boxcox, inv_boxcox1p) from numpy import (where, arange, putmask, ravel, sum, shape, log, sqrt, exp, arctanh, tan, sin, arcsin, arctan,