[Scipy-svn] r4587 - in branches/Interpolate1D: . docs

scipy-svn@scip... scipy-svn@scip...
Thu Jul 31 14:09:03 CDT 2008


Author: fcady
Date: 2008-07-31 14:09:00 -0500 (Thu, 31 Jul 2008)
New Revision: 4587

Modified:
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/interpolate1d.py
Log:
documentation improved, tweaked the api in the case of an instance of a callable class being passed

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-07-31 16:58:25 UTC (rev 4586)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-07-31 19:09:00 UTC (rev 4587)
@@ -2,25 +2,25 @@
 Overview
 ==================
 
-Often researchers need to infer from a data set the values at points where measurements
-were not taken.  Perhaps they are running a simulation that will demand data points
-never measured.  Perhaps to estimate a statistic of a function (say, the integral) they
-want to guess its value everywhere based on a few points.  Perhaps they know a function
-which must be evaluated multiple times, but evaluating it is expensive; they want to evaluate
-it several times at first, and make educated guesses in the future.  In all these cases, interpolation
-it the tool of choice.
+Interpolation is the process of using known data to guess the value of unknown data.  It turns
+a sample of a function into an approximation of that function for every value, and is one of
+the most basic mathematical tools available to a researcher.
 
 The interpolate package provides tools for interpolating and extrapolating new data points from a known set of data points.  
 Interpolate provides both a functional interface that is flexible and easy to use as well as an object oriented interface that 
 can be more efficient and flexible for some cases.  It is able to interpolate and extrapolate in 1D, 2D, and even N 
-dimensions.*[fixme: 1D only right now]*
+dimensions. *[FIXME : 1D only right now]*
 
 For 1D interpolation, it handles linear and spline(cubic, quadratic, and quintic) for both uniformly and non-uniformly spaced 
 data points "out of the box."  Users can choose the behavior when new values fall outside the range of known data, and
 with a little more work, they can incorporate interpolation methods that are specially tailored to their needs.
 
-For 2D interpolation, 
+For 2D interpolation, *[FIXME : include this]*
 
+This tutorial covers how to use the interpolate module, provides some basic examples, and shows
+them at work in realistic sample sessions.  These sessions demonstrate how to use the 
+interpolate module, but also highlight some of the uses of interpolation techniques.
+
 ================================================
 1D Interpolation with the Functional Interface
 ================================================
@@ -38,36 +38,32 @@
 The following example uses the 'interp1d' function to linearly interpolate a sin 
 curve from a sparse set of values. ::
     
-	# start up ipython for our examples.
-	$ ipython -pylab
-	
-	In [1]: from interpolate import interp1d
-	
-	# Create our "known" set of 5 points with the x values in one array and the y values in another.
-	In [2]: x = linspace(0, 2*pi, 5)
-	In [3]: y = sin(x)
+    # start up ipython for our examples.
+    $ ipython -pylab
     
-::
+    In [1]: from interpolate import interp1d
     
+    # Create our "known" set of 5 points with the x values in one array and the y values in another.
+    In [2]: x = linspace(0, 2*pi, 5)
+    In [3]: y = sin(x)
+    
     # If we only want a value at a single point, we can pass in a scalar and interp1d
     # will return a scalar
-    In [9]: interp1d(x, y, 1.2)
-    Out [10]: 0.76394372684109768
+    In []: interp1d(x, y, 1.2)
+    Out []: 0.76394372684109768
     
-::
-    
     # 0-dimensional arrays are also treated as scalars
-    In [9]: interp1d(x, y, array(1.2) )
-    Out [10]: 0.76394372684109768
+    In []: interp1d(x, y, array(1.2) )
+    Out []: 0.76394372684109768
     
     # To interpolate from these x,y values at multiple points, possibly to get a more dense set
     # of new_x, new_y values to approximate the function, pass a numpy array to interp1d, 
     # and the return type will also be a numpy array.
