Source Code for Module pylal.chebyshev_interpolation

  1  # Copyright (C) 2011  Drew Keppel 
  2  # 
  3  # This program is free software; you can redistribute it and/or modify it 
  4  # under the terms of the GNU General Public License as published by the 
  5  # Free Software Foundation; either version 2 of the License, or (at your 
  6  # option) any later version. 
  7  # 
  8  # This program is distributed in the hope that it will be useful, but 
  9  # WITHOUT ANY WARRANTY; without even the implied warranty of 
 10  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
 11  # Public License for more details. 
 12  # 
 13  # You should have received a copy of the GNU General Public License along 
 14  # with this program; if not, write to the Free Software Foundation, Inc., 
 15  # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA. 
 16   
 17  import scipy 
 18  import numpy 
 19   
 20 -def U_cheby_nodes_2D(N1, N2): 
 21          """ 
 22          A function to compute the positions (x,y) of the nodes of the 
 23          2-dimensional Chebyshev polynomials of degree (N1,N2). 
 24          """ 
 25   
 26          x_cheby = numpy.array([numpy.cos(numpy.pi*(2*idx+1)/N1/2) for idx in range(N1)]) 
 27          y_cheby = numpy.array([numpy.cos(numpy.pi*(2*idx+1)/N2/2) for idx in range(N2)]) 
 28          return scipy.meshgrid(x_cheby, x_cheby) 
 29   
 30   
 31 -def factorial(x): 
 32          """ 
 33          A wrapper for scipy's factorial. 
 34          """ 
 35   
 36          return scipy.factorial(x,exact=1) 
 37   
 38 -def T_cheby_2D(x, y, n1, N1, n2, N2): 
 39          """ 
 40          A function that returns the (n,m)th 2-dimensional Chebyshev 
 41          polynomials of degree (N,M) at the locations given in x and y. 
 42          """ 
 43   
 44          ux = numpy.zeros(scipy.shape(x)) 
 45          uy = numpy.zeros(scipy.shape(y)) 
 46          for thisn1 in range(n1/2+1): 
 47                  ux += (x**2.-1.)**(thisn1) * x**(n1-2*thisn1) * factorial(n1)/factorial(2*thisn1)/factorial(n1-2*thisn1) 
 48          for thisn2 in range(n2/2+1): 
 49                  uy += (y**2.-1.)**(thisn2) * y**(n2-2*thisn2) * factorial(n2)/factorial(2*thisn2)/factorial(n2-2*thisn2) 
 50   
 51          w = 1. 
 52          if n1: 
 53                  w *= N1/2. 
 54          else: 
 55                  w *= N1 
 56          if n2: 
 57                  w *= N2/2. 
 58          else: 
 59                  w *= N2 
 60          return ux*uy/w**.5 
 61   
 62 -def interpolation_kernel_values(x, y, N1, N2): 
 63          """ 
 64          A function that creates all of the 2-dimensional Chebyshev polynomials 
 65          of degree (N1,N2) evalulated at positions given by (x,y).  
 66          """ 
 67   
 68          Us = [] 
 69          for idx in range(N1): 
 70                  Us.append([]) 
 71                  for jdx in range(N2): 
 72                          Us[idx].append(T_cheby_2D(x, y, idx, N1, jdx, N2)) 
 73          Us = numpy.array(Us) 
 74          return Us 
 75   
 76 -def interpolation_coefficients(func_measurements, Us, N1, N2): 
 77          """ 
 78          A function that takes function measurements at the interpolation 
 79          locations and computes the interpolation coefficients associated with 
 80          the 2-dimensional Chebyshev polynomials of degree (N1,N2) (Us). 
 81          """ 
 82   
 83          coeffs = [] 
 84          for idx in range(N1): 
 85                  coeffs.append([]) 
 86                  for jdx in range(N2): 
 87                          coeffs[idx].append(sum((func_measurements*Us[idx][jdx]).flatten())) 
 88          return coeffs 
 89   
 90 -def reconstruct_interpolation(coeffs, Vs, N1, N2): 
 91          """ 
 92          A function that takes interpolation coefficients associated with 
 93          the 2-dimensional Chebyshev polynomials of degree (N1,N2) and 
 94          recontructs the interpolated function at locations where the 
 95          2-dimensional Chebyshev polynomials of degree (N1,N2) (Vs) have been sampled. 
 96          """ 
 97   
 98          func_interpolated = numpy.zeros(numpy.shape(Vs[0,0])) 
 99          for idx in range(N1): 
100                  for jdx in range(N2): 
101                          func_interpolated += coeffs[idx][jdx]*Vs[idx,jdx] 
102          return func_interpolated 
103