danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added padua.py

Browse Source
master
Per.Andreas.Brodtkorb 10 years ago
parent c93efca266
commit 0a23402fc6
2 changed files with 575 additions and 0 deletions
Show Stats Download Patch File Download Diff File

486
pywafo/src/wafo/padua.py
Unescape Escape View File

@ -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)

89
pywafo/src/wafo/test/test_padua.py
Unescape Escape View File

@ -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)
Write Preview
Loading…
Cancel
Save