-	In [4]: new_x = linspace(0, 2*pi, 21)
-	In [5]: new_y = interp1d(x, y, new_x)
-
-	# Plot the results using matplotlib. [note examples assume you are running in ipython -pylab]
-	In [6]: plot(x, y, 'ro', new_x, new_y, 'b-')
+    In []: new_x = linspace(0, 2*pi, 21)
+    In []: new_y = interp1d(x, y, new_x)
+    
+    # Plot the results using matplotlib. [note examples assume you are running in ipython -pylab]
+    In []: plot(x, y, 'ro', new_x, new_y, 'b-')
         
 .. image:: interp1d_linear_simple.png 
 
@@ -83,15 +79,15 @@
 
     # If we attempt to extrapolate values outside the interpolation range, interp1d defaults
     # to returning NaN
-    In [7]: interp1d(x, y, array([-2, -1, 1, 2]))
-    Out [8]: array([        NaN,     NaN,     0.63661977,   0.72676046])
+    In []: interp1d(x, y, array([-2, -1, 1, 2]))
+    Out []: array([        NaN,     NaN,     0.63661977,   0.72676046])
 
 If we want a type of interpolation other than linear, there is a range of options which we can specify 
-with the keyword argument kind, which is usually a string.  For example::
+with the keyword argument "kind", which is usually a string.  Continuing from the previous example,::
 
     # If we want quadratic (2nd order) spline interpolation, we can use the string 'quadratic'
-    In [7]: new_y_quadratic = interp1d(x, y, new_x, kind = 'quadratic')
-    In [8]: plot(x, y, 'r', new_x, new_y_quadratic, 'g')
+    In []: new_y_quadratic = interp1d(x, y, new_x, kind = 'quadratic')
+    In []: plot(x, y, 'r', new_x, new_y_quadratic, 'g')
     
 .. image:: interp1d_linear_and_quadratic.png
 
@@ -105,7 +101,7 @@
 #. 'quintic' : 5th order spline interpolation
 
 The same flexibility is afforded for extrapolation by the keywords low and high, which
-specify how to treat values below and above the range of know data: ::
+specify how to treat values below and above the range of known data: ::
 
     In []: z = array([ 1.0, 2.0 ])
     In []: interp1d(z, z, array([-5.0, 5.0]), low = 'linear', high = 'linear') # -5 and 5 both out of range
@@ -133,7 +129,7 @@
 Removal of Bad Datapoints
 -----------------------------
 
-Many datasets have missing or corrupt data which it is desirable to ignore when interpolating,
+Many datasets have missing or corrupt data which we want to ignore when interpolating,
 and to this end, interp1d has the keyword argument bad_data.
 
 bad_data defaults to being None.  But if it is a list, all "bad" points (x[i], y[i]) will be removed
@@ -144,7 +140,7 @@
 Note that bad_data must be either None or a list of numbers.  Including NaN or None in the list,
 for example, is not supported and will cause errors. 
 
-The following example shows how ::
+The following example demonstrates using this keyword argument ::
 
     # data will be linear, except for artificial bad points
     In []: x = arange(10.); y = arange(10.)
@@ -165,12 +161,12 @@
 
 The string interface is designed to conveniently take care of most things a user would want
 to do in a way that is easy and, when something goes wrong, informative and helpful.
-If, however, you want more direct control than is afforded by the string interface, that is also possible,
-thought it's a little trickier than using strings.  You must be very careful to have correct
-format, and failure to do so can cause a range of errors which won't necessarily result in
+If, however, you want more direct control than is afforded by the string interface, that is also possible.
+If you define your own types, you must be very careful to have correct
+format; failure to do so can cause a range of errors which won't necessarily result in
 informative error messages.
 
-kind (or, equivalently, low and high) can also be set to a function, a callable 
+kind (or, equivalently, low or high) can also be set to a function, a callable 
 class, or an instance of a callable class.
 
 If a function is passed, it will be called when interpolating.
@@ -178,36 +174,45 @@
 
         newy = interp(x, y, newx)
         
-where x, y, newx, and newy are all numpy arrays.
+where x, y, newx, and newy are all 1D numpy arrays.
             
-If a callable class is passed, it is assumed to have format::
+If a class is passed, it is assumed to have ones of two formats.
+If there is a "init_xy" or "set_xy" method, the class is instantiated
+with no argument, then the relevant method is called to initialize 
+x and y, and the class is later called with a 1D array as an argument.::
 
-        instance = Class(x, y).
-        
-which can then be called by ::
+        instance = Class().
+        instance.set_xy(x, y)
+        new_y = instance(new_x)
 
