.. highlight:: python
   :linenothreshold: 5

.. _chap7:

*********
Functions
*********

As you develop more complex computer code, it becomes increasingly important to organize your code into modular blocks.  One important means for doing so is  *user-defined* Python functions.  User-defined functions are a lot like built-in functions that we have encountered in core Python as well as in NumPy and Matplotlib.  The main difference is that user-defined functions are written by you.  The idea is to define functions to simplify your code and to allow you to reuse the same code in different contexts.

The number of ways that functions are used in programming is so varied that we cannot possibly enumerate all the possibilities.  As our use of Python functions in scientific program is somewhat specialized, we introduce only a few of the possible uses of Python functions, ones that are the most common in scientific programming.

.. index::
   single: functions; user defined

.. _userDefdFuncs:

User-defined functions
======================

The NumPy package contains a plethora of mathematical functions.  You can find a listing of the mathematical functions available through NumPy on the web page http://docs.scipy.org/doc/numpy/reference/routines.math.html.  While the list may seem pretty exhaustive, you may nevertheless find that you need a function that is not available in the NumPy Python library.  In those cases, you will want to write your own function.

In studies of optics and signal processing one often runs into the sinc function, which is defined as

.. math::

    \mathrm{sinc}\,x \equiv \frac{\sin x}{x} \;.

Let's write a Python function for the sinc function.  Here is our first attempt:

.. sourcecode:: python
    
    def sinc(x):
        y = np.sin(x)/x
        return y

Every function definition begins with the word ``def`` followed by the name you want to give to the function, ``sinc`` in this case, then a list of arguments enclosed in parentheses, and finally terminated with a colon.  In this case there is only one argument, ``x``, but in general there can be as many arguments as you want, including no arguments at all.  For the moment, we will consider just the case of a single argument.

The indented block of code following the first line defines what the function does.  In this case, the first line calculates :math:`\mathrm{sinc}\,x = \sin x/x` and sets it equal to ``y``.  The ``return`` statement of the last line tells Python to return the value of ``y`` to the user.

We can try it out in the IPython shell.  First we type in the function definition.

.. sourcecode:: ipython
    
    In [1]: def sinc(x):
       ...:     y = sin(x)/x
       ...:     return y

Because we are doing this from the IPython shell, we don't need to import NumPy; it's preloaded.  Now the function :math:`\mathrm{sinc}\,x` is available to be used from the IPython shell

.. sourcecode:: ipython
    
    In [2]: sinc(4)
    Out[2]: -0.18920062382698205
    
    In [3]: a = sinc(1.2)
    
    In [4]: a
    Out[4]: 0.77669923830602194
       
    In [5]: sin(1.2)/1.2
    Out[5]: 0.77669923830602194
    
Inputs and outputs 4 and 5 verify that the function does indeed give the same result as an explicit calculation of :math:`\sin x/x`.

You may have noticed that there is a problem with our definition of :math:`\mathrm{sinc}\,x` when ``x=0.0``.  Let's try it out and see what happens

.. sourcecode:: ipython
        
    In [6]: sinc(0.0)
    Out[6]: nan

IPython returns ``nan`` or "not a number", which occurs when Python attempts a division by zero, which is not defined.  This is not the desired response as :math:`\mathrm{sinc}\,x` is, in fact, perfectly well defined for :math:`x=0`.  You can verify this using L'Hopital's rule, which you may have learned in your study of calculus, or you can ascertain the correct answer by calculating the Taylor series for :math:`\mathrm{sinc}\,x`.  Here is what we get

.. math::

    \mathrm{sinc}\,x = \frac{\sin x}{x} 
                    = \frac{x - \frac{x^3}{3!} + \frac{x^5}{5!} + ...}{x} 
                    = 1 - \frac{x^2}{3!} + \frac{x^4}{5!} + ... \;.

From the Taylor series, it is clear that :math:`\mathrm{sinc}\,x` is well-defined at and near :math:`x=0` and that, in fact,  :math:`\mathrm{sinc}(0)=1`.  Let's modify our function so that it gives the correct value for ``x=0``.

.. sourcecode:: ipython
        
    In [7]: def sinc(x):
       ...:     if x==0.0:
       ...:         y = 1.0
       ...:     else:
       ...:         y = sin(x)/x
       ...:     return y
    
    In [8]: sinc(0)
    Out[8]: 1.0
    
    In [9]: sinc(1.2)
    Out[9]: 0.77669923830602194

Now our function gives the correct value for ``x=0`` as well as for values different from zero.

.. index::
   single: functions; looping over arrays

.. _loopingarrays:

Looping over arrays in user-defined functions
---------------------------------------------

The code for :math:`\mathrm{sinc}\,x` works just fine when the argument is a single number or a variable that represents a single number.  However, if the argument is a NumPy array, we run into a problem, as illustrated below.

.. sourcecode:: ipython
            
    In [10]: x = arange(0, 5., 0.5)
    
    In [11]: x
    Out[11]: array([ 0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ,  3.5,
                     4. ,  4.5])
    
    In [12]: sinc(x)
    ----------------------------------------------------------
    ValueError               Traceback (most recent call last)
    ----> 1 sinc(x)
    
          1 def sinc(x):
    ----> 2     if x==0.0:
          3         y = 1.0
          4     else:
          5         y = np.sin(x)/x
    
    ValueError: The truth value of an array with more than one
                element is ambiguous.
    
The ``if`` statement in Python is set up to evaluate the truth value of a single variable, not of multielement arrays.  When Python is asked to evaluate the truth value for a multi-element array, it doesn't know what to do and therefore returns an error.

An obvious way to handle this problem is to write the code so that it processes the array one element at a time, which you could do using a ``for`` loop, as illustrated below.


.. sourcecode:: python
    :linenos:

    def sinc(x):
        y = []              # creates an empty list to store results
        for xx in x:        # loops over all elements in x array
            if xx==0.0:     # adds result of 1.0 to y list if
                y += [1.0]  # xx is zero
            else:           # adds result of sin(xx)/xx to y list if
                y += [np.sin(xx)/xx]  # xx is not zero
        return np.array(y)  # converts y to array and returns array
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    x = np.linspace(-10, 10, 256)
    y = sinc(x)
    
    plt.plot(x, y)
    plt.axhline(color="gray", zorder=-1)
    plt.axvline(color="gray", zorder=-1)
    plt.show()

