From 0a23402fc65b4e20152ceaa7b1aab9dfccb8bc28 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Thu, 9 Jul 2015 16:03:56 +0000 Subject: [PATCH] Added padua.py --- pywafo/src/wafo/padua.py | 486 +++++++++++++++++++++++++++++ pywafo/src/wafo/test/test_padua.py | 89 ++++++ 2 files changed, 575 insertions(+) create mode 100644 pywafo/src/wafo/padua.py create mode 100644 pywafo/src/wafo/test/test_padua.py diff --git a/pywafo/src/wafo/padua.py b/pywafo/src/wafo/padua.py new file mode 100644 index 0000000..4aa31a2 --- /dev/null +++ b/pywafo/src/wafo/padua.py @@ -0,0 +1,486 @@ +''' + +All the software contained in this library is protected by copyright. +Permission to use, copy, modify, and distribute this software for any +purpose without fee is hereby granted, provided that this entire notice +is included in all copies of any software which is or includes a copy +or modification of this software and in all copies of the supporting +documentation for such software. + +THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED +WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY +MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", +NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY +MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF +USING THE SOFTWARE LIES WITH THE PARTY DOING SO. + +ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE +ABOVE STATEMENT. + + +AUTHORS: +Per A Brodtkorb +Python code Based on matlab code written by: + +Marco Caliari +University of Verona, Italy +E-mail: marco.caliari@univr.it + +Stefano de Marchi, Alvise Sommariva, Marco Vianello +University of Padua, Italy +E-mail: demarchi@math.unipd.it, alvise@math.unipd.it, + marcov@math.unipd.it + +Reference +--------- +Padua2DM: fast interpolation and cubature at the Padua points in Matlab/Octave +NUMERICAL ALGORITHMS, 56 (2011), PP. 45-60 + + +Padua module +------------ +In polynomial interpolation of two variables, the Padua points are the first +known example (and up to now the only one) of a unisolvent point set +(that is, the interpolating polynomial is unique) with minimal growth of their +Lebesgue constant, proven to be O(log2 n). +This module provides all the functions needed to perform interpolation and +cubature at the Padua points, together with the functions and the demos used +in the paper. + +pdint.m : main function for interpolation at the Padua points +pdcub.m : main function for cubature at the Padua points +pdpts.m : function for the computation of the Padua points +padua_fit.m : function for the computation of the interpolation + coefficients by FFT (recommended) +pdcfsMM.m : function for the computation of the interpolation + coefficients by matrix multiplications +pdval.m : function for the evaluation of the interpolation + polynomial +pdwtsFFT.m : function for the computation of the cubature + weights by FFT +pdwtsMM.m : function for the computation of the cubature + weights by matrix multiplications (recommended) +funct.m : function containing some test functions +demo_pdint.m : demo script for pdint +demo_cputime_pdint.m : demo script for the computation of CPU time for + interpolation +demo_errors_pdint.m : demo script for the comparison of interpolation with + coefficients computed by FFT or by matrix + multiplications +demo_pdcub : demo script for pdcub +demo_cputime_pdcub.m : demo script for the computation of CPU time for + cubature +demo_errors_pdcub.m : demo script for the comparison of cubature with + weights computed by FFT or by matrix multiplications +demo_errors_pdcub_gl.m : demo script for the comparison of different cubature + formulas +cubature_square.m : function for the computation of some cubature + formulas for the square +omelyan_solovyan_rule.m : function for the computation of Omelyan-Solovyan + cubature points and weights +Contents.m : Contents file for Matlab + + +''' +from __future__ import division +import numpy as np +from numpy.fft import fft + +class TestFunctions(object): + ''' + Computes test function in the points (x, y) + + Parameters + ---------- + x,y : array-like + evaluate the function in the points (x,y) + id : scalar int (default 0) + id defining which test function to use. Options are + 0: franke + 1: half_sphere + 2: poly_degree20 + 3: exp_fun1 + 4: exp_fun100 + 5: cos30 + 6: constant + 7: exp_xy + 8: runge + 9: abs_cubed + 10: gauss + 11: exp_inv + + Returns + ------- + z : array-like + value of the function in the points (x,y) + ''' + @staticmethod + def franke(x, y): + '''Franke function. + + The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 500, is + 2.1547794245591083e+000 with an estimated absolute error of 8.88e-016. + + The value of the definite integral on the square [0,1] x [0,1], + obtained using a Padua Points cubature formula of degree 500, is + 4.06969589491556e-01 with an estimated absolute error of 8.88e-016. + + Maple: 0.40696958949155611906 + ''' + exp = np.exp + return (3. / 4 * exp(-((9. * x - 2)**2 + (9. * y - 2)**2) / 4) + + 3. / 4 * exp(-(9. * x + 1)**2 / 49 - (9. * y + 1) / 10) + + 1. / 2 * exp(-((9. * x - 7)**2 + (9. * y - 3)**2) / 4) - + 1. / 5 * exp(-(9. * x - 4)**2 - (9. * y - 7)**2)) + + @staticmethod + def half_sphere(x, y): + '''The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 2000, is + 3.9129044444568244e+000 with an estimated absolute error of 3.22e-010. + ''' + return ((x - 0.5)**2 + (y - 0.5)**2)**(1. / 2) + + @staticmethod + def poly_degree20(x, y): + ''''Bivariate polynomial having moderate degree. + The value of the definite integral on the square [-1,1] x + [-1,1], up to machine precision, is 18157.16017316017 (see ref. 6). + The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 500, + is 1.8157160173160162e+004. + + 2D modification of an example by L.N.Trefethen (see ref. 7), f(x)=x^20. + ''' + return (x + y)**20 + + @staticmethod + def exp_fun1(x, y): + ''' The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 2000, + is 2.1234596326670683e+001 with an estimated absolute error of 7.11e-015. + ''' + return np.exp((x - 0.5)**2 + (y - 0.5)**2) + + @staticmethod + def exp_fun100(x, y): + '''The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 2000, + is 3.1415926535849605e-002 with an estimated absolute error of 3.47e-017. + ''' + return np.exp(-100 * ((x - 0.5)**2 + (y - 0.5)**2)) + + @staticmethod + def cos30(x, y): + ''' The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 500, + is 4.3386955120336568e-003 with an estimated absolute error of 2.95e-017. + ''' + return np.cos(30 * (x + y)) + + @staticmethod + def constant(x, y): + '''Constant. + To test interpolation and cubature at degree 0. + The value of the definite integral on the square [-1,1] x [-1,1] + is 4. + ''' + return np.ones(np.shape(x)) + + @staticmethod + def exp_xy(x, y): + '''The value of the definite integral on the square [-1,1] x [-1,1] + is up to machine precision is 5.524391382167263 (see ref. 6). + 2. The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 500, + is 5.5243913821672628e+000 with an estimated absolute error of 0.00e+000. + 2D modification of an example by L.N.Trefethen (see ref. 7), + f(x)=exp(x). + ''' + return np.exp(x + y) + + @staticmethod + def runge(x, y): + ''' Bivariate Runge function: as 1D complex function is analytic + in a neighborhood of [-1; 1] but not throughout the complex plane. + + The value of the definite integral on the square [-1,1] x [-1,1], + up to machine precision, is 0.597388947274307 (see ref. 6). + The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 500, + is 5.9738894727430725e-001 with an estimated absolute error of 0.00e+000. + + 2D modification of an example by L.N.Trefethen (see ref. 7), + f(x)=1/(1+16*x^2). + ''' + return 1. / (1 + 16 * (x**2 + y**2)) + + @staticmethod + def abs_cubed(x, y): + '''Low regular function. + The value of the definite integral on the square [-1,1] x [-1,1], + up to machine precision, is 2.508723139534059 (see ref. 6). + The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 500, + is 2.5087231395340579e+000 with an estimated absolute error of 0.00e+000. + + 2D modification of an example by L.N.Trefethen (see ref. 7), + f(x)=abs(x)^3. + ''' + return (x**2 + y**2)**(3 / 2) + + @staticmethod + def gauss(x, y): + '''Bivariate gaussian: smooth function. + The value of the definite integral on the square [-1,1] x [-1,1], + up to machine precision, is 2.230985141404135 (see ref. 6). + The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 500, + is 2.2309851414041333e+000 with an estimated absolute error of 2.66e-015. + + 2D modification of an example by L.N.Trefethen (see ref. 7), + f(x)=exp(-x^2). + ''' + return np.exp(-x**2 - y**2) + + @staticmethod + def exp_inv(x, y): + '''Bivariate example stemming from a 1D C-infinity function. + The value of the definite integral on the square [-1,1] x [-1,1], + up to machine precision, is 0.853358758654305 (see ref. 6). + The value of the definite integral on the square [-1,1] x [-1,1], + obtained using a Padua Points cubature formula of degree 2000, + is 8.5335875865430544e-001 with an estimated absolute error of 3.11e-015. + + 2D modification of an example by L.N.Trefethen (see ref. 7), + f(x)=exp(-1/x^2). + ''' + arg_z = (x**2 + y**2) + # Avoid cases in which "arg_z=0", setting only in those instances + # "arg_z=eps". + arg_z = arg_z + (1 - np.abs(np.sign(arg_z))) * 1.e-100 + arg_z = 1. / arg_z + return np.exp(-arg_z) + + def __call__(self, x, y, id=0): + s = self + test_function = [s.franke, s.half_sphere, s.poly_degree20, s.exp_fun1, + s.exp_fun100, s.cos30, s.constant, s.exp_xy, s.runge, + s.abs_cubed, s.gauss, s.exp_inv] + return test_function[id](x, y) +testfunct = TestFunctions() + + +def _find_m(n): + ix = np.r_[1:(n + 1) * (n + 2):2] + if np.mod(n, 2) == 0: + n2 = n // 2 + offset = np.array([[0, 1] * n2 + [0, ]] * (n2 + 1)) + ix = ix - offset.ravel(order='F') + return ix + + +def padua_points(n, domain=(-1, 1, -1, 1)): + ''' Return Padua points + + Parameters + ---------- + n : scalar integer + interpolation degree + domain : vector [a,b,c,d] + defining the rectangle [a,b] x [c,d]. (default domain = (-1,1,-1,1)) + + Returns + ------- + pad : array of shape (2 x (n+1)*(n+2)/2) such that + (pad[0,:], pad[1,: ]) defines the Padua points in the domain + rectangle [a,b] x [c,d]. + or + X1,Y1,X2,Y2 : arrays + Two subgrids X1,Y1 and X2,Y2 defining the Padua points + ------------------------------------------------------------------------------- + ''' + a, b, c, d = domain + t0 = [np.pi] if n == 0 else np.linspace(0, np.pi, n + 1) + t1 = np.linspace(0, np.pi, n + 2) + zn = (a + b + (b - a) * np.cos(t0)) / 2 + zn1 = (c + d + (d - c) * np.cos(t1)) / 2 + + Pad1, Pad2 = np.meshgrid(zn, zn1) + ix = _find_m(n) + return np.vstack((Pad1.ravel(order='F')[ix], + Pad2.ravel(order='F')[ix])) + + +def error_estimate(C0f): + ''' Return interpolation error estimate from Padua coefficients + ''' + n = C0f.shape[1] + C0f2 = np.fliplr(C0f) + errest = sum(np.abs(np.diag(C0f2))) + if (n >= 1): + errest = errest + sum(np.abs(np.diag(C0f2, -1))) + if (n >= 2): + errest = errest + sum(np.abs(np.diag(C0f2, -2))) + return 2 * errest + + +def padua_fit(Pad, fun, *args): + ''' + Computes the Chebyshevs coefficients + + so that f(x, y) can be approximated by: + + f(x, y) = sum cjk*Tjk(x, y) + + Parameters + ---------- + Pad : array-like + Padua points, as computed with padua_points function. + fun : function to be interpolated in the form + fun(x, y, *args), where *args are optional arguments for fun. + + Returns + ------- + coefficents: coefficient matrix + abs_err : interpolation error estimate + + ''' + + + N = np.shape(Pad)[1] + # recover the degree n from N = (n+1)(n+2)/2 + n = int(round(-3 + np.sqrt(1 + 8 * N)) / 2) + C0f = fun(Pad[0], Pad[1], *args) + if (n > 0): + ix = _find_m(n) + GfT = np.zeros((n + 2) * (n + 1)) + GfT[ix] = C0f * 2 / (n * (n + 1)) + GfT = GfT.reshape(n + 1, n + 2) + GfT = GfT.T + GfT[0] = GfT[0] / 2 + GfT[n + 1] = GfT[n + 1] / 2 + GfT[:, 0] = GfT[:, 0] / 2 + GfT[:, n] = GfT[:, n] / 2 + Gf = GfT.T + # compute the interpolation coefficient matrix C0f by FFT + Gfhat = np.real(fft(Gf, 2 * n, axis=0)) + Gfhathat = np.real(fft(Gfhat[:n + 1, :], 2 * (n + 1), axis=1)) + C0f = 2 * Gfhathat[:, 0:n + 1] + C0f[0] = C0f[0, :] / np.sqrt(2) + C0f[:, 0] = C0f[:, 0] / np.sqrt(2) + C0f = np.fliplr(np.triu(np.fliplr(C0f))) + C0f[n] = C0f[n] / 2 + + return C0f, error_estimate(C0f) + +def paduavals2coefs(f): + useFFTwhenNisMoreThan=100 + m = len(f) + n = int(round(-1.5 + np.sqrt(.25+2*m))) + x = padua_points(n) + idx = _find_m(n) + w = 0*x[0] + 1./(n*(n+1)) + idx1 = np.all(np.abs(x) == 1, axis=0) + w[idx1] = .5*w[idx1] + idx2 = np.all(np.abs(x) != 1, axis=0) + w[idx2] = 2*w[idx2] + + G = np.zeros(idx.max()+1) + G[idx] = 4*w*f + + if ( n < useFFTwhenNisMoreThan ): + t1 = np.r_[0:n+1].reshape(-1, 1) + Tn1 = np.cos(t1*t1.T*np.pi/n) + t2 = np.r_[0:n+2].reshape(-1, 1) + Tn2 = np.cos(t2*t2.T*np.pi/(n+1)); + C = np.dot(Tn2, np.dot(G,Tn1)) + else: + + # dct = @(c) chebtech2.coeffs2vals(c); + C = np.rot90(dct(dct(G.T).T), axis=1) + + C[0] = .5*C[0] + C[:,1] = .5*C[:, 1] + C[0,-1] = .5*C[0, -1] + del C[-1] + + # Take upper-left triangular part: + return np.fliplr(np.triu(np.fliplr(C))) + #C = triu(C(:,end:-1:1)); + #C = C(:,end:-1:1); + + +def padua_fit2(Pad, fun, *args): + N = np.shape(Pad)[1] + # recover the degree n from N = (n+1)(n+2)/2 + n = int(round(-3 + np.sqrt(1 + 8 * N)) / 2) + C0f = fun(Pad[0], Pad[1], *args) + return paduavals2coefs(C0f) + + + + + +def _compute_moments(n): + k = np.r_[0:n:2] + mom = 2 * np.sqrt(2) / (1 - k ** 2) + mom[0] = 2 + return mom + + +def padua_cubature(coefficients, domain=(-1, 1, -1, 1)): + ''' + Compute the integral through the coefficient matrix. + ''' + n = coefficients.shape[1] + mom = _compute_moments(n) + M1, M2 = np.meshgrid(mom, mom) + M = M1 * M2 + C0fM = coefficients[0:n:2, 0:n:2] * M + a, b, c, d = domain + integral = (b - a) * (d - c) * C0fM.sum() / 4 + return integral + + +def padua_val(X, Y, coefficients, domain=(-1, 1, -1, 1), use_meshgrid=False): + ''' + Evaluate polynomial in padua form at X, Y. + + Evaluate the interpolation polynomial defined through its coefficient + matrix coefficients at the target points X(:,1),X(:,2) or at the meshgrid(X1,X2) + + Parameters + ---------- + X, Y: array-like + evaluation points. + coefficients : array-like + coefficient matrix + domain : a vector [a,b,c,d] + defining the rectangle [a,b] x [c,d] + use_meshgrid: bool + If True interpolate at the points meshgrid(X, Y) + + Returns + ------- + fxy : array-like + evaluation of the interpolation polynomial at the target points + ''' + X, Y = np.atleast_1d(X, Y) + original_shape = X.shape + min, max = np.minimum, np.maximum + a, b, c, d = domain + n = np.shape(coefficients)[1] + + X1 = min(max(2 * (X.ravel() - a) / (b - a) - 1, -1), 1).reshape(1, -1) + X2 = min(max(2 * (Y.ravel() - c) / (d - c) - 1, -1), 1).reshape(1, -1) + tn = np.r_[0:n][:, None] + TX1 = np.cos(tn * np.arccos(X1)) + TX2 = np.cos(tn * np.arccos(X2)) + TX1[1:n + 1] = TX1[1:n + 1] * np.sqrt(2) + TX2[1:n + 1] = TX2[1:n + 1] * np.sqrt(2) + if use_meshgrid: + return np.dot(TX1.T, np.dot(coefficients, TX2)).T # eval on meshgrid points + val = np.sum(np.dot(TX1.T, coefficients) * TX2.T, axis=1) # scattered points + return val.reshape(original_shape) diff --git a/pywafo/src/wafo/test/test_padua.py b/pywafo/src/wafo/test/test_padua.py new file mode 100644 index 0000000..c9372fb --- /dev/null +++ b/pywafo/src/wafo/test/test_padua.py @@ -0,0 +1,89 @@ + + +import unittest +import numpy as np +from numpy import cos, pi +from numpy.testing import assert_array_almost_equal +from wafo.padua import (padua_points, testfunct, padua_fit, + padua_fit2, + padua_cubature, padua_val) + + +class PaduaTestCase(unittest.TestCase): + def test_padua_points_degree0(self): + pad = padua_points(0) + expected = [[-1],[-1]] + assert_array_almost_equal(pad, expected, 15) + + def test_padua_points_degree1(self): + pad = padua_points(1) + expected = [cos(np.r_[0,1,1]*pi), + cos(np.r_[1,0,2]*pi/2)] + + assert_array_almost_equal(pad, expected, 15) + + def test_padua_points_degree2(self): + pad = padua_points(2, domain=[0,1,0,2]) + expected = [(cos(np.r_[0,0,1,1,2,2]*pi/2)+1)/2, + cos(np.r_[1,3,0,2,1,3]*pi/3)+1] + + assert_array_almost_equal(pad, expected, 15) + + def test_testfunct(self): + vals = [testfunct(0, 0, id) for id in range(12)] + expected = [7.664205912849231e-01, 0.7071067811865476, 0, + 1.6487212707001282, 1.9287498479639178e-22, 1.0, + 1.0, 1.0, 1.0, 0.0, 1.0, 0.0] + assert_array_almost_equal(vals, expected, 15) + + def test_padua_fit_even_degree(self): + points = padua_points(10) + C0f, abs_error = padua_fit(points, testfunct, 6) + expected = np.zeros((11, 11)) + expected[0,0] = 1; + assert_array_almost_equal(C0f, expected, 15) + + def test_padua_fit_odd_degree(self): + points = padua_points(9) + C0f, abs_error = padua_fit(points, testfunct, 6) + expected = np.zeros((10, 10)) + expected[0,0] = 1; + assert_array_almost_equal(C0f, expected, 15) + + def test_padua_fit_odd_degree2(self): + points = padua_points(9) + C0f, abs_error = padua_fit2(points, testfunct, 6) + expected = np.zeros((10, 10)) + expected[0,0] = 1; + assert_array_almost_equal(C0f, expected, 15) + + def test_padua_cubature(self): + domain = [0,1,0,1] + points = padua_points(500, domain) + C0f, abs_error = padua_fit(points, testfunct, 0) + val = padua_cubature(C0f, domain) + expected = 4.06969589491556e-01 + assert_array_almost_equal(val, expected, 15) + + def test_padua_val_unordered(self): + domain = [0,1,0,1] + points = padua_points(20, domain) + C0f, abs_error = padua_fit(points, testfunct, 0) + X = [0,0.5,1] + + val = padua_val(X, X, C0f, domain) + expected = [7.664205912849228e-01, 3.2621734202884815e-01, + 3.587865112678535e-02] + assert_array_almost_equal(val, expected, 14) + + def test_padua_val_grid(self): + domain = [0,1,0,1] + a, b, c, d = domain + points = padua_points(21, domain) + C0f, abs_error = padua_fit(points, testfunct, 0) + X1 = np.linspace(a, b, 2) + X2 = np.linspace(c, d, 2) + val = padua_val(X1, X2, C0f,domain, use_meshgrid=True); + expected = [[7.664205912849229e-01,1.0757071952145181e-01], + [2.703371615911344e-01,3.5734971024838565e-02]] + assert_array_almost_equal(val, expected, 14)