+If the class does not have an init_xy or set_xy method, the class
+is instantiated with x and y as arguments, and passed a 1D array
+during interpolation. ::
+
+            instance = Class(x, y)
             new_y = instance(new_x)
             
-If a callable object with method "init_xy" or "set_xy" is
-passed, that method will be used to set x and y as follows: ::
+You can also pass an instance of the callable class, rather than the class
+itself.  This is useful is the class has other parameters besides x, y, and
+new_x (perhaps smoothing coefficients, orders for polynomials, etc).
 
+If the instance has a method "init_xy" or "set_xy", 
+that method will be used to set x and y, and the instance will be
+called later: ::
+
         instance.set_xy(x, y)
-        
-and the object will be called during interpolation. ::
-
         new_y = instance(new_x)
                 
-If the "init_xy" and "set_xy" are not present, it will be called as
+If the instance has no "init_xy" or "set_xy" method, it will be called as ::
 
-        new_y = argument(new_x)
-                
-A primitive type which is not a string signifies a function
-which is identically that value (e.g. val and 
-lambda x, y, newx : val are equivalent). ::
+        new_y = kind(x, y, new_x)
+        
+Failure to follow these guidelines (say, by having kind require other keyword
+arguments, having a method "initialize_xy" rather than "init_xy", etc) can result
+in cryptic errors, so be careful.  Here is a demo of how to properly use these features:
 
-    # However, this behavior can be overwritten in the same way as linear interpolation,
-    # by setting the keyword extrap_low (for values below the range of interpolation) and
-    # extrap_high (for values above that range)
+::
+
     In []: def dummy(x, y, newx):
                 # Note that dummy has acceptable form
                 return 5.7
@@ -222,10 +227,10 @@
     In []: y = arange(5.0)
     In []: new_x = np.array([ -1, 2.4, 7 ])
     In []: new_y = interp1d(x, y, 
-                                        kind = Phony, 
-                                        low = dummy,
-                                        high = dummy
-                                        )
+                            kind = Phony, 
+                            low = dummy,
+                            high = dummy
+                            )
     In []: new_y
     Out []: array([ 5.7, 4.0, 5.7 ])
 
@@ -251,7 +256,7 @@
 
     # The default behavior is virtually the same
     In []: x = linspace(0, 2*pi, 5)
-	In []: y = sin(x)
+    In []: y = sin(x)
     In []: new_x = linspace(0, 2*pi, 21)
     In []: new_y1 = interp1d(x, y, new_x)
     In []: interp_obj1 = Interpolate1d(x, y)
@@ -279,16 +284,12 @@
 Estimating Function Statistics and Displaying Data
 -----------------------------------------------------
 
-In this session, the engineer
-has a data set of geological data indicating the temperature at various
-depths in the ground.  The data set is noisy and large.  He wants to 1)
-get a feel for the data, and 2) estimate the average temperature.
+In this session, the geologist
+has a data set of data indicating the temperature at various
+depths in the ground.  He wants to 1) get a visual feel for the data, 
+and 2) estimate the average temperature.
 ::
 
-    # start up ipython for our examples.
-	$ ipython -pylab
-    
-    # load the data from a text file
     In []: data_array = loadtxt('dataset1.txt')
     In []: shape(data_array)
     Out []: (12, 2)