The ``for`` loop evaluates the elements of the ``x`` array one by one and appends the results to the list ``y`` one by one.  When it is finished, it converts the list to an array and returns the array.  The code following the function definition plots :math:`\mathrm{sinc}\,x` as a function of :math:`x`.

In the program above, you may have noticed that the NumPy library is imported *after* the ``sinc(x)`` function definition.  As the function uses the NumPy functions ``sin`` and ``array``, you may wonder how this program can work.  Doesn't the ``import numpy`` statement have to be called before any NumPy functions are used?  The answer it an emphatic "YES".  What you need to understand is that the function definition is *not executed* when it is defined, nor can it be as it has no input ``x`` data to process.  That part of the code is just a definition.  The first time the code for the ``sinc(x)`` function is actually executed is when it is called on line 14 of the program, which occurs after the NumPy library is imported in line 10.  The figure below shows the plot of the :math:`\mathrm{sinc}\,x` function generated by the above code.

.. figure:: /chap7/sinc.*
   :scale: 80 %
   :align: center
   :alt: sinc function
   
   Plot of user-defined ``sinc(x)`` function.

.. index::
   single: functions; fast array processing
   single: conditionals; applied to arrays

Fast array processing in user-defined functions
-----------------------------------------------

While using loops to process arrays works just fine, it is usually not the best way to accomplish the task in Python.  The reason is that loops in Python are executed rather slowly.  To deal with this problem, the developers of NumPy introduced a number of functions designed to process arrays quickly and efficiently.  For the present case, what we need is a conditional statement or function that can process arrays directly.  The function we want is called ``where`` and it is a part of the NumPy library.  There ``where`` function has the form

.. sourcecode:: python

    where(condition, output if True, output if False)
    
The first argument of the ``where`` function is a conditional statement involving an array.  The ``where`` function applies the condition to the array element by element, and returns the second argument for those array elements for which the condition is ``True``, and returns the third argument for those array elements that are ``False``.  We can apply it to the ``sinc(x)`` function as follows

.. sourcecode:: python
    
    def sinc(x):
        z = np.where(x==0.0, 1.0, np.sin(x)/x)
        return z

The ``where`` function creates the array ``y`` and sets the elements of ``y`` equal to 1.0 where the corresponding elements of ``x`` are zero, and otherwise sets the corresponding elements to ``sin(x)/x``.  This code executes much faster, 25 to 100 times, depending on the size of the array, than the code using a ``for`` loop.  Moreover, the new code is much simpler to write and read.  An additional benefit of the ``where`` function is that it can handle single variables and arrays equally well.  The code we wrote for the sinc function with the ``for`` loop cannot handle single variables.  Of course we could rewrite the code so that it did, but the code becomes even more clunky.  It's better just to use NumPy's ``where`` function.

The moral of the story
^^^^^^^^^^^^^^^^^^^^^^

The moral of the story is that you should avoid using ``for`` and ``while`` loops to process arrays in Python programs whenever an array-processing method is available.  As a beginning Python programmer, you may not always see how to avoid loops, and indeed, avoiding them is not always possible, but you should look for ways to avoid loops, especially loops that iterate a large number of times. As you become more experienced, you will find that using array-processing methods in Python becomes more natural. Using them can greatly speed up the execution of your code, especially when working with large arrays.

.. index::
   single: functions; multiple inputs and/or outputs

Functions with more (or less) than one input or output
------------------------------------------------------

Python functions can have any number of input arguments and can return any number of variables.  For example, suppose you want a function that outputs :math:`n` :math:`(x,y)` coordinates around a circle of radius :math:`r` centered at the point :math:`(x_0,y_0)`.  The inputs to the function would be  :math:`r`, :math:`x_0`, :math:`y_0`, and  :math:`n`.  The outputs would be the :math:`n` :math:`(x,y)` coordinates.  The following code implements this function.

.. sourcecode:: python

    def circle(r, x0, y0, n):
        theta = np.linspace(0., 2.*np.pi, n, endpoint=False)
        x = r * np.cos(theta)
        y = r * np.sin(theta)
        return x0+x, y0+y

This function has four inputs and two outputs.  In this case, the four inputs are simple numeric variables and the two outputs are NumPy arrays.  In general, the inputs and outputs can be any combination of data types: arrays, lists, strings, *etc*. Of course, the body of the function must be written to be consistent with the prescribed data types.

Functions can also return nothing to the calling program but just perform some task.  For example, here is a program that clears the terminal screen

.. sourcecode:: python

    import subprocess
    import platform
    
    def clear():
        subprocess.Popen( "cls" if platform.system() == 
                          "Windows" else "clear", shell=True)

The function is invoked by typing ``clear()``.  It has no inputs and no outputs but it performs a useful task.  This function uses two standard Python libraries, ``subprocess`` and ``platform`` that are useful for performing computer system tasks.  It's not important that you know anything about them at this point.  We simply use them here to demonstrate a useful cross-platform function that has no inputs and returns no values.

.. index::
   single: functions; arguments; keyword
   single: functions; arguments; positional

Positional and keyword arguments
--------------------------------

It is often useful to have function arguments that have some default setting.  This happens when you want an input to a function to have some standard value or setting most of the time, but you would like to reserve the possibility of giving it some value other than the default value.

For example, in the program ``circle`` from the previous section, we might decide that under most circumstances, we want ``n=12`` points around the circle, like the points on a clock face, and we want the circle to be centered at the origin.  In this case, we would rewrite the code to read

.. sourcecode:: python

    def circle(r, x0=0.0, y0=0.0, n=12):
        theta = np.linspace(0., 2.*np.pi, n, endpoint=False)
        x = r * np.cos(theta)
        y = r * np.sin(theta)
        return x0+x, y0+y
    
The default values of the arguments ``x0``, ``y0``, and ``n`` are specified in the argument of the function definition in the ``def`` line.  Arguments whose default values are specified in this manner are called *keyword arguments*, and they can be omitted from the function call if the user is content using those values.  For example, writing ``circle(4)`` is now a perfectly legal way to call the ``circle`` function and it would produce 12 :math:`(x,y)` coordinates centered about the origin :math:`(x,y)=(0,0)`.  On the other hand, if you want the values of  ``x0``, ``y0``, and ``n`` to be something different from the default values, you can specify their values as you would have before.  

If you want to change only some of the keyword arguments, you can do so by using the keywords in the function call.  For example, suppose you are content with have the circle centered on :math:`(x,y)=(0,0)` but you want only 6 points around the circle rather than 12.  Then you would call the ``circle`` function as follows:

.. sourcecode:: python

    circle(2, n=6)

The unspecified keyword arguments keep their default values of zero but the number of points ``n`` around the circle is now 6 instead of the default value of 12.

The normal arguments without keywords are called *positional arguments*; they have to appear *before* any keyword arguments and, when the function is called, must be supplied values in the same order as specified in the function definition.  The keyword arguments, if supplied, can be supplied in any order providing they are supplied with their keywords.  If supplied without their keywords, they too must be supplied in the order they appear in the function definition.  The following function calls to ``circle`` both give the same output.

.. sourcecode:: ipython

    In [13]: circle(3, n=3, y0=4, x0=-2)
    Out[13]: (array([ 1. , -3.5, -3.5]),
              array([ 4.        ,  6.59807621,  1.40192379]))
    
    In [14]: circle(3, -2, 4, 3)    # w/o keywords, arguments
                                    # supplied in order
    Out[14]: (array([ 1. , -3.5, -3.5]), array([ 4.        ,
                      6.59807621,  1.40192379]))

By now you probably have noticed that we used the keyword argument ``endpoint`` in calling ``linspace`` in our definition of the ``circle`` function.  The default value of ``endpoint`` is ``True``, meaning that ``linspace`` includes the endpoint specified in the second argument of ``linspace``.  We set it equal to ``False`` so that the last point was not included.  Do you see why?

.. index::
   single: functions; arguments; variable number
   single: functions; arguments; *args

Variable number of arguments
----------------------------

While it may seem odd, it is sometimes useful to leave the number of arguments unspecified.  A simple example is a function that computes the product of an arbitrary number of numbers:

.. sourcecode:: python

    def product(*args):
        print("args = {}".format(args))
        p = 1
        for num in args:
            p *= num
        return p

.. sourcecode:: ipython

    In [15]: product(11., -2, 3)
    args = (11.0, -2, 3)
    Out[15]: -66.0
    
    In [16]: product(2.31, 7)
    args = (2.31, 7)
    Out[16]: 16.17