@@ -305,26 +306,29 @@
     
     In []: import interpolate as I
     In []: plot( I.interp1d(depth, temp, linspace(0,20,100), bad_data = [1000])
-    In []: # much better, but he wants to see it smoother too
+    # much better, but he wants to see it smoother too
     In []: plot( I.interp1d(depth, temp, linspace(0,20,100), kind='cubic', bad_data = [1000])
     
-    # To find the average temp he can't average the data because the samples
-    # are not necessarily uniform, but he can uniformly sample the interpolated function
-    In []: average_temp = average( I.interp1d(depth, temp, linspace(0,20,100), 'cubic', bad_data=[1000]))
+    # To find the average temp he can't average the data points because the samples
+    # are not necessarily uniform, but it is easy to uniformly sample the interpolated function
+    In []: average_temp = average( I.interp1d(depth, temp, linspace(0,20,100), 'cubic', bad_data=[1000]) )
     
 ---------------------------------
-Modelling from a small dataset
+Modeling from a small dataset
 ---------------------------------
 
 This computational biologist wants to model the growth rate of 
-cancer cells in tissue.  He has measurements of the metabolic rate of cancer 
-cells at several concentrations of blood glucose.  He also has measurements
-of the growth rate of these cells as a function of their CO2 metabolic output.  Each data point represents 
+cancer cells in tissue.  For several levels of blood glucose, he has measurements 
+of the CO2 output of the cancer cells. For several different levels of CO2 ouput,
+he also has measurements of the growth rate of these cells.  Each data point represents 
 a week's work on the part of experimentalists, so though there isn't much 
 data he'll have to make due.  Now, his full simulation takes up hundreds of lines of
-code, but we show here the module object estimate_growth_rate which he wrote. 
+code, so we only show the module estimate_growth_rate.py which he wrote. 
 ::
 
+    """ Contains callable class EstimateGrowthRate, which accepts blood glucose level as
+        an argument and returns interpolated growth rate of cells.
+    """
     import numpy as np
     import interpolate as I
     
@@ -387,8 +391,6 @@
     In []: where_to_insert = max( find(thickness < guess_thick) ) +1
     In []: thickness = insert(thickness, where_to_insert, guess_thick)
     In []: peformance = insert(performance, where_to_insert, measured_perf)
-
-
-            
-
-        
\ No newline at end of file
+    
+More sophisticated optimization tools are also available from the scipy.optimize
+module.
\ No newline at end of file

Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py	2008-07-31 16:58:25 UTC (rev 4586)
+++ branches/Interpolate1D/interpolate1d.py	2008-07-31 19:09:00 UTC (rev 4587)
@@ -1,79 +1,92 @@
 # FIXME: information strings giving mathematical descriptions of the actions
 #     of the functions.
 
-from interpolate_wrapper import linear, logarithmic, block, block_average_above, atleast_1d_and_contiguous
+from interpolate_wrapper import linear, logarithmic, block, \
+                block_average_above, atleast_1d_and_contiguous
 from fitpack_wrapper import Spline
 import numpy as np
 from numpy import array, arange, empty, float64, NaN
 
-# dictionary of tuples.  First element is a callable (class, instance of a class, or function
-# second argument is dictionary of additional keywords, if any
-dict_of_interp_types = \
-                { 'linear' : (linear, {}), 
-                    'logarithmic' : (logarithmic, {}), 
-                    'block' : (block, {}),
-                    'block_average_above' : (block_average_above, {}),
-                    'Spline' : (Spline, {}), 'spline' : (Spline, {}),
-                    'Quadratic' : (Spline, {'k':2}), 'quadratic' : (Spline, {'k':2}),
-                    'Quad' : (Spline, {'k':2}), 'quad' : (Spline, {'k':2}),
-                    'Cubic' : (Spline, {'k':3}), 'cubic' : (Spline, {'k':3}),
-                    'Quartic' : (Spline, {'k':4}), 'quartic' : (Spline, {'k':4}),
-                    'Quar' : (Spline, {'k':4}), 'quar' : (Spline, {'k':4}),
-                    'Quintic' : (Spline, {'k':5}), 'quintic' : (Spline, {'k':5}),
-                    'Quin' : (Spline, {'k':5}), 'quin' : (Spline, {'k':5})
+# dictionary of interpolation functions/classes/objects
+method_register = \
+                { 'linear' : linear, 
+                    'logarithmic' : logarithmic, 
+                    'block' : block, 
+                    'block_average_above' : block_average_above, 
+                    'Spline' : Spline, 'spline' : Spline,
+                    'Quadratic' : Spline(k=2), 'quadratic' : Spline(k=2),
+                    'Quad' : Spline(k=2), 'quad' : Spline(k=2),
+                    'Cubic' : Spline(k=3), 'cubic' : Spline(k=3),
+                    'Quartic' : Spline(k=3), 'quartic' : Spline(k=3),
+                    'Quar' : Spline(k=4), 'quar' : Spline(k=4),
+                    'Quintic' : Spline(k=5), 'quintic' : Spline(k=5),
+                    'Quin' : Spline(k=5), 'quin' : Spline(k=5)
                 }
+                
+# dictionary of types for casting.  key = possible datatype, value = datatype it is cast to
+# BEWARE : if you cast things to integers, you will lose interpolation ability
+dtype_register = {np.float32 : np.float32, 
+                            np.float64 : np.float64
+                            }
+dtype_default = np.float64
 
 def interp1d(x, y, new_x, 
-                    interp = 'linear', extrap_low = NaN, extrap_high = NaN,
-                    interpkw = {}, lowkw = {}, highkw ={},
+                    kind = 'linear', low = NaN, high = NaN,
                     bad_data = None):
     # FIXME : all y to be multi-dimensional
-    # FIXME : update the doc string to match that of Interpolate1d
-    """ A function for interpolation of 1D data.
+    # NOTE : This docstring is considered suboordinate to that for Interpolate1d.
+    #       That is, update Interpolate1d and copy-and-paste
+    """ A callable class for interpolation of 1D, real-valued data.
         
         Parameters
         -----------
             
-        x -- list or NumPy array
-            x includes the x-values for the data set to
-            interpolate from.  It must be sorted in
-            ascending order
+            x -- list or 1D NumPy array
+                x includes the x-values for the data set to
+                interpolate from.  It must be sorted in
+                ascending order.
+                    
+            y -- list or 1D NumPy array
+                y includes the y-values for the data set  to
+                interpolate from.  Note that 2-dimensional
+                y is not currently supported.
                 
-        y -- list or NumPy array
-            y includes the y-values for the data set  to
-            interpolate from.  Note that y must be
-            one-dimensional.
-            
-        new_x -- list or NumPy array
-            points whose value is to be interpolated from x and y.
-            new_x must be in sorted order, lowest to highest.
-                
         Optional Arguments
         -------------------
         
-        interp -- Usu. function or string.  But can be any type.
-            Specifies the type of extrapolation to use for values within
-            the range of x.  If a string is passed, it will look for an object
-            or function with that name and call it when evaluating.  If 
-            a function or object is passed, it will be called when interpolating.
-            If nothing else, assumes the argument is intended as a value
-            to be returned for all arguments.  Defaults to linear interpolation.
+            kind -- Usually a string.  But can be any type.
+                Specifies the type of interpolation to use for values within
+                the range of x.
+                
+                By default, linear interpolation is used.
+                
+                See below for details on other options.
+                
+            low  -- same as for kind
+                How to extrapolate values for inputs below the range of x.
+                Same options as for 'kind'.  Defaults to returning numpy.NaN ('not 
+                a number') for all values below the range of x.
+                
+            high  -- same as for kind
+                How to extrapolate values for inputs above the range of x.
+                Same options as for 'kind'.  Defaults to returning numpy.NaN ('not 
+                a number') for all values above the range of x.
+                
+            bad_data -- list
+                List of numerical values (in x or y) which indicate unacceptable data. 
+                
+                If bad_data is not None (its default), all points whose x or y coordinate is in
+                bad_data, OR ones of whose coordinates is NaN, will be removed.
+                
+            interpkw -- dictionary
+                If interp is set to a function, class or callable object, this contains
+                additional keywords.
+                
+            lowkw (highkw) -- dictionary
+                like interpkw, but for extrap_low and extrap_high
+                
             
-        low (high) -- same as for interp
-            Same options as for 'interp'.  Defaults to returning numpy.NaN ('not 
-            a number') for all values outside the range of x.
-        
-        interpkw -- dictionary
-            If 
-            
-        bad_data -- list
-            List of values (in x or y) which indicate unacceptable data. All points
-            that have x or y value in missing_data will be removed before
-            any interpolatin is performed if bad_data is not None.
-            
-            numpy.NaN is always considered bad data.
-            
-        Acceptable Input Strings
+        Some Acceptable Input Strings
         ------------------------
         
             "linear" -- linear interpolation : default
@@ -82,8 +95,39 @@
             "block_average_above' -- block average above
             "Spline" -- spline interpolation.  keyword k (defaults to 3) 
                 indicates order of spline
-            numpy.NaN -- return numpy.NaN
+            "quad", "quadratic" -- spline interpolation order 2
+            "cubic" -- spline interpolation order 3
+            "quartic" -- spline interpolation order 4
+            "quintic" -- spline interpolation order 5
+            
+        Other options for interp, extrap_low, and extrap_high
+        ---------------------------------------------------
         
+            If you choose to use a non-string argument, you must
+            be careful to use correct formatting.
+            
+            If a function is passed, it will be called when interpolating.
+            It is assumed to have the form 
+                newy = interp(x, y, newx, **kw), 
+            where x, y, newx, and newy are all numpy arrays.
+            
+            If a callable class is passed, it is assumed to have format
+                instance = Class(x, y, **kw).
+            which can then be called by
+                new_y = instance(new_x)
+            
+            If a callable object with method "init_xy" or "set_xy" is
+            passed, that method will be used to set x and y as follows
+                instance.set_xy(x, y, **kw)
+            and the object will be called during interpolation.
+                new_y = instance(new_x)
+            If the "init_xy" and "set_xy" are not present, it will be called as
+                new_y = argument(new_x)
+                
+            A primitive type which is not a string signifies a function
+            which is identically that value (e.g. val and 
+            lambda x, y, newx : val are equivalent).
+            
         Examples
         ---------
         
@@ -96,12 +140,9 @@
             array([.2, 2.3, 5.6, NaN])
     """
     return Interpolate1d(x, y, 
-                                interp = interp,
-                                extrap_low = extrap_low,
-                                extrap_high = extrap_high,
-                                interpkw = interpkw,
-                                lowkw = lowkw,
-                                highkw = highkw,
+                                kind = kind,
+                                low = low,
+                                high = high,
                                 bad_data = bad_data
                                 )(new_x)
 
@@ -124,7 +165,7 @@
         Optional Arguments
         -------------------
         
-            interp -- Usually a string.  But can be any type.
+            kind -- Usually a string.  But can be any type.
                 Specifies the type of interpolation to use for values within
                 the range of x.
                 
@@ -132,12 +173,12 @@
                 
                 See below for details on other options.
                 
-            extrap_low  -- same as for kind
+            low  -- same as for kind
                 How to extrapolate values for inputs below the range of x.
                 Same options as for 'kind'.  Defaults to returning numpy.NaN ('not 
                 a number') for all values below the range of x.
                 
-            extrap_high  -- same as for kind
+            high  -- same as for kind
                 How to extrapolate values for inputs above the range of x.
                 Same options as for 'kind'.  Defaults to returning numpy.NaN ('not 
                 a number') for all values above the range of x.
@@ -148,14 +189,6 @@
                 If bad_data is not None (its default), all points whose x or y coordinate is in
                 bad_data, OR ones of whose coordinates is NaN, will be removed.
                 
-            interpkw -- dictionary
-                If interp is set to a function, class or callable object, this contains
-                additional keywords.
-                
-            lowkw (highkw) -- dictionary
-                like interpkw, but for extrap_low and extrap_high
-                
-            
         Some Acceptable Input Strings
         ------------------------
         
@@ -170,7 +203,7 @@
             "quartic" -- spline interpolation order 4
             "quintic" -- spline interpolation order 5
             
-        Other options for interp, extrap_low, and extrap_high
+        Other options for kind, low, and high
         ---------------------------------------------------
         
             If you choose to use a non-string argument, you must
@@ -178,17 +211,17 @@
             
             If a function is passed, it will be called when interpolating.
             It is assumed to have the form 
-                newy = interp(x, y, newx, **kw), 
+                newy = interp(x, y, newx), 
             where x, y, newx, and newy are all numpy arrays.
             
             If a callable class is passed, it is assumed to have format
-                instance = Class(x, y, **kw).
+                instance = Class(x, y).
             which can then be called by
                 new_y = instance(new_x)
             
             If a callable object with method "init_xy" or "set_xy" is
             passed, that method will be used to set x and y as follows
-                instance.set_xy(x, y, **kw)
+                instance.set_xy(x, y)
             and the object will be called during interpolation.
                 new_y = instance(new_x)
             If the "init_xy" and "set_xy" are not present, it will be called as
@@ -211,30 +244,32 @@
             array([.2, 2.3, 5.6, NaN])
             
     """
-    # FIXME: more informative descriptions of sample arguments
-    # FIXME: examples in doc string
-    # FIXME : Allow copying or not of arrays.  non-copy + remove_bad_data should flash 
-    #           a warning (esp if we interpolate missing values), but work anyway.
     
     def __init__(self, x, y, 
-                        interp = 'linear', 
-                        extrap_low = NaN, 
-                        extrap_high = NaN,
-                        interpkw = {},
-                        lowkw = {},
-                        highkw = {},
+                        kind = 'linear', 
+                        low = NaN, 
+                        high = NaN,
                         bad_data = None):
-        # FIXME: don't allow copying multiple times.
-        # FIXME : allow no copying, in case user has huge dataset
         
-        # remove bad data, is there is any
+        # put data into nice format and store it
+        self._init_xy(x, y, bad_data)
+        
+        # store interpolation functions for each range
+        self.interp = self._init_interp_method(kind)
+        self.extrap_low = self._init_interp_method(low)
+        self.extrap_high = self._init_interp_method(high)
+
+    def _init_xy(self, x, y, bad_data):
+        # FIXME : no-copying option, in case user has huge dataset.  non-copy + remove bad data should
+        #               flash a warning (esp if, in the future, bad values are interpolated).
+        
+        # remove bad data if applicable
         if bad_data is not None:
-            try:
+            try: # check that bad_data contains only numerical values
                 sum_of_bad_data = sum(bad_data)
             except:
                 raise TypeError, "bad_data must be either None \
-                        or a list of numerical types"
-            
+                        or a list of numbers"            
             x, y = self._remove_bad_data(x, y, bad_data)
         
         # check acceptable size and dimensions
@@ -245,20 +280,9 @@
         assert y.ndim == 1 , "y must be one-dimensional" 
         assert len(x) == len(y) , "x and y must be of the same length"
         
-        # store data, and remove bad data points is applicable
-        # FIXME : may be good to let x and y be initialized later, or changed after-the-fact
-        self._init_xy(x, y)
-        
-        # store interpolation functions for each range
-        self.interp = self._init_interp_method(interp, interpkw)
-        self.extrap_low = self._init_interp_method(extrap_low, lowkw)
-        self.extrap_high = self._init_interp_method(extrap_high, highkw)
-
-    def _init_xy(self, x, y):
-        
         # select proper dataypes and make arrays
-        self._xdtype = {np.float32 : np.float32}.setdefault(type(x[0]), np.float64) # unless data is float32,  cast to float64
-        self._ydtype = {np.float32 : np.float32}.setdefault(type(y[0]), np.float64)
+        self._xdtype = dtype_register.setdefault(type(x[0]), dtype_default)
+        self._ydtype = dtype_register.setdefault(type(y[0]), dtype_default)
         self._x = atleast_1d_and_contiguous(x, self._xdtype).copy()
         self._y = atleast_1d_and_contiguous(y, self._ydtype).copy()
 
@@ -266,9 +290,6 @@
         """ removes data points whose x or y coordinate is
             either in bad_data or is a NaN.
         """
-        # FIXME : In the future, it may be good to just replace the bad points with good guesses.
-        #       Especially in generalizing the higher dimensions
-        # FIXME : This step is very inefficient because it iterates over the array
         
         bad_data_mask = np.isnan(x) | np.isnan(y)
         for bad_num in bad_data:
@@ -278,20 +299,15 @@
         y = y[~bad_data_mask]
         return x, y
         
-    def _init_interp_method(self, interp_arg, kw):
-        """
-            returns the interpolating function specified by interp_arg.
-        """
-        # FIXME : error checking specific to interpolation method.  x and y long
-        #   enough for order-3 spline if that's indicated, etc.  Functions should throw
-        #   errors themselves, but errors at instantiation would be nice.
-        
+    def _init_interp_method(self, interp_arg):
+        """ returns the interpolating function specified by interp_arg.
+        """        
         from inspect import isclass, isfunction
         
         # primary usage : user passes a string indicating a known function
+        # pick interpolator accordingly
         if isinstance(interp_arg, basestring):
-            interpolator, kw = dict_of_interp_types.setdefault(interp_arg, (None, {}) )
-            
+            interpolator = method_register.setdefault(interp_arg, None )
             if interpolator is None: 
                 raise TypeError, "input string %s not valid" % interp_arg
         else:
@@ -301,29 +317,29 @@
         if hasattr(interpolator, '__call__'):
             # function
             if isfunction(interpolator):
-                result = lambda newx : interpolator(self._x, self._y, newx, **kw)
+                result = lambda newx : interpolator(self._x, self._y, newx)
                 
             # callable class 
             elif isclass(interpolator):
                 if hasattr(interpolator, 'set_xy'):
-                    result = interpolator(**kw)
+                    result = interpolator()
                     result.set_xy(self._x, self._y)
                 if hasattr(interpolator, 'init_xy'):
-                    result = interpolator(**kw)
+                    result = interpolator()
                     result.init_xy(self._x, self._y)
                 else:
-                    result = interpolator(self._x, self._y, **kw)
+                    result = interpolator(self._x, self._y)
                 
             # instance of callable class
             else:
                 if hasattr(interpolator, 'init_xy'):
                     result = interpolator
-                    result.init_xy(self._x, self._y, **kw)
+                    result.init_xy(self._x, self._y)
                 elif hasattr(interpolator, 'set_xy'):
                     result = interpolator
-                    result.set_xy(self._x, self._y, **kw)
+                    result.set_xy(self._x, self._y)
                 else:
-                    result = interpolator
+                    result = lambda new_x : interpolator(self._x, self._y, new_x)
             
         # non-callable : user has passed a default value to always be returned
         else:
@@ -332,49 +348,38 @@
         return result
 
     def __call__(self, newx):
-        """
-            Input x must be a list or NumPy array in sorted order.
+        """ Input x must be a list or NumPy array in sorted order.
             
             Breaks x into pieces in-range, below-range, and above range.
             Performs appropriate operation on each and concatenates results.
         """
-        # FIXME : atleast_1d_and_contiguous may also be called within the interpolation technique.
-        #   waste of time, but ok for the time being.
         
-        # if input is scalar or 0-dimemsional array, output will be scalar
+        # record if input is scalar or 0-dimemsional array, in which case output will be scalar
         input_is_scalar = np.isscalar(newx) or \
                                     (
                                         isinstance(  newx , np.ndarray  ) and 
                                         np.shape(newx) == ()
                                     )
         
-        # make 
-        newx_array = atleast_1d_and_contiguous(newx)
+        # make input into a nice 1d, contiguous array
+        newx_array = atleast_1d_and_contiguous(newx, dtype=self._xdtype)
+        assert newx_array.ndim == 1, "new_x can be at most 1-dimensional"
         
         # masks indicate which elements fall into which interpolation region
         low_mask = newx_array<self._x[0]
         high_mask = newx_array>self._x[-1]
         interp_mask = (~low_mask) & (~high_mask)
                 
-        type(newx_array[low_mask])
-                
-                
-        # use correct function for x values in each region
-        if len(newx_array[low_mask]) == 0: new_low=np.array([])  # FIXME : remove need for if/else.
-                                                                            # if/else is a hack, since vectorize is failing
-                                                                            # to work on lists/arrays of length 0
-                                                                            # on the computer where this is being
-                                                                            # developed
+        # use correct function for x values in each region and create output as an array
+        if len(newx_array[low_mask]) == 0: new_low=np.array([])  # FIXME : remove need for if/else hack.
+                                                                            # it's there since vectorize is failing on arrays of zero length
         else: new_low = self.extrap_low(newx_array[low_mask])
         if len(newx_array[interp_mask])==0: new_interp=np.array([])
         else: new_interp = self.interp(newx_array[interp_mask])
         if len(newx_array[high_mask]) == 0: new_high = np.array([])
         else: new_high = self.extrap_high(newx_array[high_mask])
+        result_array = np.concatenate((new_low, new_interp, new_high))
         
-        result_array = np.concatenate((new_low, new_interp, new_high)) # FIXME : deal with mixed datatypes
-                                                                                          # Would be nice to say result = zeros(dtype=?)
-                                                                                          # and fill in
-        
         # convert to scalar if scalar was passed in
         if input_is_scalar:
             result = float(result_array)



More information about the Scipy-svn mailing list