The ``print("args...)`` statement in the function definition is not necessary, of course, but is put in to show that the argument ``args`` is a tuple inside the function.  Here it used because one does not know ahead of time how many numbers are to be multiplied together.

The ``*args`` argument is also quite useful in another context: when passing the name of a function as an argument in another function.  In many cases, the function name that is passed may have a number of parameters that must also be passed but aren't known ahead of time.  If this all sounds a bit confusing---functions calling other functions---a concrete example will help you understand.

Suppose we have the following function that numerically computes the value of the derivative of an arbitrary function :math:`f(x)`:

.. sourcecode:: python

    def deriv(f, x, h=1.e-9, *params):
        return (f(x+h, *params)-f(x-h, *params))/(2.*h)

The argument ``*params`` is an optional positional argument.  We begin by demonstrating the use of the function ``deriv`` without using the optional ``*params`` argument.  Suppose we want to compute the derivative of the function :math:`f_0(x)=4x^5`.  First, we define the function

.. sourcecode:: python

    def f0(x):
        return 4.*x**5

Now let's find the derivative of :math:`f_0(x)=4x^5` at :math:`x=3` using the function ``deriv``:

.. sourcecode:: ipython

    In [17]: deriv(f0, 3)
    Out[17]: 1620.0001482502557

The exact result, given by evaluating :math:`f_0^\prime(x)=20x^4` at :math:`x=3` is 1620, so our function to numerically calculate the derivative works pretty well.

Suppose we had defined a more general function :math:`f_1(x)=ax^p` as follows:

.. sourcecode:: python

    def f1(x, a, p):
        return a*x**p
        
Suppose we want to calculate the derivative of this function for a particular set of parameters :math:`a` and :math:`p`.  Now we face a problem, because it might seem that there is no way to pass the parameters :math:`a` and :math:`p` to the ``deriv`` function.  Moreover, this is a generic problem for functions such as ``deriv`` that use a function as an input, because different functions you want to use as inputs generally come with different parameters. Therefore, we would like to write our program ``deriv`` so that it works, irrespective of how many parameters are needed to specify a particular function.

This is what the optional positional argument ``*params`` defined in ``deriv`` is for: to pass parameters of ``f1``, like :math:`a` and :math:`b`, through ``deriv``.  To see how this works, let's set :math:`a` and :math:`b` to be 4 and 5, respectively, the same values we used in the definition of ``f0``, so that we can compare the results:

.. sourcecode:: ipython

    In [16]: deriv(f1, 3, 1.e-9, 4, 5)
    Out[16]: 1620.0001482502557

We get the same answer as before, but this time we have used ``deriv`` with a more general form of the function :math:`f_1(x)=ax^p`.  

The order of the parameters is important.  The function ``deriv`` uses ``x``, the first argument of ``f1``, as its principal argument, and then uses ``a`` and ``p``, in the same order that they are defined in the function ``f1``, to fill in the additional arguments---the parameters---of the function ``f1``.

.. index::
   single: functions; arguments; **kwargs   

Optional arguments must appear after the regular positional and keyword arguments in a function call. The order of the arguments must adhere to the following convention:

.. sourcecode:: python

    def func(pos1, pos2, ..., keywd1, keywd2, ..., *args, **kwargs):

That is, the order of arguments is: positional arguments first, then keyword arguments, then optional positional arguments (``*args``), then optional keyword arguments (``**kwargs``).  Note that to use the ``*params`` argument, we had to explicitly include the keyword argument ``h`` even though we didn't need to change it from its default value.

Python also allows for a variable number of keyword arguments---``**kwargs``---in a function call.  While ``*args`` is a tuple, ``kwargs`` is a dictionary, so that the value of an optional keyword argument is accessed through its dictionary key.

Passing data to and from functions
----------------------------------

Functions are like mini-programs within the larger programs that call them.  Each function has a set of variables with certain names that are to some degree or other isolated from the calling program.  We shall get more specific about just how isolated those variables are below, but before we do, we introduce the concept of a *namespace*.  Each function has its own namespace, which is essentially a mapping of variable names to objects, like numerics, strings, lists, and so forth.  It's a kind of dictionary.  The calling program has its own namespace, distinct from that of any functions it calls.  The distinctiveness of these namespaces plays an important role in how functions work, as we shall see below.

Variables and arrays created entirely within a function
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

An important feature of functions is that variables and arrays created *entirely within* a function cannot be seen by the program that calls the function unless the variable or array is explicitly passed to the calling program in the ``return`` statement.  This is important because it means you can create and manipulate variables and arrays, giving them any name you please, without affecting any variables or arrays outside the function, even if the variables and arrays inside and outside a function share the same name.

To see what how this works, let's rewrite our program to plot the sinc function using the sinc function definition that uses the ``where`` function.

.. sourcecode:: python
    :linenos:

    def sinc(x):
        z = np.where(x==0.0, 1.0, np.sin(x)/x)
        return z
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    x = np.linspace(-10, 10, 256)
    y = sinc(x)
    
    plt.plot(x, y)
    plt.axhline(color="gray", zorder=-1)
    plt.axvline(color="gray", zorder=-1)
    plt.show()

Running this program produces a plot like the plot of sinc shown in the previous section.  Notice that the array variable ``z`` is only defined within the function definition of sinc.  If we run the program from the IPython terminal, it produces the plot, of course.  Then if we ask IPython to print out the arrays, ``x``, ``y``, and ``z``, we get some interesting and informative results, as shown below.

.. sourcecode:: ipython

    In [15]: run sinc3.py
    
    In [16]: x
    Out[16]: array([-10.        ,  -9.99969482,  -9.99938964, ...,
             9.99938964,   9.99969482,  10.        ])
    
    In [17]: y
    Out[17]: array([-0.05440211, -0.05437816, -0.0543542 , ...,
                    -0.0543542 , -0.05437816, -0.05440211])
    
    In [18]: z
    ---------------------------------------------------------
    NameError               Traceback (most recent call last)
    
    NameError: name 'z' is not defined

When we type in ``x`` at the ``In [16]:`` prompt, IPython prints out the array ``x`` (some of the output is suppressed because the array ``x`` has many elements); similarly for ``y``.  But when we type ``z`` at the ``In [18]:`` prompt, IPython returns a ``NameError`` because ``z`` is not defined.  The IPython terminal is working in the same *namespace* as the program.  But the namespace of the sinc function is isolated from the namespace of the program that calls it, and therefore isolated from IPython.  This also means that when the sinc function ends with ``return z``, it doesn't return the name ``z``, but instead assigns the values in the array ``z`` to the array ``y``, as directed by the main program in line 9.

Passing variables and arrays to functions: mutable and immutable objects
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

What happens to a variable or an array passed to a function when the variable or array is *changed* within the function?  It turns out that the answers are different depending on whether the variable passed is a simple numeric variable, string, or tuple, or whether it is an array or list.  The program below illustrates the different ways that Python handles single variables *vs* the way it handles lists and arrays.

.. sourcecode:: python
    :linenos:

    def test(s, v, t, l, a):
        s = "I am doing fine"
        v = np.pi**2
        t = (1.1, 2.9)
        l[-1] = 'end'
        a[0] = 963.2
        return s, v, t, l, a
        
    import numpy as np
    
    s = "How do you do?"
    v = 5.0
    t = (97.5, 82.9, 66.7)
    l = [3.9, 5.7, 7.5, 9.3]
    a = np.array(l)
    
    print('*************')
    print("s = {0:s}".format(s))
    print("v = {0:5.2f}".format(v))
    print("t = {0:s}".format(t))
    print("l = {0:s}".format(l))
    print("a = "),              # comma suppresses line feed
    print(a)
    print('*************')
    print('*call "test"*')
    
    s1, v1, t1, l1, a1 = test(s, v, t, l, a)
    
    print('*************')
    print("s1 = {0:s}".format(s1))
    print("v1 = {0:5.2f}".format(v1))
    print("t1 = {0:s}".format(t1))
    print("l1 = {0:s}".format(l1))
    print("a1 = "),
    print(a1)
    print('*************')
    print("s = {0:s}".format(s))
    print("v = {0:5.2f}".format(v))
    print("t = {0:s}".format(t))
    print("l = {0:s}".format(l))
    print("a = "),              # comma suppresses line feed
    print(a)
    print('*************')

The function ``test`` has five arguments, a string ``s``, a numerical variable ``v``, a tuple ``t``, a list ``l``, and a NumPy array ``a``.  ``test`` modifies each of these arguments and then returns the modified ``s``, ``v``, ``t``, ``l``, ``a``.  Running the program produces the following output.

.. sourcecode:: ipython

    In [17]: run passingVars.py
    *************
    s = How do you do?
    v =  5.00
    t = (97.5, 82.9, 66.7)
    l = [3.9, 5.7, 7.5, 9.3]
    a =  [ 3.9  5.7  7.5  9.3]
    *************
    *call "test"*
    *************
    s1 = I am doing fine
    v1 =  9.87
    t1 = (1.1, 2.9)
    l1 = [3.9, 5.7, 7.5, 'end']
    a1 =  [ 963.2    5.7    7.5    9.3]
    *************
    s = How do you do?
    v =  5.00
    t = (97.5, 82.9, 66.7)
    l = [3.9, 5.7, 7.5, 'end']
    a =  [ 963.2    5.7    7.5    9.3]
    *************

The program prints out three blocks of variables separated by asterisks.  The first block merely verifies that the contents of ``s``, ``v``, ``t``, ``l``, and ``a`` are those assigned in lines 10-13.  Then the function ``test`` is called.  The next block prints the output of the call to the function ``test``, namely the variables  ``s1``, ``v1``, ``t1``, ``l1``, and ``a1``.  The results verify that the function modified the inputs as directed by the ``test`` function.

The third block prints out the variables ``s``, ``v``, ``t``, ``l``, and ``a`` from the calling program *after* the function ``test`` was called.  These variables served as the inputs to the function ``test``.  Examining the output from the third printing block, we see that the values of the string ``s``, the numeric variable ``v``, and the contents of ``t`` are unchanged after the function call.  This is probably what you would expect.  On the other hand, we see that the list ``l`` and the array ``a`` are changed after the function call.  This might surprise you!  But these are important points to remember, so important that we summarize them in two bullet points here:

    * Changes to string, variable, and tuple arguments of a function within the function do not affect their values in the calling program.

    * Changes to values of elements in list and array arguments of a function within the function are reflected in the values of the same list and array elements in the calling function.
    
The point is that simple numerics, strings and tuples are immutable while lists and arrays are mutable.  Because immutable objects can't be changed, changing them within a function creates new objects with the same name inside of the function, but the old immutable objects that were used as arguments in the function call remain unchanged in the calling program.  On the other hand, if elements of mutable objects like those in lists or arrays are changed, then those elements that are changed inside the function are also changed in the calling program.

Methods and attributes
======================

You have already encountered quite a number of functions that are part of either NumPy or Python or Matplotlib.  But there is another way in which Python implements things that act like functions.  To understand what they are, you need to understand that variables, strings, arrays, lists, and other such data structures in Python are not merely the numbers or strings we have defined them to be.  They are *objects*.  In general, an object in Python has associated with it a number of *attributes* and a number of specialized functions called *methods* that act on the object.  How attributes and methods work with objects is best illustrated by example.

Let's start with the NumPy array.  A NumPy array is a Python object and therefore has associated with it a number of attributes and methods. Suppose, for example, we write ``a = random.random(10)``, which creates an array of 10 uniformly distributed random numbers between 0 and 1.  An example of an attribute of an array is the size or number of elements in the array.  An attribute of an object in Python is accessed by typing the object name followed by a period followed by the attribute name.  The code below illustrates how to access two different attributes of an array, it's size and its data type.

.. sourcecode:: ipython

    In [18]: a = random.random(10)
    
    In [19]: a.size
    Out[19]: 10
    
    In [20]: a.dtype
    Out[20]: dtype('float64')

Any object in Python can and in general does have a number of attributes that are accessed in just the way demonstrated above, with a period and the attribute name following the name of the particular object.  In general, attributes involve properties of the object that are stored by Python with the object and require no computation.   Python just looks up the attribute and returns its value.

Objects in Python also have associated with them a number of specialized functions called *methods* that act on the object.  In contrast to attributes, methods generally involve Python performing some kind of computation.  Methods are accessed in a fashion similar to attributes, by appending a period followed the method's name, which is followed by a pair of open-close parentheses, consistent with methods being a kind of function that acts on the object.  Often methods are used with no arguments, as methods by default act on the object whose name they follow.  In some cases. however, methods can take arguments.  Examples of methods for NumPy arrays are sorting, calculating the mean, or standard deviation of the array.  The code below illustrates a few array methods.

.. sourcecode:: ipython

    In [21]: a
    Out[21]: 
    array([ 0.859057  ,  0.27228037,  0.87780026,  0.14341207,
            0.05067356,  0.83490135,  0.54844515,  0.33583966,
            0.31527767,  0.15868803])
    
    In [22]: a.sum()                # sum
    Out[22]: 4.3963751104791005
    
    In [23]: a.mean()               # mean or average
    Out[23]: 0.43963751104791005
    
    In [24]: a.var()                # variance
    Out[24]: 0.090819477333711512
    
    In [25]: a.std()                # standard deviation
    Out[25]: 0.30136270063448711
    
    In [26]: a.sort()               # sort small to large
    
    In [27]: a
    Out[27]: 
    array([ 0.05067356,  0.14341207,  0.15868803,  0.27228037,
            0.31527767,  0.33583966,  0.54844515,  0.83490135,
            0.859057  ,  0.87780026])

    In [28]: a.clip(0.3, 0.8)
    Out[29]: 
    array([ 0.3       ,  0.3       ,  0.3       ,  0.3       ,
            0.31527767,  0.33583966,  0.54844515,  0.8       ,
            0.8       ,  0.8       ])

The ``clip()`` method provides an example of a method that takes an argument, in this case the arguments are the lower and upper values to which array elements are cutoff if their values are outside the range set by these values.

.. index::
   single: curve fitting; linear

.. _linfitfunc:

Example: linear least squares fitting
=====================================

In this section we illustrate how to use functions and methods in the context of modeling experimental data.

In science and engineering we often have some theoretical curve or *fitting function* that we would like to fit to some experimental data.  In general, the fitting function is of the form :math:`f(x; a, b, c, ...)`, where :math:`x` is the independent variable and :math:`a`, :math:`b`, :math:`c`, ... are parameters to be adjusted so that the function :math:`f(x; a, b, c, ...)` best fits the experimental data.  For example, suppose we had some data of the velocity *vs* time for a falling mass.  If the mass falls only a short distance such that its velocity remains well below its terminal velocity, we can ignore air resistance.  In this case, we expect the acceleration to be constant and the velocity to change linearly in time according to the equation

.. math::
   :label: eq:veltime
   
   v(t) = v_{0} - g t \;,

where :math:`g` is the local gravitational acceleration.  We can fit the data graphically, say by plotting it as shown below in Fig. :ref:`4.6<fig:FallingMassDataPlot>` and then drawing a line through the data.  When we draw a straight line through a data, we try to minimize the distance between the points and the line, globally averaged over the whole data set.  

.. _fig:FallingMassDataPlot:

.. figure:: /chap7/VelocityVsTimePlot.*
   :scale: 80 %
   :align: center
   :alt: Velocity *vs* time for falling mass.
   
   Velocity *vs* time for falling mass.

While this can give a reasonable estimate of the best fit to the data, the procedure is rather *ad hoc*.  We would prefer to have a more well-defined analytical method for determining what constitutes a "best fit".  One way to do that is to consider the sum

.. math::
   :label: eq:lsqrsum

   S = \sum_{i}^{n} [y_{i} - f(x_{i}; a, b, c, ...)]^2 \;,

where :math:`y_{i}` and :math:`f(x_{i}; a, b, c, ...)` are the values of the experimental data and the fitting function, respectively, at :math:`x_{i}`, and :math:`S` is the square of their difference summed over all :math:`n` data points.  The quantity :math:`S` is a sort of global measure of how much the the fit :math:`f(x_{i}; a, b, c, ...)` differs from the experimental data :math:`y_{i}`.  

Notice that for a given set of data points :math:`\{x_i, y_i\}`, :math:`S` is a function only of the fitting parameters :math:`a, b, ...`, that is, :math:`S=S(a, b, c, ...)`.  One way of defining a *best* fit, then, is to find the values of the fitting parameters :math:`a, b, ...` that minimize the :math:`S`.

In principle, finding the values of the fitting parameters :math:`a, b, ...` that minimize the :math:`S` is a simple matter.  Just set the partial derivatives of :math:`S` with respect to the fitting parameter equal to zero and solve the resulting system of equations:

.. math::
   :label: eq:sysSzero
   
   \frac{\partial S}{\partial a} = 0 \;, \quad
   \frac{\partial S}{\partial b} = 0 \;, ...
   
Because there are as many equations as there are fitting paramters, we should be able to solve the system of equations and find the values of the fitting parameters that minimize :math:`S`.  Solving those systems of equations is straightforward if the fitting function :math:`f(x; a, b, ...)` is linear in the fitting parameters.  Some examples of fitting functions linear in the fitting parameters are:

.. math::
   :label: eq:fitfuncs
   
   f(x; a, b) &= a + bx \\
   f(x; a, b, c) &= a + bx + cx^2 \\
   f(x; a, b, c) &= a \sin x + b e^x + c e^{-x^2} \;.
   
For fitting functions such as these, taking the partial derivatives with respect to the fitting parameters, as proposed in :eq:`eq:sysSzero`, results in a set of algebraic equations that are linear in the fitting paramters :math:`a, b, ...` Because they are linear, these equations can be solved in a straightforward manner.

For cases in which the fitting function is not linear in the fitting parameters, one can generally still find the values of the fitting parameters that minimize :math:`S` but finding them requires more work, which goes beyond our immediate interests here.

Linear regression
-----------------
    
We start by considering the simplest case, fitting a straight line to a data set, such as the one shown in Fig. :ref:`4.6 <fig:FallingMassDataPlot>` above.  Here the fitting function is :math:`f(x) = a + bx`, which is linear in the fitting parameters :math:`a` and :math:`b`.  For a straight line, the sum in :eq:`eq:lsqrsum` becomes

.. math::
   :label: eq:linreg1

   S(a,b) = \sum_{i} (y_{i} - a - bx_{i})^2 \;.

Finding the best fit in this case corresponds to finding the values of the fitting parameters :math:`a` and :math:`b` for which :math:`S(a,b)` is a minimum.  To find the minimum, we set the derivatives of :math:`S(a,b)` equal to zero:

.. math::
   :label: eq:linreg2

   \frac{\partial S}{\partial a} &= \sum_{i}-2(y_{i}-a-bx_{i}) = 2 \left(na + b\sum_{i}x_{i} - \sum_{i}y_{i} \right) = 0 \\
   \frac{\partial S}{\partial b} &= \sum_{i}-2(y_{i}-a-bx_{i})\,x_{i} = 2 \left(a\sum_{i}x_{i} + b\sum_{i}x_{i}^2 - \sum_{i}x_{i}y_{i} \right) = 0

Dividing both equations by :math:`2n` leads to the equations

.. math::
   :label: eq:ablinreg

   a + b\bar{x} &= \bar{y}\\
   a\bar{x} + b\frac{1}{n}\sum_{i}x_{i}^2  &=  \frac{1}{n}\sum_{i}x_{i}y_{i}

where

.. math::
   :label: eq:linreg3

   \bar{x} &= \frac{1}{n}\sum_{i}x_{i}\\
   \bar{y} &= \frac{1}{n}\sum_{i}y_{i}\;.
   
Solving Eq. :eq:`eq:ablinreg` for the fitting parameters gives

.. math::
   :label: eq:b1
   
   b &= \frac{\sum_{i}x_{i}y_{i} - n\bar{x}\bar{y}} {\sum_{i}x_{i}^2 - n \bar{x}^2}\\
   a &= \bar{y} - b\bar{x} \;.

Noting that :math:`n\bar{y}=\sum_{i}y` and :math:`n\bar{x}=\sum_{i}x`, the results can be written as
    
.. math:: 
   :label: eq:b2
   
   b &= \frac{\sum_{i}(x_{i}- \bar{x})\,y_{i}} {\sum_{i}(x_{i}- \bar{x})\,x_{i}} \\
   a &= \bar{y} - b\bar{x} \;.
   
While Eqs. :eq:`eq:b1` and :eq:`eq:b2` are equivalent analytically, Eq. :eq:`eq:b2` is preferred for numerical calculations because Eq. :eq:`eq:b2` is less sensitive to roundoff errors.  Here is a Python function implementing this algorithm::

  def LineFit(x, y):
      ''' Returns slope and y-intercept of linear fit to (x,y)
      data set'''
      xavg = x.mean()
      slope = (y*(x-xavg)).sum()/(x*(x-xavg)).sum()
      yint = y.mean()-slope*xavg
      return slope, yint

It's hard to imagine a simpler implementation of the linear regression algorithm.

.. index::
   single: curve fitting; linear; with weighting

Linear regression with weighting: :math:`\chi^2`
------------------------------------------------
    
The linear regression routine of the previous section weights all data points equally.  That is fine if the absolute uncertainty is the same for all data points.  In many cases, however, the uncertainty is different for different points in a data set.  In such cases, we would like to weight the data that has smaller uncertainty more heavily than those data that have greater uncertainty.  For this case, there is a standard method of weighting and fitting data that is known as :math:`\chi^2`  (or *chi-squared*) fitting.  In this method we suppose that associated with each :math:`(x_{i},y_{i})` data point is an uncertainty in the value of :math:`y_{i}` of :math:`\pm\sigma_{i}`.  In this case, the "best fit" is defined as the the one with the set of fitting parameters that minimizes the sum

.. math::
   :label: eq:chisq
    
   \chi^2 = \sum_{i} \left(\frac{y_{i} - f(x_{i})} {\sigma_{i}}\right)^2 \;.

Setting the uncertainties :math:`\sigma_{i}=1` for all data points yields the same sum :math:`S` we introduced in the previous section.  In this case, all data points are weighted equally.  However, if :math:`\sigma_{i}` varies from point to point, it is clear that those points with large :math:`\sigma_{i}` contribute less to the sum than those with small :math:`\sigma_{i}`.  Thus, data points with large :math:`\sigma_{i}` are weighted less than those with small :math:`\sigma_{i}`.

To fit data to a straight line, we set :math:`f(x) = a + bx` and write

.. math::
   :label: eq:chisqlin
   
   \chi^2(a,b) = \sum_{i} \left(\frac{y_{i} - a -bx_{i}} {\sigma_{i}}\right)^2 \;.
    
Finding the minimum for :math:`\chi^2(a,b)` follows the same procedure used for finding the minimum of :math:`S(a,b)` in the previous section.  The result is
 
.. math::
   :label: eq:abwchisq
   
   b &= \frac{\sum_{i}(x_{i} - \hat{x})\,y_{i}/\sigma_{i}^2} {\sum_{i}(x_{i} - \hat{x})\,x_{i}/\sigma_{i}^2}\\
   a &= \hat{y} - b\hat{x} \;.

where

.. math::
   :label: eq:xychisq
   
   \hat{x} &= \frac{\sum_{i}x_{i}/\sigma_{i}^2} {\sum_{i}1/\sigma_{i}^2}\\
   \hat{y} &= \frac{\sum_{i}y_{i}/\sigma_{i}^2} {\sum_{i}1/\sigma_{i}^2}\;.
    
For a fit to a straight line, the overall quality of the fit can be measured by the reduced chi-squared parameter

.. math::
   :label: :eq:linreg13
   
   \chi_{r}^2 = \frac{\chi^2}{n-2}
   
where :math:`\chi^2` is given by Eq. :eq:`eq:chisq` evaluated at the optimal values of :math:`a` and :math:`b` given by Eq. :eq:`eq:abwchisq`.  A good fit is characterized by :math:`\chi_{r}^2 \approx 1`.  This makes sense because if the uncertainties :math:`\sigma_{i}` have been properly estimated, then :math:`[y_{i}-f(x_{i})]^2` should on average be roughly equal to :math:`\sigma_{i}^2`, so that the sum in Eq. :eq:`eq:chisq` should consist of :math:`n` terms approximately equal to 1.  Of course, if there were only 2 terms (`n=2`), then :math:`\chi^2` would be zero as the best straight line fit to two points is a perfect fit.  That is essentially why :math:`\chi_{r}^2` is normalized using :math:`n-2` instead of :math:`n`.  If :math:`\chi_{r}^2` is significantly greater than 1, this indicates a poor fit to the fitting function (or an underestimation of the uncertainties :math:`\sigma_{i}`).  If :math:`\chi_{r}^2` is significantly less than 1, then it indicates that the uncertainties were probably overestimated (the fit and fitting function may or may not be good).

.. _fig:LinRegErrBars:

.. figure:: /chap7/VelocityVsTimeFit.*
   :scale: 80 %
   :align: center
   :alt: Least square linear.
   
   Fit using :math:`\chi^2` least squares fitting routine with data weighted by error bars.
    
We can also get estimates of the uncertainties in our determination of the fitting parameters :math:`a` and :math:`b`, although deriving the formulas is a bit more involved that we want to get into here.  Therefore, we just give the results:

.. math::
   :label: eq:absigma

   \sigma_{b}^2 &= \frac{1} {\sum_{i}(x_{i} - \hat{x})\,x_{i}/\sigma_{i}^2}\\
   \sigma_{a}^2 &= \sigma_{b}^2 \frac{\sum_{i}x_{i}^2/\sigma_{i}^2} {\sum_{i}1/\sigma_{i}^2}\;.
   
The estimates of uncertainties in the fitting parameters depend explicitly on :math:`\{\sigma_{i}\}` and will only be meaningful if (*i*) :math:`\chi_{r}^2 \approx 1` and (*ii*) the estimates of the uncertainties :math:`\sigma_{i}` are accurate.
    
You can find more information, including a derivation of Eq. :eq:`eq:absigma`, in *Data Reduction and Error Analysis for the Physical Sciences, 3rd ed* by P. R. Bevington & D. K. Robinson, McGraw-Hill, New York, 2003.

.. index::
   single: anonymous functions
   single: anonymous functions; lambda expressions
   single: lambda expressions

.. _lambda:

Anonymous functions (lambda)
============================

Python provides another way to generate functions called *lambda* expressions.  A lambda expression is a kind of in-line function that can be generated on the fly to accomplish some small task.  You can assign lambda functions a name, but you don't need to; hence, they are often called *anonymous* functions.  A lambda uses the keyword ``lambda`` and has the general form

.. sourcecode:: ipython

    lambda arg1, arg2, ... : output

The arguments ``arg1, arg2, ...`` are inputs to a lambda, just as for a functions, and the output is an expression using the arguments.

While lambda expressions need not be named, we illustrate their use by comparing a conventional Python function definition to a lambda expression to which we give a name.  First, we define a conventional python function

.. sourcecode:: ipython

    In [1]: def f(a, b):
       ...:     return 3*a+b**2

    In [2]: f(2,3)
    Out[2]: 15

Next, we define a lambda that does the same thing

.. sourcecode:: ipython

    In [3]: g = lambda a, b : 3*a+b**2

    In [4]: g(2,3)
    Out[4]: 15

The ``lambda`` defined by ``g`` does the same thing as the function ``f``.  Such ``lambda`` expressions are useful when you need a very short function definition, usually to be used locally only once or a few times.

Sometimes lambda expressions are used in function arguments that call for a function *name*, as opposed to the function itself.  Moreover, in cases where a the function to be integrated is already defined but is a function one independent variable and several parameters, the lambda expression can be a convenient way of fashioning a single variable function.  Don't worry if this doesn't quite make sense to you right now.  You will see examples of lambda expressions used in just this way in the section :ref:`numericalIntegration`.

There are also a number of nifty programming tricks that can be implemented using ``lambda`` expressions, but we will not go into them here.  Look up ``lambdas`` on the web if you are curious about their more exotic uses.

.. raw:: latex

    \newpage

Exercises
=========

1.  Write a function that can return each of the first three spherical Bessel functions :math:`j_n(x)`:

    .. math::
        :label: eq:bessel
    
        j_0(x) &= \frac{\sin x}{x}\\
        j_1(x) &= \frac{\sin x}{x^2} - \frac{\cos x}{x}\\
        j_2(x) &= \left(\frac{3}{x^2}-1\right)\frac{\sin x}{x} - \frac{3\cos x}{x^2}
        
    Your function should take as arguments a NumPy array :math:`x` and the order :math:`n`, and should return an array of the designated order :math:`n` spherical Bessel function.  Take care to make sure that your functions behave properly at :math:`x=0`.
    
    Demonstrate the use of your function by writing a Python routine that plots the three Bessel functions for :math:`0 \le x \le 20`.  Your plot should look like the one below.  Something to think about: You might note that :math:`j_1(x)` can be written in terms of :math:`j_0(x)`, and that :math:`j_2(x)` can be written in terms of :math:`j_1(x)` and :math:`j_0(x)`.  Can you take advantage of this to write a more efficient function for the calculations of :math:`j_1(x)` and :math:`j_2(x)`?

    .. _fig:besselSph:

    .. figure:: /chap7/besselSph.*
       :scale: 80 %
       :align: center
       :alt: Spherical Bessel functions.

#.  (a) Write a function that simulates the rolling of :math:`n` dice.  Use the NumPy function ``random.random_integers(6)``, which generates a random integer between 1 and 6 with equal probability (like rolling fair dice).  The input of your function should be the number of dice thrown each roll and the output should be the sum of the :math:`n` dice.

    (b) "Roll" 2 dice 10,000 times keeping track of all the sums of each set of rolls in a list.  Then use your program to generate a histogram  summarizing the rolls of two dice 10,000 times.  The result should look like the histogram plotted below.  Use the MatPlotLib function ``hist`` (see http://matplotlib.org/api/pyplot_summary.html) and set the number of bins in the histogram equal to the number of different possible outcomes of a roll of your dice.  For example, the sum of two dice can be anything between 2 and 12, which corresponds to 11 possible outcomes.  You should get a histogram that looks like the one below.
    
    (c) "Repeat part (b) using 3 dice and plot the resulting histogram.

    .. _fig:diceRolln:

    .. figure:: /chap7/diceRoll2.*
       :scale: 80 %
       :align: center
       :alt: Rolling dice.

#.  Write a function to draw a circular smiley face with eyes, a nose, and a mouth.  One argument should set the overall size of the face (the circle radius).  Optional arguments should allow the user to specify the :math:`(x,y)` position of the face, whether the face is smiling or frowning, and the color of the lines.  The default should be a smiling blue face centered at :math:`(0,0)`.  Once you write your function, write a program that calls it several times to produce a plot like the one below (creative improvisation is encouraged!).  In producing your plot, you may find the call ``plt.axes().set_aspect(1)`` useful so that circles appear as circles and not ovals.  You should only use MatPlotLib functions introduced in this text.  To create a circle you can create an array of angles that goes from 0 to :math:`2\pi` and then produce the :math:`x` and :math:`y` arrays for your circle by taking the cosine and sine, respectively, of the array.  Hint: You can use the same :math:`(x,y)` arrays to make the smile and frown as you used to make the circle by plotting appropriate slices of those arrays.  You do not need to create new arrays.

    .. _fig:Faces:

    .. figure:: /chap7/smiley.*
       :scale: 80 %
       :align: center
       :alt: Faces.

#.  In the section :ref:`linfitfunc`, we showed that the best fit of a line :math:`y = a + bx` to a set of data :math:`\{(x_i,y_i)\}` is obtained for the values of :math:`a` and :math:`b` given by Eq. :eq:`eq:b2`.  Those formulas were obtained by finding the values of :math:`a` and :math:`b` that minimized the sum in Eq. :eq:`eq:linreg1`. This approach and these formulas are valid when the uncertainties in the data are the same for all data points.  The Python function ``LineFit(x, y)`` in the section :ref:`linfitfunc` implements Eq. :eq:`eq:b2`.
    
    (a) Write a new fitting function ``LineFitWt(x, y)`` that implements the formulas given in Eq. :eq:`eq:xychisq` that minimize the :math:`\chi^2` function give by Eq. :eq:`eq:chisqlin`.  This more general approach is valid when the individual data points have different weightings *or* when they all have the same weighting.  You should also write a function to calculate the reduced chi-squared :math:`\chi_r^2` defined by Eq. :eq:`eq:chisqlin`.
    
    (b) Write a Python program that reads in the data below, plots it, and fits it using the two fitting functions  ``LineFit(x, y)`` and ``LineFitWt(x, y)``.  Your program should plot the data with error bars and with *both* fits with and without weighting, that is from ``LineFit(x, y)`` and ``LineFitWt(x, y, dy)``.  It should also report the results for both fits on the plot, similar to the output of the supplied program above, as well as the values of :math:`\chi_r^2`, the reduce chi-squared value, for both fits.  Explain why weighting the data gives a steeper or less steep slope than the fit without weighting.
    
        ::
        
            Velocity vs time data
            for a falling mass
            time (s)   velocity (m/s)   uncertainty (m/s)
              2.23          139               16
              4.78          123               16
              7.21          115                4
              9.37           96                9
             11.64           62               17
             14.23           54               17
             16.55           10               12
             18.70           -3               15
             21.05          -13               18
             23.21          -55               10

#.  Modify the function ``LineFitWt(x, y)`` you wrote in Exercise 4 above so that in addition to returning the fitting parameters :math:`a` and :math:`b`, it also returns the uncertainties in the fitting parameters :math:`\sigma_a` and :math:`\sigma_b` using the formulas given by Eq. :eq:`eq:absigma`.  Use your new fitting function to find the uncertainties in the fitted slope and :math:`y`-intercept for the data provided with Exercise 4.