## 14.2.1 Example

The following C-extension module, make the cos function from the standard math library available to Python:

```c
/* Example of wrapping cos function from math.h with the Python-C-API. */

#include <Python.h>
#include <math.h>

/* wrapped cosine function */
static PyObject* cos_func(PyObject* self, PyObject* args)
{
    double value;
    double answer;

    /* parse the input, from python float to c double */
    if (!PyArg_ParseTuple(args, "d", &value))
        return NULL;
    /* if the above function returns -1, an appropriate Python exception will
     * have been set, and the function simply returns NULL
     */

    /* call cos from libm */
    answer = cos(value);

    /* construct the output from cos, from c double to python float */
    return Py_BuildValue("f", answer);
}

/* define functions in module */
static PyMethodDef CosMethods[] =
{
    {"cos_func", cos_func, METH_VARARGS, "evaluate the cosine"},
    {NULL, NULL, 0, NULL}
};

/* module initialization */
PyMODINIT_FUNC

initcos_module(void)
{
    (void) Py_InitModule("cos_module", CosMethods);
}
```

As you can see, there is much boilerplate, both to «massage» the arguments and return types into place and for the module initialisation. Although some of this is amortised, as the extension grows, the boilerplate required for each function(s) remains.

The standard python build system distutils supports compiling C-extensions from a setup.py, which is rather convenient:

```python
from distutils.core import setup, Extension

# define the extension module
cos_module = Extension('cos_module', sources=['cos_module.c'])

# run the setup
setup(ext_modules=[cos_module])
```

This can be compiled:

```bash
$ cd advanced/interfacing_with_c/python_c_api
$ ls
```
```python
cos_module.c setup.py

$ python setup.py build_ext --inplace
running build_ext
building 'cos_module' extension
creating build
creating build/temp.linux-x86_64-2.7
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib -lpython2.7 -o /home/esc/anaconda/lib/libcos_module.so
$ ls
build/ cos_module.c cos_module.so setup.py
```

* `build_ext` is to build extension modules

* `--inplace` will output the compiled extension module into the current directory

The file `cos_module.so` contains the compiled extension, which we can now load in the IPython interpreter:

```python
In [1]: import cos_module

In [2]: cos_module?
Type:         module
String Form: <module 'cos_module' from 'cos_module.so'>
File:         /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/python_c_api/cos_module.so
Docstring:    <no docstring>

In [3]: dir(cos_module)
Out[3]: ['__doc__', '__file__', '__name__', '__package__', 'cos_func']

In [4]: cos_module.cos_func(1.0)
Out[4]: 0.5403023058681398

In [5]: cos_module.cos_func(0.0)
Out[5]: 1.0

In [6]: cos_module.cos_func(3.14159265359)
Out[7]: -1.0
```

Now let's see how robust this is:

```python
In [10]: cos_module.cos_func('foo')
```
```
TypeError: Traceback (most recent call last)
<ipython-input-10-11bee483665d> in <module>()
---> 1 cos_module.cos_func('foo')
TypeError: a float is required
```

## 14.2.2 Numpy Support

Analog to the Python-C-API, Numpy, which is itself implemented as a C-extension, comes with the Numpy-C-API. This API can be used to create and manipulate Numpy arrays from C, when writing a custom C-extension.
See also: `:ref:'advanced numpy'_`.

If you do ever need to use the Numpy C-API refer to the documentation about Arrays and Iterators.

The following example shows how to pass Numpy arrays as arguments to functions and how to iterate over Numpy arrays using the (old) Numpy-C-API. It simply takes an array as argument applies the cosine function from the `math.h` and returns a resulting new array.

```c
/* Example of wrapping the cos function from math.h using the Numpy-C-API. */
```
```c
#include <Python.h>
#include <numpy/arrayobject.h>
#include <math.h>

/* wrapped cosine function */
static PyObject* cos_func NPY NPY NPY NPY
{
    PyArrayObject *in_array;
    PyObject *out_array;
    NpyIter *in_iter;
    NpyIter *out_iter;
    NpyIter_IterNextFunc *in_iternext;
    NpyIter_IterNextFunc *out_iternext;

    /* parse single numpy array argument */
    if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &in_array))
        return NULL;

    /* construct the output array, like the input array */
    out_array = PyArray_NewLikeArray(in_array, NPY_ANYORDER, NULL, 0);
    if (out_array == NULL)
        return NULL;

    /* create the iterators */
    in_iter = NpyIter_New(in_array, NPY ITER_READONLY, NPY_KEEPORDER,
                          NPY_NO_CASTING, NULL);
    if (in_iter == NULL)
        goto fail;

    out_iter = NpyIter_New((PyArrayObject *)out_array, NPY_iter_READWRITE,
                           NPY_KEEPORDER, NPY_NO_CASTING, NULL);
    if (out_iter == NULL) {
        NpyIter_Deallocate(in_iter);
        goto fail;
    }

    in_iternext = NpyIter_GetIterNext(in_iter, NULL);
    out_iternext = NpyIter_GetIterNext(out_iter, NULL);
    if (in_iternext == NULL || out_iternext == NULL) {
        NpyIter_Deallocate(in_iter);
        NpyIter_Deallocate(out_iter);
        goto fail;
    }

    double ** in_dataptr = (double **) NpyIter_GetDataPtrArray(in_iter);
    double ** out_dataptr = (double **) NpyIter_GetDataPtrArray(out_iter);

    /* iterate over the arrays */
    do {
        **out_dataptr = cos(**in_dataptr);
    } while(in_iternext(in_iter) && out_iternext(out_iter));

    /* clean up and return the result */
    NpyIter_Deallocate(in_iter);
    NpyIter_Deallocate(out_iter);
    Py_INCREF(out_array);
    return out_array;

    /* in case bad things happen */
    fail:
        Py_XDECREF(out_array);
        return NULL;
}
```
```c
/* define functions in module */
static PyMethodDef CosMethods[] =
{
    {"cos_func np", cos_func np, METH_VARARGS,
     "evaluate the cosine on a numpy array"},
    {NULL, NULL, 0, NULL}
};

/* module initialization */
PyMODINIT_FUNC
initcos_module np(void)
{
    (void) Py_InitModule("cos_module np", CosMethods);
    /* IMPORTANT: this must be called */
    import_array();
}
```

To compile this we can use `distutils` again. However we need to be sure to include the Numpy headers by using `:func:numpy.get_include`.

```python
from distutils.core import setup, Extension
import numpy

# define the extension module
cos_module np = Extension('cos_module np', sources=['cos_module np.c'],
                         include_dirs=[numpy.get_include()])

# run the setup
setup(ext_modules=[cos_module np])
```

To convince ourselves if this does actually works, we run the following test script:

```python
import cos_module np
import numpy as np
import pylab

x = np.arange(0, 2 * np.pi, 0.1)
y = cos_module np.cos_func np(x)
pylab .plot(x, y)
pylab.show()
```

And this should result in the following figure:

<figure><img src="image_213.png" /></figure>

## 14.3 Ctypes

Ctypes is a *foreign function library* for Python. It provides C compatible data types, and allows calling functions in DLLs or shared libraries. It can be used to wrap these libraries in pure Python.

### Advantages

* Part of the Python standard library

* Does not need to be compiled

* Wrapping code entirely in Python
## Disadvantages

* Requires code to be wrapped to be available as a shared library (roughly speaking *.dll in Windows *.so in Linux and *.dylib in Mac OSX.)

* No good support for C++

### 14.3.1 Example

As advertised, the wrapper code is in pure Python.

```python
""" Example of wrapping cos function from math.h usingctypes. """

importctypes
fromctypes.utilimportfind_library

# find and load the library
libm =ctypes.cdll LoadLibrary (find_library('m'))
# set the argument type
libm .cos .argtypes =[ctypes.c_double]
# set the return type
libm .cos .restype =ctypes.c_double

def cos_func (arg):
    '''Wrapper for cos from math.h '''
    return libm .cos (arg)
```

* Finding and loading the library may vary depending on your operating system, check the documentation for details

* This may be somewhat deceptive, since the math library exists in compiled form on the system already. If you were to wrap a in-house library, you would have to compile it first, which may or may not require some additional effort.

We may now use this, as before:

```
In [1]: import cos_module

In [2]: cos_module?
Type: module
String Form:<module 'cos_module' from 'cos_module.py'>
File: /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/ctypes/cos_module.py
Docstring: <no docstring>

In [3]: dir(cos_module)
Out[3]:
['__builtins__',
'__doc__',
'__file__',
'__name__',
'__package__',
'cos_func',
'ctypes',
'find_library',
'libm']

In [4]: cos_module .cos_func(1.0)
Out[4]: 0.5403023058681398

In [5]: cos_module .cos_func(0.0)
Out[5]: 1.0
```
In [6]: cos_module.cos_func(3.14159265359)
Out[6]: -1.0

As with the previous example, this code is somewhat robust, although the error message is not quite as helpful, since it does not tell us what the type should be.

In [7]: cos_module.cos_func('foo')
------------------------------------------------------------------------
ArgumentError
    Traceback (most recent call last)
  <ipython-input-7-11bee483665d> in <module>()
  ----> 1 cos_module.cos_func('foo')

/home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/ctypes/cos_module.py in cos_func(arg)
  12 def cos_func(arg):
  13     ''' Wrapper for cos from math.h '''
---> 14     return libm.cos(arg)

ArgumentError: argument 1: <type 'exceptions TypeError'>: wrong type

### 14.3.2 Numpy Support

Numpy contains some support for interfacing withctypes. In particular there is support for exporting certain attributes of a Numpy array asctypes data-types and there are functions to convert from C arrays to Numpy arrays and back.

For more information, consult the corresponding section in the Numpy Cookbook and the API documentation for numpy.ndarray.ctypes and numpy.ctypeslib.

For the following example, let's consider a C function in a library that takes an input and an output array, computes the cosine of the input array and stores the result in the output array.

The library consists of the following header file (although this is not strictly needed for this example, we list it for completeness):

void cos_doubles(double *in_array, double *out_array, int size);

The function implementation resides in the following C source file:

# include <math.h>
/* Compute the cosine of each element in in_array, storing the result in
 * out_array. */
void cos_doubles(double *in_array, double *out_array, int size){
    int i;
    for(i=0;i (<size;i++)){
        out_array[i] = cos(in_array[i]);
    }
}

And since the library is pure C, we can't use distutils to compile it, but must use a combination of make and gcc:

m.PHONY : clean

libcos_doubles.so : cos_doubles.o
    gcc -shared -Wl,-soname,libcos_doubles.so -o libcos_doubles.so cos_doubles.o

cos_doubles.o : cos_doubles.c
    gcc -c -fPIC cos_doubles.c -o cos_doubles.o

clean :
    -rm -vf libcos_doubles.so cos_doubles.o cos_doubles.pyc
We can then compile this (on Linux) into the shared library `libcos_doubles.so`:

```bash
$ ls
cos_doubles.c cos_doubles.h cos_doubles.py makefile test_cos_doubles.py
$ make
gcc -c -fPIC cos_doubles.c -o cos_doubles.o
gcc -shared -Wl,-soname,libcos_doubles.so -o libcos_doubles.so cos_doubles.o
$ ls
cos_doubles.c cos_doubles.o libcos_doubles.so* test_cos_doubles.py
cos_doubles.h cos_doubles.py makefile
```

Now we can proceed to wrap this library via `ctypes` with direct support for (certain kinds of) Numpy arrays:

```python
""" Example of wrapping a C library function that accepts a C double array as input using the numpy.ctypeslib. """

import numpy as np
import numpy.ctypeslib as npct
fromctypesimportc_int

# input type for the cos_doubles function
# must be a double array, with single dimension that is contiguous
array_1d_double = npct.ndpointer(dtype = np double, ndim = 1, flags = 'CONTIGUOUS')

# load the library, using numpy mechanisms
libcd = npct.load_library("libcos_doubles", ".")

# setup the return types and argument types
libcd.cos_doubles.re种 = None
libcd.cos_doubles.argtypes = [array_1d_double, array_1d_double, c_int]

def cos_doubles_func(in_array, out_array):
    return libcd.cos_doubles(in_array, out_array, len(in_array))
```

* Note the inherent limitation of contiguous single dimensional Numpy arrays, since the C functions requires this kind of buffer.

* Also note that the output array must be preallocated, for example with `numpy(zeros())` and the function will write into it's buffer.

* Although the original signature of the `cos_doubles` function is `ARRAY`, `ARRAY`, `int` the final `cos_doubles_func` takes only two Numpy arrays as arguments.

And, as before, we convince ourselves that it worked:

```python
import numpy as np
importpylab
importcos_doubles

x = np arange(0, 2 * np pi, 0.1)
y = np empty_like(x)

cos_doubles.cos_doubles_func(x, y)
pylab plot(x, y)
pylab show()
```

<figure><img src="image_214.png" /></figure>
## 14.4 SWIG

SWIG, the Simplified Wrapper Interface Generator, is a software development tool that connects programs written in C and C++ with a variety of high-level programming languages, including Python. The important thing with SWIG is, that it can autogenerate the wrapper code for you. While this is an advantage in terms of development time, it can also be a burden. The generated file tend to be quite large and may not be too human readable and the multiple levels of indirection which are a result of the wrapping process, may be a bit tricky to understand.

The autogenerated C code uses the Python-C-Api.

### Advantages

* Can automatically wrap entire libraries given the headers

* Works nicely with C++

### Disadvantages

* Autogenerates enormous files

* Hard to debug if something goes wrong

* Steep learning curve

### 14.4.1 Example

Let's imagine that our `cos` function lives in a `cos_module` which has been written in C and consists of the source file `cos_module.c`:

```c
#include <math.h>

double cos_func(double arg){
    return cos(arg);
}
```

and the header file `cos_module.h`:

```c
double cos_func(double arg);
```

And our goal is to expose the `cos_func` to Python. To achieve this with SWIG, we must write an *interface file* which contains the instructions for SWIG.

```c
/* Example of wrapping cos function from math.h using SWIG. */

%module cos_module
%{
    /* the resulting C file should be built as a python extension */
    #define SWIG_FILE_WITH_INIT
    /* Includes the header in the wrapper code */
    #include "cos_module.h"
%}

/* Parse the header file to generate wrappers */
%include "cos_module.h"
```

As you can see, not too much code is needed here. For this simple example it is enough to simply include the header file in the interface file, to expose the function to Python. However, SWIG does allow for more fine grained inclusion/exclusion of functions found in header files, check the documentation for details.

Generating the compiled wrappers is a two stage process:

1. Run the `swig` executable on the interface file to generate the files `cos_module-wrap.c`, which is the source file for the autogenerated Python C-extension and `cos_module.py`, which is the autogenerated pure python module.
Python there is no double type):

```python
In [7]: cos_module.cos_func('foo')
```
```
TypeError
```
```
Traceback (most recent call last)
  <ipython-input-7-11bee483665d> in <module>(*)
  ---> 1 cos_module.cos_func('foo')
```
```
TypeError: in method 'cos_func', argument 1 of type 'double'
```

## 14.4.2 Numpy Support

Numpy provides support for SWIG with the ` numpy. i` file. This interface file defines various so-called *typemaps* which support conversion between Numpy arrays and C-Arrays. In the following example we will take a quick look at how such typemaps work in practice.

We have the same `cos_doubles` function as in the `ctypes` example:

```c
void cos_doubles(double * in_array, double * out_array, int size);
```

```c
#include <math.h>

/* Compute the cosine of each element in in_array, storing the result in
 * out_array. */
void cos_doubles(double * in_array, double * out_array, int size){
    int i;
    for(i=0;i (< size; i++){
        out_array[i] = cos(in_array[i]);
    }
}
```

This is wrapped as `cos_doubles_func` using the following SWIG interface file:

```c
/* Example of wrapping a C function that takes a C double array as input using
 * numpy typemaps for SWIG. */

%module cos_doubles
%{
    /* the resulting C file should be built as a python extension */
    #define SWIG_FILE_WITH_INIT
    /* Includes the header in the wrapper code */
    #include "cos_doubles.h"
%}

/* include the numpy typemaps */
%include "numpy.i"
/* need this for correct module initialization */
%init %{
    import_array();
%}

/* typemaps for the two arrays, the second will be modified in-place */
%apply (double* IN_ARRAY1, int DIM1) {(double * in_array, int size_in)}
%apply (double* INPLACE_ARRAY1, int DIM1) {(double * out_array, int size_out)}

/* Wrapper for cos_doubles that masses the types */
%inline %{
    /* takes as input two numpy arrays */
    void cos_doubles_func(double * in_array, int size_in, double * out_array, int size_out) {
        /* calls the original function, providing only the size of the first */
        cos_doubles(in_array, out_array, size_in);
```
```python
}
```

* To use the Numpy typemaps, we need include the ` numpy`. `i` file.

* Observe the call to `import_array()` which we encountered already in the Numpy-C-API example.

* Since the type maps only support the signature `ARRAY`, `SIZE` we need to wrap the `cos_doubles` as `cos_doubles_func` which takes two arrays including sizes as input.

* As opposed to the simple SWIG example, we don't include the `cos_doubles.h` header, There is nothing there that we wish to expose to Python since we expose the functionality through `cos_doubles_func`.

And, as before we can use `distutils` to wrap this:

```python
from distutils.core import setup, Extension
import numpy

setup(ext_modules=[Extension("_cos_doubles",
                             sources=["cos_doubles.c", "cos_doubles.i"],
                             include_dirs=[numpy.get_include()])])
```

As previously, we need to use `include_dirs` to specify the location.

```bash
$ ls
cos_doubles.c cos_doubles.h cos_doubles.i numpy.i setup.py test_cos_doubles.py
$ python setup.py build_ext -i
running build_ext
building '_cos_doubles' extension
swigging cos_doubles.i to cos_doubles-wrap.c
swig -python -o cos_doubles-wrap.c cos_doubles.i
cos_doubles.i:24: Warning(490): Fragment 'NumPy_Backward_Compatibility' not found.
cos_doubles.i:24: Warning(490): Fragment 'NumPy_Backward_Compatibility' not found.
cos_doubles.i:24: Warning(490): Fragment 'NumPy_Backward_Compatibility' not found.
creating build
creating build/temp.linux-x86_64-2.7
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:17:10: warning: 'from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:15' is here
from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:15, 10: warning: 'from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:15' is here
from cos_doublesWrap.c:2706:
/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_dePRECATED_api.h:11:2: warning: #w
gcc -pthread -shared build/temp.linux-x86_64-2.7/cos_doubles.o build/temp.linux-x86_64-2.7/cos_doubles-wrap.o -lpython3.7 -o ls
build/ cos_doubles.h cos_doubles.py cos_doubles-wrap.c setup.py
cos_doubles.c cos_doubles.i _cos_doubles.so* numpy.i test_cos_doubles.py
```

And, as before, we convince ourselves that it worked:

```python
import numpy as np
import pylab
import cos_doubles

x = np.arange(0, 2 * np.pi, 0.1)
y = np.empty_like(x)

cos_doubles.cos_doubles_func(x, y)
pylab.plot(x, y)
pylab.show()
```
<figure><img src="image_215.png" /></figure>

## 14.5 Cython

Cython is both a Python-like language for writing C-extensions and an advanced compiler for this language. The Cython *language* is a superset of Python, which comes with additional constructs that allow you call C functions and annotate variables and class attributes with c types. In this sense one could also call it a *Python with types*.

In addition to the basic use case of wrapping native code, Cython supports an additional use-case, namely interactive optimization. Basically, one starts out with a pure-Python script and incrementally adds Cython types to the bottleneck code to optimize only those code paths that really matter.

In this sense it is quite similar to SWIG, since the code can be autogenerated but in a sense it also quite similar to ctypes since the wrapping code can (almost) be written in Python.

While others solutions that autogenerate code can be quite difficult to debug (for example SWIG) Cython comes with an extension to the GNU debugger that helps debug Python, Cython and C code.

The autogenerated C code uses the Python-C-Api.

### Advantages

* Python like language for writing C-extensions

* Autogenerated code

* Supports incremental optimization

* Includes a GNU debugger extension

* Support for C++ (Since version 0.13)

### Disadvantages

* Must be compiled

* Requires an additional library (but only at build time, at this problem can be overcome by shipping the generated C files)

### 14.5.1 Example

The main Cython code for our `cos_module` is contained in the file `cos_module.pyx`:

```c
""" Example of wrapping cos function from math.h using Cython. """

cdef extern from "math.h":
    double cos(double arg)

def cos_func(arg):
    return cos(arg)
```

Note the additional keywords such as `cdef` and `extern`. Also the `cos_func` is then pure Python.

Again we can use the standard `distutils` module, but this time we need some additional pieces from the `Cython.Distutils`:
```python
from distutils.core import setup, Extension
from Cython.Distutils import build_ext

setup(
    cmdclass={'build_ext': build_ext},
    ext_modules=[Extension("cos_module", ["cos_module.pyx"])]
)
```

### Compiling this:

```bash
$ cd advanced/interfacing_with_c/cython
$ ls
cos_module.pyx      setup.py
$ python setup.py build_ext --inplace
running build_ext
cythoning cos_module.pyx to cos_module.c
building 'cos_module' extension
creating build
creating build/temp.linux-x86_64-2.7
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib -lpython2.7 -o /home/esc/anaconda/lib/cpython3.7-cos_module.so*
$ ls
build/ cos_module.c cos_module.pyx cos_module.so* setup.py
```

### And running it:

```python
In [1]: import cos_module

In [2]: cos_module?
Type:          module
String Form:   <module 'cos_module' from 'cos_module.so'>
File:          /home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/cython/cos_module.so
Docstring:     <no docstring>

In [3]: dir(cos_module)
Out[3]:
['__builtins__',
'__doc__',
'__file__',
'__name__',
'__package__',
'__test__',
'cos_func']

In [4]: cos_module.cos_func(1.0)
Out[4]: 0.5403023058681398

In [5]: cos_module.cos_func(0.0)
Out[5]: 1.0

In [6]: cos_module.cos_func(3.14159265359)
Out[6]: -1.0
```

### And, testing a little for robustness, we can see that we get good error messages:

```python
In [7]: cos_module.cos_func('foo')
```
```
TypeError                                 Traceback (most recent call last)
<ipython-input-7-11bee483665d> in <module>()
---> 1 cos_module.cos_func('foo')
```

```
/home/esc/git-working/scipy-lecture-notes/advanced/interfacing_with_c/cython/cos_module.so in cos_module.cos_func
TypeError: a float is required
```
Additionally, it is worth noting that Cython ships with complete declarations for the C math library, which simplifies the code above to become:

```c
""" Simpler example of wrapping cos function from math.h using Cython. """

from libc.math cimport cos

def cos_func(arg):
    return cos(arg)
```

In this case the `cimport` statement is used to import the `cos` function.

### 14.5.2 Numpy Support

Cython has support for Numpy via the `numpy.pyx` file which allows you to add the Numpy array type to your Cython code. I.e. like specifying that variable `i` is of type `int`, you can specify that variable `a` is of type `numpy.ndarray` with a given `dtype`. Also, certain optimizations such as bounds checking are supported. Look at the corresponding section in the *Cython documentation*. In case you want to pass Numpy arrays as C arrays to your Cython wrapped C functions, there is a section about this in the *Cython wiki*.

In the following example, we will show how to wrap the familiar `cos_doubles` function using Cython.

```c
void cos_doubles(double * in_array, double * out_array, int size);

#include <math.h>

/* Compute the cosine of each element in in_array, storing the result in
 * out_array. */
void cos_doubles(double * in_array, double * out_array, int size){
    int i;
    for(i=0;i (<size; i++)){
        out_array[i] = cos(in_array[i]);
    }
}
```

This is wrapped as `cos_doubles_func` using the following Cython code:

```c
""" Example of wrapping a C function that takes C double arrays as input using
the Numpy declarations from Cython """

# cimport the Cython declarations for numpy
cimport numpy as np

# if you want to use the Numpy-C-API from Cython
# (not strictly necessary for this example, but good practice)
np.import_array()

# cdefine the signature of our c function
cdef extern from "cos_doubles.h":
    void cos_doubles (double * in_array, double * out_array, int size)

# create the wrapper code, with numpy type annotations
def cos_doubles_func(np.ndarray(double, ndim=1, mode="c") in_array not None,
                    np.ndarray(double, ndim=1, mode="c") out_array not None):
    cos_doubles(<double*> np.PyArray_DATA(in_array),
                 <double*> np.PyArray_DATA(out_array),
                 in_array.shape[0])
```

And can be compiled using `distutils`:

```python
from distutils.core import setup, Extension
import numpy
from Cython.Distutils import build_ext
```
```python
setup(
    cmdclass={ 'build_ext': build_ext},
    ext_modules=[Extension("cos_doubles",
                          sources=["_cos_doubles.pyx", "cos_doubles.c"],
                          include_dirs=[numpy.get_include()])],
)
```

* As with the previous compiled Numpy examples, we need the `include_dirs` option.

```bash
$ ls
cos_doubles.c cos_doubles.h _cos_doubles.pyx setup.py test_cos_doubles.py
$ python setup.py build_ext -i
running build_ext
cythoning _cos_doubles.pyx to _cos_doubles.c
building 'cos_doubles' extension
creating build
creating build/temp.linux-x86_64-2.7
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:17:10: warning: _cos_doubles.c:253: In file included from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_de_api.h:11:2: warning: #walled
/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/_ufunc_api.h:236: warning: 'import npy_deungest' is not a function
gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:15:10: warning: _cos_doubles.c:253: In file included from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_deungest_api.h:11:2: warning: #walled
/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/_ufunc_api.h:236: warning: 'import npy_deungest' is not a function
gcc -pthread -shared build/temp.linux-x86_64-2.7/_cos_doubles.o build/temp.linux-x86_64-2.7/cos_doubles.o -L/home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17:10: warning: _cos_doubles.c:253: In file included from /home/esc/anaconda/lib/python2.7/site-packages/numpy/core/include/numpy/npy_deungest_api.h:11:2: warning: #walled
$ ls
build/ _cos_doubles.c cos_doubles.c cos_doubles.h _cos_doubles.pyx cos_doubles.so* setup.py test_cos_doubles.py
```

And, as before, we convince ourselves that it worked:

```python
import numpy as np
import pylab
import cos_doubles

x = np.arange(0, 2 * np.pi, 0.1)
y = np.empty_like(x)

cos_doubles.cos_doubles_func(x, y)
pylab.plot(x, y)
pylab.show()
```

<figure><img src="image_216.png" /></figure>

## 14.6 Summary

In this section four different techniques for interfacing with native code have been presented. The table below roughly summarizes some of the aspects of the techniques.

<table><thead><tr><th>x</th><th>Part of CPython</th><th>Compiled</th><th>Autogenerated</th><th>Numpy Support</th></tr></thead><tbody><tr><td>Python-C-API</td><td>True</td><td>True</td><td>False</td><td>True</td></tr><tr><td>Ctypes</td><td>True</td><td>False</td><td>False</td><td>True</td></tr><tr><td>Swig</td><td>False</td><td>True</td><td>True</td><td>True</td></tr><tr><td>Cython</td><td>False</td><td>True</td><td>True</td><td>True</td></tr></tbody></table>
Of all three presented techniques, Cython is the most modern and advanced. In particular, the ability to optimize code incrementally by adding types to your Python code is unique.

## 14.7 Further Reading and References

* Gaël Varoquaux's blog post about avoiding data copies provides some insight on how to handle memory management cleverly. If you ever run into issues with large datasets, this is a reference to come back to for some inspiration.

## 14.8 Exercises

Since this is a brand new section, the exercises are considered more as pointers as to what to look at next, so pick the ones that you find more interesting. If you have good ideas for exercises, please let us know!

1. Download the source code for each example and compile and run them on your machine.

2. Make trivial changes to each example and convince yourself that this works. (E.g. change `cos` for `sin`.)

3. Most of the examples, especially the ones involving Numpy may still be fragile and respond badly to input errors. Look for ways to crash the examples, figure what the problem is and devise a potential solution. Here are some ideas:

(a) Numerical overflow.

(b) Input and output arrays that have different lengths.

(c) Multidimensional array.

(d) Empty array

(e) Arrays with non-double types

4. Use the `%timeit` IPython magic to measure the execution time of the various solutions

### 14.8.1 Python-C-API

1. Modify the Numpy example such that the function takes two input arguments, where the second is the preallocated output array, making it similar to the other Numpy examples.

2. Modify the example such that the function only takes a single input array and modifies this in place.

3. Try to fix the example to use the new Numpy iterator protocol. If you manage to obtain a working solution, please submit a pull-request on github.

4. You may have noticed, that the Numpy-C-API example is the only Numpy example that does not wrap `cos_doubles` but instead applies the `cos` function directly to the elements of the Numpy array. Does this have any advantages over the other techniques.

5. Can you wrap `cos_doubles` using only the Numpy-C-API. You may need to ensure that the arrays have the correct type, are one dimensional and contiguous in memory.

### 14.8.2 Ctypes

1. Modify the Numpy example such that `cos_doubles_func` handles the preallocation for you, thus making it more like the Numpy-C-API example.
### 14.8.3 SWIG

1. Look at the code that SWIG autogenerates, how much of it do you understand?

2. Modify the Numpy example such that `cos_doubles_func` handles the preallocation for you, thus making it more like the Numpy-C-API example.

3. Modify the `cos_doubles` C function so that it returns an allocated array. Can you wrap this using SWIG typemaps? If not, why not? Is there a workaround for this specific situation? (Hint: you know the size of the output array, so it may be possible to construct a Numpy array from the returned `double *`.)

### 14.8.4 Cython

1. Look at the code that Cython autogenerates. Take a closer look at some of the comments that Cython inserts. What do you see?

2. Look at the section *Working with Numpy* from the Cython documentation to learn how to incrementally optimize a pure python script that uses Numpy.

3. Modify the Numpy example such that `cos_doubles_func` handles the preallocation for you, thus making it more like the Numpy-C-API example.
Part III

**Packages and applications**
This part of the *Scipy lecture notes* is dedicated to various scientific packages useful for extended needs.
# Statistics in Python

**Author:** Gaël Varoquaux

## Requirements

* Standard scientific Python environment (numpy, scipy, matplotlib)

* Pandas

* Statsmodels

* Seaborn

To install Python and these dependencies, we recommend that you download Anaconda Python or Enthought Canopy, or preferably use the package manager if you are under Ubuntu or other linux.

## See also:

### Bayesian statistics in Python

This chapter does not cover tools for Bayesian statistics. Of particular interest for Bayesian modelling is PyMC, which implements a probabilistic programming language in Python.

### Why Python for statistics?

R is a language dedicated to statistics. Python is a general-purpose language with statistics modules. R has more statistical analysis features than Python, and specialized syntaxes. However, when it comes to building complex analysis pipelines that mix statistics with e.g. image analysis, text mining, or control of a physical experiment, the richness of Python is an invaluable asset.
## Contents

* Data representation and interaction

- Data as a table

- The `panda` data-frame

* Hypothesis testing: comparing two groups

- Student's t-test: the simplest statistical test

- Paired tests: repeated measurements on the same individuals

* Linear models, multiple factors, and analysis of variance

- “formulas” to specify statistical models in Python

- Multiple Regression: including multiple factors

- Post-hoc hypothesis testing: analysis of variance (ANOVA)

* More visualization: `seaborn` for statistical exploration

- `Pairplot`: scatter matrices

- `Implot`: plotting a univariate regression

* Testing for interactions

In this document, the Python inputs are represented with the sign `>>>`.

## Disclaimer: Gender questions

Some of the examples of this tutorial are chosen around gender questions. The reason is that on such questions controlling the truth of a claim actually matters to many people.

# 15.1 Data representation and interaction

## 15.1.1 Data as a table

The setting that we consider for statistical analysis is that of multiple *observations* or *samples* described by a set of different *attributes* or *features*. The data can than be seen as a 2D table, or matrix, with columns giving the different attributes of the data, and rows the observations. For instance, the data contained in `examples/brain_size.csv`:

```
"";"Gender";"FSIQ";"VIQ";"PIQ";"Weight";"Height";"MRI_Count"
"1";"Female";133;132;124;"118";"64.5";816932
"2";"Male";140;150;124;".";"72.5";1001121
"3";"Male";139;123;150;"143";"73.3";1038437
"4";"Male";133;129;128;"172";"68.8";965353
"5";"Female";137;132;134;"147";"65.0";951545
```

## 15.1.2 The `panda` data-frame

We will store and manipulate this data in a `pandas Framed`, from the `pandas` module. It is the Python equivalent of the spreadsheet table. It is different from a 2D `numpy` array as it has named columns, can contain a mixture of different data types by column, and has elaborate selection and pivotal mechanisms.

### Creating dataframes: reading data files or converting arrays

#### Separator

It is a CSV file, but the separator is ";"
**Reading from a CSV file:** Using the above CSV file that gives observations of brain size and weight and IQ (Willerman et al. 1991), the data are a mixture of numerical and categorical values:

```python
>>> import pandas
>>> data = pandas.read_csv('examples/brain_size.csv', sep=' ', na_values=".")
>>> data
  Unnamed: 0  Gender  FSIQ  VIQ  PIQ  Weight  Height  MRI_Count
0      1  Female  133  132  124   118   64.5   816932
1      2  Male   140  150  124   NaN   72.5  1001121
2      3  Male   139  123  150   143   73.3  1038437
3      4  Male   133  129  128   172   68.8  965353
4      5  Female  137  132  134   147   65.0  951545
...
```

<figure><img src="image_217.png" /></figure>

### Missing values

The weight of the second individual is missing in the CSV file. If we don't specify the missing value (NA = not available) marker, we will not be able to do statistical analysis.

**Creating from arrays:** A `pandas Frame` can also be seen as a dictionary of 1D 'series', eg arrays or lists. If we have 3 numpy arrays:

```python
>>> import numpy as np
>>> t = np.linspace(-6, 6, 20)
>>> sin_t = np.sin(t)
>>> cos_t = np.cos(t)
```

We can expose them as a `pandas Frame`:

```python
>>> pandas Frame({'t': t, 'sin': sin_t, 'cos': cos_t})
  cos      sin      t
0  0.960170  0.279415 -6.000000
1  0.609977  0.792419 -5.368421
2  0.024451  0.999701 -4.736842
3  -0.570509  0.821291 -4.105263
4  -0.945363  0.326021 -3.473684
5  -0.955488 -0.295030 -2.842105
6  -0.596979 -0.802257 -2.210526
7  -0.008151 -0.999967 -1.578947
8  0.583822 -0.811882 -0.947368
...
```

**Other inputs:** `pandas` can input data from SQL, excel files, or other formats. See the `pandas documentation`.

### Manipulating data

data is a `pandas Frame`, that resembles R's dataframe:

```python
>>> data.shape  # 40 rows and 8 columns
(40, 8)
```
```python
>>> data.columns  # It has columns
Index([u'Unnamed: 0', u'Gender', u'FSIQ', u'VIQ', u'PIQ', u'Weight', u'Height', u'MRI_Count'], dtype='object')
>>> print(data['Gender'])  # Columns can be addressed by name
0      Female
1      Male
2      Male
3      Male
4      Female
...
>>> # Simpler selector
>>> data(data['Gender'] == 'Female')['VIQ'].mean()
109.45
```

For a quick view on a large dataframe, use its `describe()` method: `pandas Frame describe()`.

**groupby**: splitting a dataframe on values of categorical variables:

```python
>>> groupbyGender = data.groupby('Gender')
>>> for gender, value in groupbyGender['VIQ']:
...     print((gender, value.mean()))
('Female', 109.45)
('Male', 115.25)
```

`groupbyGender` is a powerful object that exposes many operations on the resulting group of dataframes:

```python
>>> groupbyGender.mean()
             Unnamed: 0   FSIQ    VIQ    PIQ    Weight   Height   MRI_Count
Gender
Female     19.65   111.9   109.45  110.45  137.200000  65.765000  862654.6
Male       21.35   115.0   115.25  111.60  166.444444  71.431579  954855.4
```

Use tab-completion on `groupbyGender` to find more. Other common grouping functions are median, count (useful for checking to see the amount of missing values in different subsets) or sum. Groupby evaluation is lazy, no work is done until an aggregation function is applied.

<figure><img src="image_218.png" /></figure>
## Exercise

* What is the mean value for VIQ for the full population?

* How many males/females were included in this study?

**Hint** use 'tab completion' to find out the methods that can be called, instead of 'mean' in the above example.

* What is the average value of MRI counts expressed in log units, for males and females?

`groupby gender.boxplot` is used for the plots above (see *this example*).

## Plotting data

Pandas comes with some plotting tools (`pandas.tools plotting`, using ` matplotlib` behind the scene) to display statistics of the data in `dataframes`:

### Scatter matrices:

```python
>>> from pandas.tools import plotting
>>> plotting.scatter_matrix(data[["Weight", "Height", "MRI_Count"]])
```

<figure><img src="image_219.png" /></figure>

```python
>>> plotting.scatter_matrix(data[["PIQ", "VIQ", "FSIQ"]])
```

### Two populations

The IQ metrics are bimodal, as if there are 2 sub-populations.
<figure><img src="image_220.png" /></figure>

### Exercise

Plot the scatter matrix for males only, and for females only. Do you think that the 2 sub-populations correspond to gender?

## 15.2 Hypothesis testing: comparing two groups

For simple statistical tests, we will use the `scipy.stats` sub-module of `scipy`:

```python
>>> from scipy import stats
```

### See also:

Scipy is a vast library. For a quick summary to the whole library, see the *scipy* chapter.
### 15.2.1 Student's t-test: the simplest statistical test

#### 1-sample t-test: testing the value of a population mean

<figure><img src="image_221.png" /></figure>

`scipy.stats.ttest_1samp()` tests if the population mean of data is likely to be equal to a given value (technically if observations are drawn from a Gaussian distributions of given population mean). It returns the T statistic, and the p-value (see the function's help):

```python
>>> stats.ttest_1samp(data['VIQ'], 0)
(...30.088099970..., 1.32891964...e-28)
```

With a p-value of $10^{-28}$ we can claim that the population mean for the IQ (VIQ measure) is not 0.

#### 2-sample t-test: testing for difference across populations

We have seen above that the mean VIQ in the male and female populations were different. To test if this is significant, we do a 2-sample t-test with `scipy.stats.ttest_ind()`:

```python
>>> female_viq = data(data['Gender'] == 'Female']['VIQ']
>>> male_viq = data(data['Gender'] == 'Male']['VIQ']
>>> stats.ttest_ind(female_viq, male_viq)
(...-0.77261617232..., 0.4445287677858...)
```
### 15.2.2 Paired tests: repeated measurements on the same individuals

<figure><img src="image_222.png" /></figure>

PIQ, VIQ, and FSIQ give 3 measures of IQ. Let us test if FISQ and PIQ are significantly different. We can use a 2 sample test:

```python
>>> stats.ttest_ind(data['FSIQ'], data['PIQ'])
(...0.46563759638..., 0.64277250...)
```

The problem with this approach is that it forgets that there are links between observations: FSIQ and PIQ are measured on the same individuals. Thus the variance due to inter-subject variability is confounding, and can be removed, using a "paired test", or "repeated measures test":

```python
>>> stats.ttest_rel(data['FSIQ'], data['PIQ'])
(...1.784201940..., 0.082172638183...)
```

<figure><img src="image_223.png" /></figure>

This is equivalent to a 1-sample test on the difference:

```python
>>> stats.ttest_1samp(data['FSIQ'] - data['PIQ'], 0)
(...1.784201940..., 0.082172638...)
```

T-tests assume Gaussian errors. We can use a Wilcoxon signed-rank test, that relaxes this assumption:

```python
>>> stats.wilcoxon(data['FSIQ'], data['PIQ'])
(274.5, 0.106594927...)
```

The corresponding test in the non paired case is the Mann-Whitney U test, `scipy.stats.mannwhitneyu()`.

#### Exercise

* Test the difference between weights in males and females.

* Use non parametric statistics to test the difference between VIQ in males and females.

**Conclusion:** we find that the data does not support the hypothesis that males and females have different VIQ.
## 15.3 Linear models, multiple factors, and analysis of variance

### 15.3.1 "formulas" to specify statistical models in Python

#### A simple linear regression

<figure><img src="image_224.png" /></figure>

Given two set of observations, x and y, we want to test the hypothesis that y is a linear function of x. In other terms:

$$y = x * \text{coef} + \text{intercept} + e$$

where e is observation noise. We will use the `statmodels` module to:

1. Fit a linear model. We will use the simplest strategy, ordinary least squares (OLS).

2. Test that coef is non zero.

First, we generate simulated data according to the model:

```python
>>> import numpy as np
>>> x = np.linspace(-5, 5, 20)
>>> np.random.seed(1)
>>> # normal distributed noise
>>> y = -5 + 3*x + 4 * np.random.normal(size=x.shape)
>>> # Create a data frame containing all the relevant variables
>>> data = pandas.DataFrame({'x': x, 'y': y})
```

#### "formulas" for statistics in Python

See the `statsmodels` documentation

Then we specify an OLS model and fit it:
```python
>>> from statsmodels.formula api import ols
>>> model = ols("y ~ x", data).fit()
```

We can inspect the various statistics derived from the fit:

```python
>>> print(model summary())
```

OLS Regression Results

========================...

<table><tr><td>Dep. Variable:</td><td>y</td><td>R-squared:</td><td>0.804</td></tr><tr><td>Model:</td><td>OLS</td><td>Adj. R-squared:</td><td>0.794</td></tr><tr><td>Method:</td><td>Least Squares</td><td>F-statistic:</td><td>74.03</td></tr><tr><td>Date:</td><td>...</td><td>Prob (F-statistic):</td><td>8.56e-08</td></tr><tr><td>Time:</td><td>...</td><td>Log-Likelihood:</td><td>-57.988</td></tr><tr><td>No. Observations:</td><td>20</td><td>AIC:</td><td>120.0</td></tr><tr><td>Df Residuals:</td><td>18</td><td>BIC:</td><td>122.0</td></tr><tr><td>Df Model:</td><td>1</td><td></td><td></td></tr></table>

========================...

<table><thead><tr><th></th><th>coef</th><th>std err</th><th>t</th><th>P>|t|</th><th colspan="2">[95.0% Conf. Int.]</th></tr></thead><tbody><tr><td>Intercept</td><td>-5.5335</td><td>1.036</td><td>-5.342</td><td>0.000</td><td>-7.710</td><td>-3.357</td></tr><tr><td>x</td><td>2.9369</td><td>0.341</td><td>8.604</td><td>0.000</td><td>2.220</td><td>3.654</td></tr></tbody></table>

========================...

<table><tr><td>Omnibus:</td><td>0.100</td><td>Durbin-Watson:</td><td>2.956</td></tr><tr><td>Prob(Omnibus):</td><td>0.951</td><td>Jarque-Bera (JB):</td><td>0.322</td></tr><tr><td>Skew:</td><td>-0.058</td><td>Prob(JB):</td><td>0.851</td></tr><tr><td>Kurtosis:</td><td>2.390</td><td>Cond. No.</td><td>3.03</td></tr></table>

## Terminology

Statsmodel uses a statistical terminology: the y variable in statsmodel is called 'endogenous' while the x variable is called exogenous. This is discussed in more detail here.

To simplify, y (endogenous) is the value you are trying to predict, while x (exogenous) represents the features you are using to make the prediction.

## Exercise

Retrieve the estimated parameters from the model above. **Hint**: use tab-completion to find the relevent attribute.

## Categorical variables: comparing groups or multiple categories

Let us go back the data on brain size:

```python
>>> data = pandas.read_csv('examples/brain_size.csv', sep='; ', na_values=".")
```

We can write a comparison between IQ of male and female using a linear model:

```python
>>> model = ols("VIQ ~ Gender + 1", data).fit()
```

```python
>>> print(model summary())
```

OLS Regression Results

========================...

<table><tr><td>Dep. Variable:</td><td>VIQ</td><td>R-squared:</td><td>0.015</td></tr><tr><td>Model:</td><td>OLS</td><td>Adj. R-squared:</td><td>-0.010</td></tr><tr><td>Method:</td><td>Least Squares</td><td>F-statistic:</td><td>0.5969</td></tr></table>
Date: ... Prob (F-statistic): 0.445
Time: ... Log-Likelihood: -182.42
No. Observations: 40 AIC: 368.8
Df Residuals: 38 BIC: 372.2
Df Model: 1

========================...

coef      std err      t      P>|t|      [95.0% Conf. Int.]
-----------------------------------------------------------------...

Intercept      109.4500      5.308      20.619      0.000      98.704      120.196
Gender[T.Male]  5.8000       7.507      0.773       0.445      -9.397       20.997

========================...

Omnibus:                 26.188  Durbin-Watson:          1.709
Prob(Omnibus):           0.000  Jarque-Bera (JB):       3.703
Skew:                    0.010  Prob(JB):               0.157
Kurtosis:                1.510  Cond. No.               2.62

========================...

## Tips on specifying model

**Forcing categorical:** the 'Gender' is automatical detected as a categorical variable, and thus each of its different values are treated as different entities.

An integer column can be forced to be treated as categorical using:

```r
>>> model = ols('VIQ ~ C(Gender)', data).fit()
```

**Intercept:** We can remove the intercept using -1 in the formula, or force the use of an intercept using +1.

By default, statsmodel treats a categorical variable with K possible values as K-1 'dummy' boolean variables (the last level being absorbed into the intercept term). This is almost always a good default choice - however, it is possible to specify different encodings for categorical variables (http://statsmodels.sourceforge.net/devel/contrasts.html).
## Link to t-tests between different FSIQ and PIQ

To compare different type of IQ, we need to create a "long-form" table, listing IQs, where the type of IQ is indicated by a categorical variable:

```python
>>> data_fisq = pandas.DataFrame({'iq': data['FSIQ'], 'type': 'fsiq'})
```

```python
>>> data_piq = pandas.DataFrame({'iq': data['PIQ'], 'type': 'piq'})
```

```python
>>> data_long = pandas.concat((data_fisq, data_piq))
```

```python
>>> print(data_long)
```

<table><thead><tr><th></th><th>type</th></tr></thead><tbody><tr><td>0</td><td>133 fsiq</td></tr><tr><td>1</td><td>140 fsiq</td></tr><tr><td>2</td><td>139 fsiq</td></tr><tr><td>...</td><td></td></tr><tr><td>31</td><td>137 piq</td></tr><tr><td>32</td><td>110 piq</td></tr><tr><td>33</td><td>86 piq</td></tr><tr><td>...</td><td></td></tr></tbody></table>

```python
>>> model = ols("iq ~ type", data_long).fit()
```

```python
>>> print(model-summary())
```

OLS Regression Results

...

============================...

<table><thead><tr><th></th><th>coef</th><th>std err</th><th>t</th><th>P>|t|</th><th colspan="2">[95.0% Conf. Int.]</th></tr></thead><tbody><tr><td>Intercept</td><td>113.4500</td><td>3.683</td><td>30.807</td><td>0.000</td><td>106.119</td><td>120.781</td></tr><tr><td>type[T.piq]</td><td>-2.4250</td><td>5.208</td><td>-0.466</td><td>0.643</td><td>-12.793</td><td>7.943</td></tr><tr><td>...</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

We can see that we retrieve the same values for t-test and corresponding p-values for the effect of the type of iq than the previous t-test:

```python
>>> stats.ttest_ind(data['FSIQ'], data['PIQ'])
(...0.46563759638..., 0.64277250...)
```

## 15.3.2 Multiple Regression: including multiple factors

<figure><img src="image_225.png" /></figure>

Consider a linear model explaining a variable z (the dependent variable) with 2 variables x and y:

$$z = xc_1 + yc_2 + i + e$$

Such a model can be seen in 3D as fitting a plane to a cloud of (x, y, z) points.
### Example: the iris data (examples/iris.csv)

Sepal and petal size tend to be related: bigger flowers are bigger! But is there in addition a systematic effect of species?

<figure><img src="image_226.png" /><figcaption>blue: setosa, green: versicolor, red: virginica</figcaption></figure>

```python
>>> data = pandas.read_csv('examples/iris.csv')
```

```python
>>> model = ols('sepal_width ~ name + petal_length', data).fit()
```

```python
>>> print(model summary())
```

OLS Regression Results

```
============================...
```

<table><tr><td>Dep. Variable:</td><td>sepal_width</td><td>R-squared:</td><td>0.478</td></tr><tr><td>Model:</td><td>OLS</td><td>Adj. R-squared:</td><td>0.468</td></tr><tr><td>Method:</td><td>Least Squares</td><td>F-statistic:</td><td>44.63</td></tr><tr><td>Date:</td><td>...</td><td>Prob (F-statistic):</td><td>1.58e-20</td></tr><tr><td>Time:</td><td>...</td><td>Log-Likelihood:</td><td>-38.185</td></tr><tr><td>No. Observations:</td><td>150</td><td>AIC:</td><td>84.37</td></tr><tr><td>Df Residuals:</td><td>146</td><td>BIC:</td><td>96.41</td></tr><tr><td>Df Model:</td><td>3</td><td></td><td></td></tr></table>

```
============================...
```

<table><thead><tr><th></th><th>coef</th><th>std err</th><th>t</th><th>P>|t|</th><th colspan="2">[95.0% Conf. Int.]</th></tr></thead><tbody><tr><td>Intercept</td><td>2.9813</td><td>0.099</td><td>29.989</td><td>0.000</td><td>2.785</td><td>3.178</td></tr><tr><td>name[T.versicolor]</td><td>-1.4821</td><td>0.181</td><td>-8.190</td><td>0.000</td><td>-1.840</td><td>-1.124</td></tr><tr><td>name[T.virginica]</td><td>-1.6635</td><td>0.256</td><td>-6.502</td><td>0.000</td><td>-2.169</td><td>-1.158</td></tr><tr><td>petal_length</td><td>0.2983</td><td>0.061</td><td>4.920</td><td>0.000</td><td>0.178</td><td>0.418</td></tr></tbody></table>

```
============================...
```

<table><tr><td>Omnibus:</td><td>2.868</td><td>Durbin-Watson:</td><td>1.753</td></tr><tr><td>Prob(Omnibus):</td><td>0.238</td><td>Jarque-Bera (JB):</td><td>2.885</td></tr><tr><td>Skew:</td><td>-0.082</td><td>Prob(JB):</td><td>0.236</td></tr><tr><td>Kurtosis:</td><td>3.659</td><td>Cond. No.</td><td>54.0</td></tr></table>
### 15.3.3 Post-hoc hypothesis testing: analysis of variance (ANOVA)

In the above iris example, we wish to test if the petal length is different between versicolor and virginica, after removing the effect of sepal width. This can be formulated as testing the difference between the coefficient associated to versicolor and virginica in the linear model estimated above (it is an Analysis of Variance, ANOVA). For this, we write a vector of 'contrast' on the parameters estimated: we want to test `"name[T.versicolor] - name[T.virginica]"`, with an F-test:

```python
>>> print(model.f_test([0, 1, -1, 0]))
<F test: F=Marray([ [3.24533535] ]), p=[ [0.07369059] ], df_denom=146, df_num=1>
```

Is this difference significant?

#### Exercise

Going back to the brain size + IQ data, test if the VIQ of male and female are different after removing the effect of brain size, height and weight.

## 15.4 More visualization: seaborn for statistical exploration

Seaborn combines simple statistical fits with plotting on pandas dataframes.

Let us consider a data giving wages and many other personal information on 500 individuals (Berndt, ER. The Practice of Econometrics. 1991. NY: Addison-Wesley).

The full code loading and plotting of the wages data is found in *corresponding example*.

```python
>>> print(data
         EDUCATION SOUTH SEX EXPERIENCE UNION WAGE AGE RACE \
0         8         0    1       21      0  0.707570  35    2
1         9         0    1       42      0  0.694605  57    3
2        12         0    0        1      0  0.824126  19    3
3        12         0    0        4      0  0.602060  22    3
...
```

### 15.4.1 Pairplot: scatter matrices

We can easily have an intuition on the interactions between continuous variables using `seaborn pairplot()` to display a scatter matrix:

```python
>>> import seaborn
>>> seaborn pairplot(data, vars=['WAGE', 'AGE', 'EDUCATION'],
...                  kind='reg')
```
<figure><img src="image_227.png" /></figure>

Categorical variables can be plotted as the

hue:

```python
>>> seaborn pairplot(data, vars=['WAGE', 'AGE', 'EDUCATION'],
...                  kind='reg', hue='SEX')
```

<figure><img src="image_228.png" /></figure>

## Look and feel and matplotlib settings

Seaborn changes the default of matplotlib figures to achieve a more “modern”, “excel-like” look. It does that upon import. You can reset the default using:

```python
>>> from matplotlib import plot as plt
>>> plt.rcdefaults()
```

To switch back to seaborn settings, or understand better styling in seaborn, see the relevant section of the seaborn documentation.
### 15.4.2 Implot: plotting a univariate regression

A regression capturing the relation between one variable and another, eg wage and education, can be plotted using `seaborn.lmplot()`:

```python
>>> seaborn.lmplot(y='WAGE', x='EDUCATION', data=data)
```

<figure><img src="image_229.png" /></figure>

**Robust regression**

Given that, in the above plot, there seems to be a couple of data points that are outside of the main cloud to the right, they might be outliers, not representative of the population, but driving the regression.

To compute a regression that is less sensitive to outliers, one must use a *robust model*. This is done in `seaborn` using `robust=True` in the plotting functions, or in `statsmodels` by replacing the use of the OLS by a “Robust Linear Model”, `statsmodels.formula.api.rlm()`.

## 15.5 Testing for interactions

<figure><img src="image_230.png" /></figure>

males than females?

Do wages increase more with education for
The plot above is made of two different fits. We need to formulate a single model that tests for a variance of slope across the population. This is done via an "interaction".

```python
>>> result = sm.ols(formula='wage ~ education + gender + education * gender',
...                 data=data).fit()
>>> print(result summary())
...

| coef | std err | t | P>|t|  [95.0% Conf. Int.] |
|------|---------|---|------|------------------:|
| Intercept | 0.2998 | 0.072 | 4.173 | 0.000 | 0.159 | 0.441 |
| gender[T.male] | 0.2750 | 0.093 | 2.972 | 0.003 | 0.093 | 0.457 |
| education | 0.0415 | 0.005 | 7.647 | 0.000 | 0.031 | 0.052 |
| education:gender[T.male] | -0.0134 | 0.007 | -1.919 | 0.056 | -0.027 | 0.000 |
...
```

Can we conclude that education benefits males more than females?

## Take home messages

* Hypothesis testing and p-value give you the **significance** of an effect / difference

* **Formulas** (with categorical variables) enable you to express rich links in your data

* **Visualizing** your data and simple model fits matters!

* **Conditionning** (adding factors that can explain all or part of the variation) is important modeling aspect that changes the interpretation.
# Sympy : Symbolic Mathematics in Python

Author: Fabian Pedregosa

## Objectives

1. Evaluate expressions with arbitrary precision.

2. Perform algebraic manipulations on symbolic expressions.

3. Perform basic calculus tasks (limits, differentiation and integration) with symbolic expressions.

4. Solve polynomial and transcendental equations.

5. Solve some differential equations.

**What is SymPy?** SymPy is a Python library for symbolic mathematics. It aims to be an alternative to systems such as Mathematica or Maple while keeping the code as simple as possible and easily extensible. SymPy is written entirely in Python and does not require any external libraries.

Sympy documentation and packages for installation can be found on http://www.sympy.org/

## Chapters contents

* First Steps with SymPy

- Using SymPy as a calculator

- Exercises

- Symbols

* Algebraic manipulations

- Expand

- Simplify

* Calculus

- Limits

- Differentiation

- Series expansion

- Integration

- Exercises

* Equation solving

- Exercises

* Linear Algebra

- Matrices

- Differential Equations
# 16.1 First Steps with SymPy

## 16.1.1 Using SymPy as a calculator

SymPy defines three numerical types: Real, Rational and Integer.

The Rational class represents a rational number as a pair of two Integers: the numerator and the denominator, so Rational(1,2) represents 1/2, Rational(5,2) 5/2 and so on:

```python
>>> from sympy import *
>>> a = Rational(1,2)

>>> a
1/2

>>> a*2
1
```

SymPy uses mpmath in the background, which makes it possible to perform computations using arbitrary-precision arithmetic. That way, some special constants, like e, pi, oo (Infinity), are treated as symbols and can be evaluated with arbitrary precision:

```python
>>> pi**2
2
pi

>>> pi.evalf()
3.14159265358979

>>> (pi + exp(1)).evalf()
5.85987448204884
```

as you see, evalf evaluates the expression to a floating-point number.

There is also a class representing mathematical infinity, called oo:

```python
>>> oo > 99999
True
>>> oo + 1
oo
```

## 16.1.2 Exercises

1. Calculate $\sqrt{2}$ with 100 decimals.

2. Calculate $1/2 + 1/3$ in rational arithmetic.

## 16.1.3 Symbols

In contrast to other Computer Algebra Systems, in SymPy you have to declare symbolic variables explicitly:

```python
>>> from sympy import *
>>> x = Symbol('x')
>>> y = Symbol('y')
```

Then you can manipulate them:

```python
>>> x + y + x - y
2*x

>>> (x + y)**2
```
2
(x + y)

Symbols can now be manipulated using some of python operators: +, -, *, ** (arithmetic), &, |, ~, >>, << (boolean).

Printing

Here we use the following setting for printing

>>> sympy.init_printing( use_unicode=False, wrap_line=True)

16.2 Algebraic manipulations

SymPy is capable of performing powerful algebraic manipulations. We'll take a look into some of the most frequently used: expand and simplify.

16.2.1 Expand

Use this to expand an algebraic expression. It will try to denest powers and multiplications:

>>> expand((x + y)**3)
3    2      2    3
x + 3*x *y + 3*x*y + y
>>> 3*x*y**2 + 3*y*x**2 + x**3 + y**3
3    2      2    3
x + 3*x *y + 3*x*y + y

Further options can be given in form on keywords:

>>> expand(x + y, complex=True)
re(x) + re(y) + I*im(x) + I*im(y)
>>> I*im(x) + I*im(y) + re(x) + re(y)
re(x) + re(y) + I*im(x) + I*im(y)
>>> expand(cos(x + y), trig=True)
-sin(x)*sin(y) + cos(x)*cos(y)
>>> cos(x)*cos(y) - sin(x)*sin(y)
-sin(x)*sin(y) + cos(x)*cos(y)

16.2.2 Simplify

Use simplify if you would like to transform an expression into a simpler form:

>>> simplify((x + x*y) / x)
y + 1

Simplification is a somewhat vague term, and more precises alternatives to simplify exists: powsimp (simplification of exponents), trigsimp (for trigonometric expressions), logcombine, radsimp, together.

Exercises

1. Calculate the expanded form of (x + y)⁶.

2. Simplify the trigonometric expression sin(x)/cos(x)
## 16.3 Calculus

### 16.3.1 Limits

Limits are easy to use in SymPy, they follow the syntax `limit(function, variable, point)`, so to compute the limit of $f(x)$ as $x \to 0$, you would issue `limit(f, x, 0):

```python
>>> limit(sin(x)/x, x, 0)
1
```

you can also calculate the limit at infinity:

```python
>>> limit(x, x, oo)
oo
>>> limit(1/x, x, oo)
0
>>> limit(x**x, x, 0)
1
```

### 16.3.2 Differentiation

You can differentiate any SymPy expression using `diff(func, var)`. Examples:

```python
>>> diff(sin(x), x)
cos(x)
>>> diff(sin(2*x), x)
2*cos(2*x)
>>> diff(tan(x), x)
2
tan(x) + 1
```

You can check, that it is correct by:

```python
>>> limit((tan(x+y) - tan(x))/y, y, 0)
2
tan(x) + 1
```

Higher derivatives can be calculated using the `diff(func, var, n)` method:

```python
>>> diff(sin(2*x), x, 1)
2*cos(2*x)
>>> diff(sin(2*x), x, 2)
-4*sin(2*x)
>>> diff(sin(2*x), x, 3)
-8*cos(2*x)
```

### 16.3.3 Series expansion

SymPy also knows how to compute the Taylor series of an expression at a point. Use `series(expr, var)`:

```python
>>> series(cos(x), x)
    2     4
    x     x      / 6\n
1  -  -  +  -  +  0\x /
    2   24
```
```python
>>> series(1/cos(x), x)
    2      4
      x   5*x    / 6\n
1 + -- + ---- + 0\x /
    2   24
```

### Exercises

1. Calculate $\lim_{x \to 0} \sin(x)/x$

2. Calculate the derivative of $\log(x)$ for $x$.

### 16.3.4 Integration

SymPy has support for indefinite and definite integration of transcendental elementary and special functions via `integrate()` facility, which uses powerful extended Risch-Norman algorithm and some heuristics and pattern matching. You can integrate elementary functions:

```python
>>> integrate(6*x**5, x)
6
x
>>> integrate(sin(x), x)
-cos(x)
>>> integrate(log(x), x)
x*log(x) - x
>>> integrate(2*x + sinh(x), x)
2
x + cosh(x)
```

Also special functions are handled easily:

```python
>>> integrate(exp(-x**2)*erf(x), x)
---- 2
\pi *erf(x)
----------
4
```

It is possible to compute definite integral:

```python
>>> integrate(x**3, (x, -1, 1))
0
>>> integrate(sin(x), (x, 0, pi/2))
1
>>> integrate(cos(x), (x, -pi/2, pi/2))
2
```

Also improper integrals are supported as well:

```python
>>> integrate(exp(-x), (x, 0, oo))
1
>>> integrate(exp(-x**2), (x, -oo, oo))
----
\pi
```

### 16.3.5 Exercises

## 16.4 Equation solving

SymPy is able to solve algebraic equations, in one and several variables:
In [7]: solve(x**4 - 1, x)
Out[7]: [1, 1, -1, -1]

As you can see it takes as first argument an expression that is supposed to be equaled to 0. It is able to solve a large part of polynomial equations, and is also capable of solving multiple equations with respect to multiple variables giving a tuple as second argument:

In [8]: solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
Out[8]: {y: 1, x: -3}

It also has (limited) support for trascendental equations:

In [9]: solve(exp(x) + 1, x)
Out[9]: [pi*I]

Another alternative in the case of polynomial equations is `factor`. `factor` returns the polynomial factorized into irreducible terms, and is capable of computing the factorization over various domains:

In [10]: f = x**4 - 3*x**2 + 1
In [11]: factor(f)
Out[11]: (1 + x - x**2)*(1 - x - x**2)
In [12]: factor(f, modulus=5)
Out[12]: (2 + x)**2*(2 - x)**2

SymPy is also able to solve boolean equations, that is, to decide if a certain boolean expression is satisfiable or not. For this, we use the function `satisfiable`:

In [13]: satisfiable(x & y)
Out[13]: {x: True, y: True}

This tells us that `(x & y)` is True whenever `x` and `y` are both True. If an expression cannot be true, i.e. no values of its arguments can make the expression True, it will return False:

In [14]: satisfiable(x & ~x)
Out[14]: False

### 16.4.1 Exercises

1. Solve the system of equations $x + y = 2$, $2 \cdot x + y = 0$

2. Are there boolean values `x`, `y` that make `(~x | y) & (~y | x)` true?

## 16.5 Linear Algebra

### 16.5.1 Matrices

Matrices are created as instances from the `Matrix` class:

>>> from sympy import Matrix
>>> Matrix([[1,0], [0,1]])
[1 0]
[ ]
[0 1]

unlike a NumPy array, you can also put Symbols in it:

>>> x = Symbol('x')
>>> y = Symbol('y')
>>> A = Matrix([[1,x], [y,1]])
```python
>>> A
[1 x]
[   ]
[y 1]

>>> A**2
[x*y + 1    2*x ]
[          ]
[  2*y     x*y + 1]
```

## 16.5.2 Differential Equations

SymPy is capable of solving (some) Ordinary Differential. To solve differential equations, use dsolve. First, create an undefined function by passing cls=Function to the symbols function:

```python
>>> f, g = symbols('f g', cls=Function)
```

f and g are now undefined functions. We can call f(x), and it will represent an unknown function:

```python
>>> f(x)
f(x)

>>> f(x).diff(x, x) + f(x)
2
d
f(x) + --(f(x))
2
dx

>>> dsolve(f(x).diff(x, x) + f(x), f(x))
f(x) = C1*sin(x) + C2*cos(x)
```

Keyword arguments can be given to this function in order to help if find the best possible resolution system. For example, if you know that it is a separable equations, you can use keyword hint='separable' to force dsolve to resolve it as a separable equation:

```python
>>> dsolve(sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x), f(x), hint='separable')
/           ____________\   /           ____________\
|           / C1       |   |           / C1       |
[f(x) = - asin|  /  ------ + 1 | + pi, f(x) = asin|  /  ------ + 1 | + pi,
|           2       |   |           2       |
\\/        cos (x)  |   |           \\/        cos (x)  |

/           ------------\   /           ------------\-
|           / C1       |   |           / C1       |
f(x) = -asin|  /  ------ + 1 |, f(x) = asin|  /  ------ + 1 |
|           2       |   |           2       |
\\/        cos (x)  |   |           \\/        cos (x)  |
```

### Exercises

1. Solve the Bernoulli differential equation

$$x \frac{df(x)}{x} + f(x) - f(x)^2 = 0$$

2. Solve the same equation using hint='Bernoulli'. What do you observe ?
### filters functions (numpy, skimage, scipy)

```python
>>> import numpy as np
>>> check = np(zeros((9, 9)))
>>> check[:2, 1::2] = 1
>>> check[1::2, ::2] = 1
>>> import matplotlib.pyplot as plt
>>> plt.imshow(check, cmap='gray', interpolation='nearest')
```

<figure><img src="image_231.png" /></figure>

## 17.1.1 scikit-image and the SciPy ecosystem

Recent versions of `scikit-image` is packaged in most Scientific Python distributions, such as Anaconda or Enthought Canopy. It is also packaged for Ubuntu/Debian.

```python
>>> import skimage
>>> from skimage import data  # most functions are in subpackages
```

Most `scikit-image` functions take NumPy ndarrays as arguments

```python
>>> camera = data camera()
>>> camera.dtype
dtype('uint8')
>>> camera.shape
(512, 512)
>>> from skimage import restoration
>>> filteredCamera = restoration.denoise_bilateralcamera)
>>> type(filteredCamera)
<type 'numpy.ndarray'>
```

Other Python packages are available for image processing and work with NumPy arrays:

* `scipy.ndimage`: for nd-arrays. Basic filtering, mathematical morphology, regions properties

* Mahotas

Also, powerful image processing libraries have Python bindings:

* OpenCV (computer vision)

* ITK (3D images and registration)
* and many others

(but they are less Pythonic and NumPy friendly, to a variable extent).

### 17.1.2 What's to be found in scikit-image

* Website: http://scikit-image.org/

* Gallery of examples (as in ` matplotlib` or ` scikit-learn`): http://scikit-image.org/docs/stable/auto_examples/

Different kinds of functions, from boilerplate utility functions to high-level recent algorithms.

* Filters: functions transforming images into other images.

- NumPy machinery

- Common filtering algorithms

* Data reduction functions: computation of image histogram, position of local maxima, of corners, etc.

* Other actions: I/O, visualization, etc.

## 17.2 Input/output, data types and colorspaces

I/O: `skimage.io`

```python
>>> from skimage import io
```

Reading from files: `skimage.io.imread()`

```python
>>> import os
>>> filename = os.path.join(skimage.data_dir, 'camera.png')
>>> camera = io.imread(filename)
```

<figure><img src="image_232.png" /></figure>

Works with all data formats supported by the Python
Imaging Library (or any other I/O plugin provided to `imread` with the `plugin` keyword argument).

Also works with URL image paths:

```python
>>> logo = io.imread('http://scikit-image.org/_static/img/logo.png')
```

Saving to files:
```python
>>> io.imsave('localLogo.png', logo)
```

(imsave also uses an external plugin such as PIL)

### 17.2.1 Data types

<figure><img src="image_233.png" /></figure>

<figure><img src="image_234.png" /></figure>

Image ndarrays can be represented either by integers

(signed or unsigned) or floats.

Careful with overflows with integer data types

```python
>>> camera = data Spray()
>>> camera.dtype
dtype('uint8')
>>> camera_multiplier = 3 * camera
```

Different integer sizes are possible: 8-, 16- or 32-bytes, signed or unsigned.

<figure><img src="image_235.png" /></figure>

An important (if questionable) `skimage` convention: float images are supposed to lie in [-1, 1] (in order to have comparable contrast for all float images)

```python
>>> from skimage import img_as_float
>>> camera_float = img_as_float(camera)
>>> camera.max(), camera_float.max()
(255, 1.0)
```

Some image processing routines need to work with float arrays, and may hence output an array with a different type and the data range from the input array

```python
>>> try:
...     from skimage import filters
... except ImportError:
...     from skimage import filter as filters
>>> camera_sobel = filters.sobelcamera)
>>> camera_sobel.max()
0.591502...
```

<figure><img src="image_236.png" /></figure>

In the example above, we use the `filters` submodule of `scikit-image`, that has been renamed from `filter` to `filters` between versions 0.10 and 0.11, in order to avoid a collision with Python's built-in name `filter`.

Utility functions are provided in `skimage` to convert both the `dtype` and the data range, following `skimage`'s conventions: `util.img_as_float`, `util.img_as_ubyte`, etc.

See the user guide for more details.

### 17.2.2 Colorspaces

Color images are of shape (N, M, 3) or (N, M, 4) (when an alpha channel encodes transparency)

```python
>>> lena = data lena()
>>> lena.shape
(512, 512, 3)
```
Routines converting between different colorspaces (RGB, HSV, LAB etc.) are available in `skimage.color`: `colorIFT2HSV`, `color.LAB2RGB`, etc. Check the docstring for the expected `dtype` (and data range) of input images.

### 3D images

Most functions of `skimage` can take 3D images as input arguments. Check the docstring to know if a function can be used on 3D images (for example MRI or CT images).

### Exercise

Open a color image on your disk as a NumPy array.

Find a `skimage` function computing the histogram of an image and plot the histogram of each color channel

Convert the image to grayscale and plot its histogram.

## 17.3 Image preprocessing / enhancement

Goals: denoising, feature (edges) extraction, ...

### 17.3.1 Local filters

Local filters replace the value of pixels by a function of the values of neighboring pixels. The function can be linear or non-linear.

Neighbourhood: square (choose size), disk, or more complicated *structuring element*.

<figure><img src="image_237.png" /></figure>

Example : horizontal Sobel filter

```python
>>> text = data.text()
>>> hsobel_text = filters.hsobel(text)
```

Uses the following linear kernel for computing horizontal gradients:

```
1   2   1
0   0   0
-1  -2  -1
```

<figure><img src="image_238.png" /></figure>
### 17.3.2 Non-local filters

Non-local filters use a large region of the image (or all the image) to transform the value of one pixel:

```python
>>> from skimage import exposure
>>> camera = data Spray()
>>> camera_equalized = exposure equalizeHist (camera)
```

<figure><img src="image_239.png" /></figure>

<figure><img src="image_240.png" /></figure>

Enhances contrast in large almost uniform regions.

### 17.3.3 Mathematical morphology

See wikipedia for an introduction on mathematical morphology.

Probe an image with a simple shape (a **structuring element**), and modify this image according to how the shape locally fits or misses the image.

Default structuring element: 4-connectivity of a pixel

```python
>>> from skimage import morphology
>>> morphology.diamond(1)
array([[0, 1, 0],
      [1, 1, 1],
      [0, 1, 0]], dtype=uint8)
```

<figure><img src="image_241.png" /></figure>

**Erosion** = minimum filter. Replace the value of a pixel by the minimal value covered by the structuring element.:

```python
>>> a = np(zeros((7,7), dtype(np.int))
>>> a[1:6, 2:5] = 1
>>> a
array([[0, 0, 0, 0, 0, 0],
      [0, 0, 1, 1, 1, 0],
      [0, 0, 1, 1, 1, 0],
      [0, 0, 1, 1, 1, 0],
      [0, 0, 1, 1, 1, 0],
      [0, 0, 1, 1, 1, 0]]
```
```python
>>> morphology-binary_erosion(a, morphology-diamond(1)).astype(npuint8)
array([0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0],
      [0, 0, 0, 1, 0, 0],
      [0, 0, 0, 1, 0, 0],
      [0, 0, 0, 1, 0, 0],
      [0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0]], dtype=uint8)
>>> #Erosion removes objects smaller than the structure
>>> morphology-binary_erosion(a, morphology-diamond(2)).astype(npuint8)
array([0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0]], dtype=uint8)
```

**Dilation:** maximum filter:

```python
>>> a = np(zeros((5, 5)))
>>> a[2, 2] = 1
>>> a
array([[0., 0., 0., 0., 0.],
      [0., 0., 0., 0., 0.],
      [0., 0., 1., 0., 0.],
      [0., 0., 0., 0., 0.],
      [0., 0., 0., 0., 0.]])
>>> morphology-binary_dilation(a, morphology-diamond(1)).astype(npuint8)
array([[0, 0, 0, 0, 0],
      [0, 1, 0, 0],
      [0, 1, 1, 0],
      [0, 0, 1, 0],
      [0, 0, 0, 0, 0]], dtype=uint8)
```

**Opening:** erosion + dilation:

```python
>>> a = np(zeros((5,5)), dtype(np.int))
>>> a[1:4, 1:4] = 1; a[4, 4] = 1
>>> a
array([[0, 0, 0, 0, 0],
      [0, 1, 1, 1, 0],
      [0, 1, 1, 1, 0],
      [0, 1, 1, 1, 0],
      [0, 0, 0, 0, 1]])
>>> morphology-binary_opening(a, morphology-diamond(1)).astype(npuint8)
array([[0, 0, 0, 0, 0],
      [0, 1, 1, 1, 0],
      [0, 0, 1, 0, 0],
      [0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0]], dtype=uint8)
```

Opening removes small objects and smoothes corners.

**Grayscale mathematical morphology**

Mathematical morphology operations are also available for (non-binary) grayscale images (int or float type). Erosion and dilation correspond to minimum (resp. maximum) filters.

Higher-level mathematical morphology are available: tophat, skeletonization, etc.

**See also:**
Basic mathematical morphology is also implemented in `scipy.ndimage.morphology`. The `scipy.ndimage` implementation works on arbitrary-dimensional arrays.

### Example of filters comparison: image denoising

```python
>>> from skimage.morphology import disk
>>> coins = data.coins()
>>> coins_Zoom = coins[10:80, 300:370]
>>> median_coins = filters-rank-median(coins_Zoom, disk(1))
>>> from skimage import restoration
>>> tv_coins = restoration.denoise_tv_chambolle(coins_Zoom, weight=0.1)
>>> gaussian_coins = filters.gaussian_filter(coins, sigma=2)
```

<figure><img src="image_242.png" /></figure>

## 17.4 Image segmentation

Image segmentation is the attribution of different labels to different regions of the image, for example in order to extract the pixels of an object of interest.

### 17.4.1 Binary segmentation: foreground + background

**Histogram-based method: Otsu thresholding**

The Otsu method is a simple heuristic to find a threshold to separate the foreground from the background.

**Earlier scikit-image versions**

`skimage.filters` is called `skimage.filter` in earlier versions of scikit-image

```python
from skimage import data
from skimage import filters
camera = dataordevice
val = filters-threshold_otsucamera)
mask = camera < val
```
<figure><img src="image_243.png" /></figure>

<figure><img src="image_244.png" /></figure>

<figure><img src="image_245.png" /></figure>

## Labeling connected components of a discrete image

Once you have separated foreground objects, it is use to separate them from each other. For this, we can assign a different integer labels to each one.

Synthetic data:

```python
>>> n = 20
>>> l = 256
>>> im = np(zeros((1, 1)))
>>> points = l * np.random.random((2, n ** 2))
>>> im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
>>> im = filters.gaussian_filter(im, sigma=1 / (4. * n))
>>> blobs = im > im.mean()
```

Label all connected components:

```python
>>> from skimage import measure
>>> all labels = measure.label(blobs)
```

Label only foreground connected components:

```python
>>> blobs labels = measure.label(blobs, background=0)
```

<figure><img src="image_246.png" /></figure>

<figure><img src="image_247.png" /></figure>

<figure><img src="image_248.png" /></figure>

**See also:**

`scipy.ndimage.find_objects()` is useful to return slices on object in an image.

## 17.4.2 Marker based methods

If you have markers inside a set of regions, you can use these to segment the regions.
## Watershed segmentation

The Watershed (`skimage.morphology.watershed()`) is a region-growing approach that fills "basins" in the image

```python
>>> from skimage.morphology import watershed
>>> from skimage feature import peak_local_max
>>> 
>>> # Generate an initial image with two overlapping circles
>>> x, y = np indices((80, 80))
>>> x1, y1, x2, y2 = 28, 28, 44, 52
>>> r1, r2 = 16, 20
>>> mask_circle1 = (x - x1) ** 2 + (y - y1) ** 2 < r1 ** 2
>>> mask_circle2 = (x - x2) ** 2 + (y - y2) ** 2 < r2 ** 2
>>> image = np logical_or (mask_circle1, mask_circle2)
>>> # Now we want to separate the two objects in image
>>> # Generate the markers as local maxima of the distance
>>> # to the background
>>> from scipy import ndimage
>>> distance = ndimage distance_transform_edt (image)
>>> local_maxi = peak_local_max (distance, indices = False, footprint = np ones((3, 3)), labels = image)
>>> markers = morphology .label(local_maxi)
>>> labelsWs = watershed(-distance, markers, mask = image)
```

## Random walker segmentation

The random walker algorithm (`skimage.segmentation.random_walker()`) is similar to the Watershed, but with a more "probabilistic" approach. It is based on the idea of the diffusion of labels in the image:

```python
>>> from skimage import segmentation
>>> # Transform markers image so that 0-valued pixels are to
>>> # be labelled, and -1-valued pixels represent background
>>> markers[-image] = -1
>>> labels_rwlock = segmentation .random_walker(image, markers)
```

<figure><img src="image_249.png" /></figure>

## Postprocessing label images

`skimage` provides several utility functions that can be used on label images (ie images where different discrete values identify different regions). Functions names are often self-explaining: `skimage.segmentation.clear边界()`, `skimage.segmentation.relabel_from_one()`, `skimage.morphology.remove_small_objects()`, etc.
### Exercise

* Load the coins image from the data submodule.

* Separate the coins from the background by testing several segmentation methods: Otsu thresholding, adaptive thresholding, and watershed or random walker segmentation.

* If necessary, use a postprocessing function to improve the coins / background segmentation.

## 17.5 Measuring regions' properties

```python
>>> from skimage import measure
```

Example: compute the size and perimeter of the two segmented regions:

```python
>>> properties = measure.regionprops labels_rwlock
>>> [prop.area for prop in properties]
[770.0, 1168.0]
>>> [prop.perimeter for prop in properties]
[100.91..., 126.81...]
```

### See also:

for some properties, functions are available as well in `scipy.ndimage.measurements` with a different API (a list is returned).

### Exercise (continued)

* Use the binary image of the coins and background from the previous exercise.

* Compute an image of labels for the different coins.

* Compute the size and eccentricity of all coins.

## 17.6 Data visualization and interaction

Meaningful visualizations are useful when testing a given processing pipeline.

Some image processing operations:

```python
>>> coins = data.coins()
>>> mask = coins > filters thresholds otsu(coins)
>>> clean边界 = segmentation.clear边界 mask
```

Visualize binary result:

```python
>>> plt.figure()
< matplotlib .figure .Figure object at 0x... >
>>> plt .imshow( clean边界 ,cmap='gray')
< matplotlib .image .AxesImage object at 0x... >
```

Visualize contour

```python
>>> plt .figure()
< matplotlib .figure .Figure object at 0x... >
>>> plt .imshow( coins, cmap='gray')
< matplotlib .image .AxesImage object at 0x... >
>>> plt .contour( clean边界 ,[0.5])
< matplotlib .contour .QuadContourSet ... >
```
Use `skimage` dedicated utility function:

```python
>>> coins_edges = segmentation.mark_boundaries(coins, cleanBorder.astype(np.int))
```

<figure><img src="image_250.png" /></figure>

<figure><img src="image_251.png" /></figure>

The (ex-

## perimental) scikit-image viewer

`skimage viewer` = ` matplotlib-based canvas for displaying images + experimental Qt-based GUI-toolkit`

```python
>>> from skimage import viewer
>>> new viewer = viewer . ImageViewer(coins)
>>> new viewer.show()
```

Useful for displaying pixel values.

For more interaction, plugins can be added to the viewer:

```python
>>> new viewer = viewer . ImageViewer(coins)
>>> from skimage . viewer . plugins import lineprofile
>>> new viewer += lineprofile . LineProfile()
>>> new viewer.show()
```
Line Profile

<figure><img src="image_252.png" /></figure>

<figure><img src="image_253.png" /></figure>

17.7 Feature extraction for computer vision

Geometric or textural descriptor can be extracted from images in order to

* classify parts of the image (e.g. sky vs. buildings)

• match parts of different images (e.g. for object detection)

* and many other applications of Computer Vision
```python
>>> from skimage import feature
```

Example: detecting corners using Harris detector

```python
from skimage feature import corner_harris, corner_subpix, corner_peaks
from skimage transform import warp, AffineTransform

tform = AffineTransform(scale=(1.3, 1.1), rotation=1, shear=0.7,
                         translation=(210, 50))
image = warp(data.checkerboard(), tformpelverse, outputShape=(350, 350))

coords = cornerPeaks corner_harris(image), min_distance=5)
coords_subpix = corner_subpix(image, coordsm, window_size=13)
```

<figure><img src="image_254.png" /></figure>

(this ex-

ample is taken from the `plotCorner` example in `scikit-image`)

Points of interest such as corners can then be used to match objects in different images, as described in the `plotMatching` example of `scikit-image`.
# Traits: building interactive dialogs

**Author: Didrik Pinte**

The Traits project allows you to simply add validation, initialization, delegation, notification and a graphical user interface to Python object attributes.

In this tutorial we will explore the Traits toolset and learn how to dramatically reduce the amount of boilerplate code you write, do rapid GUI application development, and understand the ideas which underlie other parts of the Enthought Tool Suite.

Traits and the Enthought Tool Suite are open source projects licensed under a BSD-style license.

## Intended Audience

Intermediate to advanced Python programmers

## Requirements

- Either `wxPython`, `PyQt` or `PySide`

- `Numpy` and `Scipy`

- `Enthought Tool Suite`

- All required software can be obtained by installing the EPD Free

## Tutorial content

- Introduction

- Example

- What are Traits

- Initialization

- Validation

- Documentation

- Visualization: opening a dialog

- Deferral

- Notification

- Some more advanced traits
## 18.1 Introduction

The Enthought Tool Suite enable the construction of sophisticated application frameworks for data analysis, 2D plotting and 3D visualization. These powerful, reusable components are released under liberal BSD-style licenses.

<figure><img src="image_255.png" /></figure>

The main packages of the Enthought Tool Suite are:

* Traits - component based approach to build our applications.

* Kiva - 2D primitives supporting path based rendering, affine transforms, alpha blending and more.

* Enable - object based 2D drawing canvas.

* Chaco - plotting toolkit for building complex interactive 2D plots.

* Mayavi - 3D visualization of scientific data based on VTK.

* Envisage - application plugin framework for building scriptable and extensible applications

In this tutorial, we will focus on Traits.

## 18.2 Example

Throughout this tutorial, we will use an example based on a water resource management simple case. We will try to model a dam and reservoir system. The reservoir and the dams do have a set of parameters :

* Name

* Minimal and maximal capacity of the reservoir [hm3]

* Height and length of the dam [m]

* Catchment area [km2]

* Hydraulic head [m]
* Power of the turbines [MW]

* Minimal and maximal release [m3/s]

* Efficiency of the turbines

The reservoir has a known behaviour. One part is related to the energy production based on the water released. A simple formula for approximating electric power production at a hydroelectric plant is $P = \rho hrgk$, where:

* $P$ is Power in watts,

* $\rho$ is the density of water (~1000 kg/m3),

* $h$ is height in meters,

* $r$ is flow rate in cubic meters per second,

* $g$ is acceleration due to gravity of 9.8 m/s2,

* $k$ is a coefficient of efficiency ranging from 0 to 1.

Annual electric energy production depends on the available water supply. In some installations the water flow rate can vary by a factor of 10:1 over the course of a year.

The second part of the behaviour is the state of the storage that depends on controlled and uncontrolled parameters :

$$ storage_{t+1} = storage_{t} + inflows - release - spillage - irrigation $$

<figure><img src="image_256.png" /></figure>

The data used in this tutorial are not real and might even not have sense in the reality.

## 18.3 What are Traits

A trait is a type definition that can be used for normal Python object attributes, giving the attributes some additional characteristics:

* **Standardization:**

- Initialization

- Validation

- Deferral

* Notification

* Visualization

* Documentation

A class can freely mix trait-based attributes with normal Python attributes, or can opt to allow the use of only a fixed or open set of trait attributes within the class. Trait attributes defined by a class are automatically inherited by any subclass derived from the class.

The common way of creating a traits class is by extending from the **HasTraits** base class and defining class traits :

```python
from traits.api import HasTraits, Str, Float
```

```python
class Reservoir(HasTraits):
```

```python
name = Str
```

```python
max_storage = Float
```
<figure><img src="image_257.png" /></figure>

For Traits 3.x users

If using Traits 3.x, you need to adapt the namespace of the traits packages:

* `traits.api` should be `enthought.traits api`

* `traitsuiapi` should be `enthought.traits.uiapi`

Using a `traits` class like that is as simple as any other Python class. Note that the trait value are passed using keyword arguments:

```python
reservoir = Reservoir(name='Lac de Vouglans', max_storage=605)
```

### 18.3.1 Initialisation

All the traits do have a default value that initialise the variables. For example, the basic python types do have the following trait equivalents:

<table><thead><tr><td>Trait</td><td>Python Type</td><td>Built-in Default Value</td></tr></thead><tbody><tr><td>Bool</td><td>Boolean</td><td>False</td></tr><tr><td>Complex</td><td>Complex number</td><td>0+0j</td></tr><tr><td>Float</td><td>Floating point number</td><td>0.0</td></tr><tr><td>Int</td><td>Plain integer</td><td>0</td></tr><tr><td>Long</td><td>Long integer</td><td>0L</td></tr><tr><td>Str</td><td>String</td><td>''</td></tr><tr><td>Unicode</td><td>Unicode</td><td>u''</td></tr></tbody></table>

A number of other predefined trait type do exist: `Array`, `Enum`, `Range`, `Event`, `Dict`, `List`, `Color`, `Set`, `Expression`, `Code`, `Callable`, `Type`, `Tuple`, etc.

Custom default values can be defined in the code:

```python
from traits.api import HasTraits, Str, Float
```

```python
class Reservoir(HasTraits):
```

```python
name = Str
```

```python
max_storage = Float(100)
```

```python
reservoir = Reservoir(name='Lac de Vouglans')
```

#### Complex initialisation

When a complex initialisation is required for a trait, a `_XXX_default` magic method can be implemented. It will be lazily called when trying to access the XXX trait. For example:

```python
def _name_default(self):
    """ Complex initialisation of the reservoir name. """
    return 'Undefined'
```

### 18.3.2 Validation

Every trait does validation when the user tries to set its content:

```python
reservoir = Reservoir(name='Lac de Vouglans', max_storage=605)
```

```python
reservoir.max_storage = '230'
```

```
Traceback (most recent call last)
  .../scipy-lecture-notes/advanced/traits/<ipython-input-7-979bdff9974a> in <module>()
  ...
---> 1 reservoir.max_storage = '230'
```
.../traits/trait_handlers.pyc in error(self, object, name, value)
    166          """
    167          raise TraitError( object, name, self.full_info( object, name, value ),
    168                          value )
    169
    170          def arg_error ( self, method, arg_num, object, name, value ):

TraitError: The 'max_storage' trait of a Reservoir instance must be a float, but a value of '23'

<type 'str'>

### 18.3.3 Documentation

By essence, all the traits do provide documentation about the model itself. The declarative approach to the creation of classes makes it self-descriptive:

```python
from traits.api import HasTraits, Str, Float

class Reservoir(HasTraits):

    name = Str
    max_storage = Float(100)
```

The `desc` metadata of the traits can be used to provide a more descriptive information about the trait:

```python
from traits.api import HasTraits, Str, Float
class Reservoir(HasTraits):

    name = Str
    max_storage = Float(100, desc='Maximal storage [hm3]')
```

Let's now define the complete reservoir class:

```python
from traits.api import HasTraits, Str, Float, Range
class Reservoir(HasTraits):
    name = Str
    max_storage = Float(1e6, desc='Maximal storage [hm3]')
    max_release = Float(10, desc='Maximal release [m3/s]')
    head = Float(10, desc='Hydraulic head [m]')
    efficiency = Range(0, 1.)

    def energy_production(self, release):
        """ Returns the energy production [Wh] for the given release [m3/s] """
        power = 1000 * 9.81 * self-head * release * self-efficiency
        return power * 3600

if __name__ == '__main__':
    reservoir = Reservoir(
        name = 'Project A',
        max_storage = 30,
        max_release = 100.0,
        head = 60,
        efficiency = 0.8
    )

    release = 80
    print 'Releasing {} m3/s produces {} kWh'.format(
        release, reservoir.energy_production(release)
    )
```
### 18.3.4 Visualization: opening a dialog

The Traits library is also aware of user interfaces and can pop up a default view for the Reservoir class:

```python
reservoir1 = Reservoir()
```

```python
reservoir1.edit traits()
```

<figure><img src="image_258.png" /></figure>

TraitsUI simplifies the way user interfaces are created. Every trait on a HasTraits class has a default editor that will manage the way the trait is rendered to the screen (e.g. the Range trait is displayed as a slider, etc.).

In the very same vein as the Traits declarative way of creating classes, TraitsUI provides a declarative interface to build user interfaces code:

```python
from traits.api import HasTraits, Str, Float, Range
```

```python
from traitsuiapidoc import View
```

```python
class Reservoir(HasTraits):
```

```python
    name = Str
    max_storage = Float(1e6, desc='Maximal storage [hm3]')
    max_release = Float(10, desc='Maximal release [m3/s]')
    head = Float(10, desc='Hydraulic head [m]')
    efficiency = Range(0, 1.)
```

```python
    traits_view = View(
        'name', 'max_storage', 'max_release', 'head', 'efficiency',
        title = 'Reservoir',
        resizable = True,
    )
```

```python
def energy_production(self, release):
    ''' Returns the energy production [Wh] for the given release [m3/s] '''
    power = 1000 * 9.81 * self-head * release * self.efficiency
    return power * 3600
```

```python
if __name__ == '__main__':
```

```python
    reservoir = Reservoir(
        name = 'Project A',
        max_storage = 30,
        max_release = 100.0,
        head = 60,
        efficiency = 0.8
    )
```
```python
reservoir.configure()
```

<figure><img src="image_259.png" /></figure>

### 18.3.5 Deferral

Being able to defer the definition of a trait and its value to another object is a powerful feature of Traits.

```python
from traits.api import HasTraits, Instance, DelegatesTo, Float, Range
```

```python
from reservoir import Reservoir
```

```python
class ReservoirState(HasTraits):
    """Keeps track of the reservoir state given the initial storage.
    """
```

```python
reservoir = Instance(Reservoir,())
min_storage = Float
max_storage = DelegatesTo('reservoir')
min_release = Float
max_release = DelegatesTo('reservoir')
```

```python
# state attributes
storage = Range(low='min_storage', high='max_storage')
```

```python
# control attributes
inflows = FloatInflows [hm3]')
release = Range(low='min_release', high='max_release')
spillage = Float desc='Spillage [hm3]')
```

```python
def print_state(self):
    print 'Storage\tRelease\tInflows\tSpillage'
    str_format = '\t'.join([':{7.2f}'for i in range(4)])
    print str_format.format(self.storage, self.release, self.inflows,
                             self.spillage)
    print '-' * 79
```

```python
if __name__ == '__main__':
```

```python
projectA = Reservoir(
    name = 'Project A',
    max_storage = 30,
    max_release = 100.0,
    hydraulic_head = 60,
    efficiency = 0.8
)

state = ReservoirState(reservoir=projectA, storage=10)
state.release = 90
state.inflows = 0
state.print_state()
```
print 'How do we update the current storage ?'

A special trait allows to manage events and trigger function calls using the magic `xxxx_fired` method:

```python
from traits.api import HasTraits, Instance, DelegatesTo, Float, Range, Event

from reservoir import Reservoir

class ReservoirState(HasTraits):
    """Keeps track of the reservoir state given the initial storage.

    For the simplicity of the example, the release is considered in
    hm3/timestep and not in m3/s.
    """
    reservoir = Instance(Reservoir,())
    min_storage = Float
    max_storage = DelegatesTo('reservoir')
    min_release = Float
    max_release = DelegatesTo('reservoir')

    # state attributes
    storage = Range(low='min_storage', high='max_storage')

    # control attributes
    inflows = FloatInflows [hm3]')
    release = Range(low='min_release', high='max_release')
    spillage = Float desc='Spillage [hm3]')

    update_storage = Event desc='Updates the storage to the next time step')

    def _update_storage_fired(self):
        """update storage state
        new_storage = self.storage - self.release + self.inflows
        self.storage = min(new_storage, self.max_storage)
        overflow = new_storage - self.max_storage
        self.spillage = max(overflow, 0)
        """
    def print_state(self):
        print 'Storage \tRelease \tInflows \tSpillage'
        str_format = '\t'.join([':{7.2f}'for i in range(4)])
        print str_format.format(self.storage, self.release, self.inflows,
                                 self.spillage)
        print '-' * 79

if __name__ == '__main__':
    projectA = Reservoir(
        name='Project A',
        max_storage=30,
        max_release=5.0,
        hydraulic_head=60,
        efficiency=0.8
    )

    state = ReservoirState(reservoir=projectA, storage=15)
    state.release = 5
    state.inflows = 0

    # release the maximum amount of water during 3 time steps
    state.update_storage = True
    state.update_storage = True
    state.print_state()
```
```python
state.update_storage = True
state.print_state()
```

Dependency between objects can be made automatic using the trait **Property**. The **depends_on** attribute expresses the dependency between the property and other traits. When the other traits gets changed, the property is invalidated. Again, Traits uses magic method names for the property:

* `_get_XXX` for the getter of the XXX Property trait

* `_set_XXX` for the setter of the XXX Property trait

```python
from traits.api import HasTraits, Instance, DelegatesTo, Float, Range
from traits.api import Property

from reservoir import Reservoir

class ReservoirState(HasTraits):
    """Keeps track of the reservoir state given the initial storage.

    For the simplicity of the example, the release is considered in
    hm3/timestep and not in m3/s.
    """
    reservoir = Instance(Reservoir,())
    max_storage = DelegatesTo('reservoir')
    min_release = Float
    max_release = DelegatesTo('reservoir')

    # state attributes
    storage = Property depends_on='inflows, release')

    # control attributes
    inflows = Float (desc='Inflows [hm3]')
    release = Range(low='min_release', high='max_release')
    spillage = Property(
        desc='Spillage [hm3]', depends_on=['storage', 'inflows', 'release']
    )

    ### Private traits.
    _storage = Float

    ### Traits property implementation.
    def _get_storage(self):
        new_storage = self._storage - self.release + self.inflows
        return min(new_storage, self.max_storage)

    def _set_storage(self, storage_value):
        self._storage = storage_value

    def _get_spillage(self):
        new_storage = self._storage - self.release + self.inflows
        overflow = new_storage - self.max_storage
        return max(overflow, 0)

    def print_state(self):
        print 'Storage\tRelease\tInflows\tSpillage'
        str_format = '\t'.join([':{:.2f}' for i in range(4)])
        print str_format.format(self.storage, self.release, self.inflows,
                                 self.spillage)
        print '-' * 79

if __name__ == '__main__':
    projectA = Reservoir(
        name='Project A',
        max_storage = 30,
```
```python
        max_release = 5,
        hydraulic_head = 60,
        efficiency = 0.8
    )

state = ReservoirState(reservoir=projectA, storage=25)
state.release = 4
state.inflows = 0

state.print_state()
```

### Caching property

Heavy computation or long running computation might be a problem when accessing a property where the inputs have not changed. The `@cached_property` decorator can be used to cache the value and only recompute them once invalidated.

Let's extend the `TraitsUI` introduction with the `ReservoirState` example:

```python
from traits.api import HasTraits, Instance, DelegatesTo, Float, Range, Property
from traitsui-api import View, Item, Group, VGroup

from reservoir import Reservoir

class ReservoirState(HasTraits):
    """Keeps track of the reservoir state given the initial storage.

    For the simplicity of the example, the release is considered in
    hm3/timestep and not in m3/s.
    """
    reservoir = Instance(Reservoir,())
    name = DelegatesTo('reservoir')
    max_storage = DelegatesTo('reservoir')
    max_release = DelegatesTo('reservoir')
    min_release = Float

    # state attributes
    storage = Property depends_on='inflows, release')

    # control attributes
    inflows = Float (desc='Inflows [hm3]')
    release = Range(low='min_release', high='max_release')
    spillage = Property(
        desc='Spillage [hm3]', depends_on=['storage', 'inflows', 'release']
    )

### Traits view
traits_view = View(
    Group(
        VGroup (Item('name'), Item('storage'), Item('spillage'),
            label='State', style='readonly'
        ),
        VGroup (Item('inflows'), Item('release'), label='Control'),
    )
)

### Private traits.
_storage = Float

### Traits property implementation.
def _get_storage(self):
    new_storage = self._storage - self.release + self.inflows
    return min(new_storage, self.max_storage)
```
```python
def _set_storage(self, storage_value):
    self._storage = storage_value
```

```python
def _get_spillage(self):
    new_storage = self._storage - self.release + self.inflows
    overflow = new_storage - self.max_storage
    return max(overflow, 0)
```

```python
def print_state(self):
    print 'Storage\tRelease\tInflows\tSpillage'
    str_format = '\t'.join(['{:.2f}' for i in range(4)])
    print str_format.format(self.storage, self.release, self.inflows,
                             self.spillage)
    print '-' * 79
```

```python
if __name__ == '__main__':
    projectA = Reservoir(
        name = 'Project A',
        max_storage = 30,
        max_release = 5,
        hydraulic_head = 60,
        efficiency = 0.8
    )
```

```python
state = ReservoirState(reservoir=projectA, storage=25)
```

```python
state.release = 4
```

```python
state.inflows = 0
```

```python
state.print_state()
```

```python
state.configure(traits())
```

<figure><img src="image_260.png" /></figure>

Some use cases need the delegation mechanism to be broken by the user when setting the value of the trait. The **PrototypeFrom** trait implements this behaviour.

```python
from traits.api import HasTraits, Str, Float, Range, PrototypedFrom, Instance
```

```python
class Turbine(HasTraits):
    turbine_type = Str
    power = Float(1.0, desc='Maximal power delivered by the turbine [Mw]')
```

```python
class Reservoir(HasTraits):
    name = Str
    max_storage = Float(1e6, desc='Maximal storage [hm3]')
```
```python
max_release = Float(10, desc='Maximal release [m3/s]')
head = Float(10, desc='Hydraulic head [m]')
efficiency = Range(0, 1.)
```

```python
turbine = Instance(Turbine)
installed_capacity = PrototypedFrom('turbine', 'power')
```

```python
if __name__ == '__main__':
    turbine = Turbine(turbine_type='type1', power=5.0)
```

```python
reservoir = Reservoir(
    name = 'Project A',
    max_storage = 30,
    max_release = 100.0,
    head = 60,
    efficiency = 0.8,
    turbine = turbine,
)
```

```python
print 'installed capacity is initialised with turbine.power'
print reservoir installed_capacity
```

```python
print '- ' * 15
print 'updating the turbine power updates the installed capacity'
turbine.power = 10
print reservoir installed_capacity
```

```python
print '- ' * 15
print 'setting the installed capacity breaks the link between turbine.power'
print 'and the installed_capacity trait'
```

```python
reservoir installed_capacity = 8
print turbine.power, reservoir installed_capacity
```

### 18.3.6 Notification

Traits implements a Listener pattern. For each trait a list of static and dynamic listeners can be fed with call-backs. When the trait does change, all the listeners are called.

Static listeners are defined using the `_XXX_changed` magic methods:

```python
from traits.api import HasTraits, Instance, DelegatesTo, Float, Range
```

```python
from reservoir import Reservoir
```

```python
class ReservoirState(HasTraits):
    """Keeps track of the reservoir state given the initial storage.
    """
```

```python
reservoir = Instance(Reservoir,())
min_storage = Float
max_storage = DelegatesTo('reservoir')
min_release = Float
max_release = DelegatesTo('reservoir')
```

```python
# state attributes
storage = Range(low='min_storage', high='max_storage')
```

```python
# control attributes
inflows = Float desc='Inflows [hm3]')
release = Range(low='min_release', high='max_release')
spillage = Float desc='Spillage [hm3]')
```
```python
def print_state(self):
    print 'Storage\tRelease\tInflows\tSpillage'
    str_format = '\t'.join(['{:.7.2f}' for i in range(4)])
    print str_format.format(self.storage, self.release, self.inflows,
                             self.spillage)
    print '-' * 79

### Traits listeners ################
def _release_changed(self, new):
    """When the release is higher than zero, warn all the inhabitants of
    the valley.
    """
    if new > 0:
        print 'Warning, we are releasing {} hm3 of water'.format(new)

if __name__ == '__main__':
    projectA = Reservoir(
        name = 'Project A',
        max_storage = 30,
        max_release = 100.0,
        hydraulic_head = 60,
        efficiency = 0.8
    )

    state = ReservoirState(reservoir=projectA, storage=10)
    state.release = 90
    state.inflows = 0
    state.print_state()
```

The static trait notification signatures can be:

* `def_release_changed(self):` pass

* `def_release_changed(self, new):` pass

* `def_release_changed(self, old, new):` pass

* `def_release_changed(self, name, old, new)` pass

## Listening to all the changes

To listen to all the changes on a `HasTraits` class, the magic `_any trait changed` method can be implemented.

In many situations, you do not know in advance what type of listeners need to be activated. `Traits` offers the ability to register listeners on the fly with the dynamic listeners

```python
from reservoir import Reservoir
from reservoir_state_property import ReservoirState

def wake_up watchman if_spillage(new_value):
    if new_value > 0:
        print 'Wake up watchman! Spilling {} hm3'.format(new_value)

if __name__ == '__main__':
    projectA = Reservoir(
        name = 'Project A',
        max_storage = 30,
        max_release = 100.0,
        hydraulic_head = 60,
        efficiency = 0.8
```
```python
)

state = ReservoirState(reservoir=projectA, storage=10)

# register the dynamic listener
state.on trait_change(wake_up watchman if spillage, name='spillage')

state.release = 90
state.inflows = 0
state.print_state()

print 'Forcing spillage'
state.inflows = 100
state.release = 0

print 'Why do we have two executions of the callback ?'
```

The dynamic trait notification signatures are not the same as the static ones :

* `def wake_up watchman()`: pass

* `def wake_up watchman(new)`: pass

* `def wake_up watchman(name, new)`: pass

* `def wake_up watchman(object, name, new)`: pass

* `def wake_up watchman(object, name, old, new)`: pass

Removing a dynamic listener can be done by:

* calling the `remove trait listener` method on the trait with the listener method as argument,

* calling the `on trait change` method with listener method and the keyword `remove=True`,

* deleting the instance that holds the listener.

Listeners can also be added to classes using the `on trait change` decorator:

```python
from traits.api import HasTraits, Instance, DelegatesTo, Float, Range
from traits.api import Property, on trait_change

from reservoir import Reservoir

class ReservoirState(HasTraits):
    """Keeps track of the reservoir state given the initial storage.

    For the simplicity of the example, the release is considered in
    hm3/timestep and not in m3/s.
    """
    reservoir = Instance(Reservoir,())
    max_storage = DelegatesTo('reservoir')
    min_release = Float
    max_release = DelegatesTo('reservoir')

    # state attributes
    storage = Property depends_on='inflows, release')

    # control attributes
    inflows = Float desc='Inflows [hm3]')
    release = Range(low='min_release', high='max_release')
    spillage = Property(
        desc='Spillage [hm3]', depends_on=['storage', 'inflows', 'release']
    )

### Private traits. ########
```
`_storage` = Float

### Traits property implementation.

def _get_storage(self):
    new_storage = self._storage - self.release + self.inflows
    return min(new_storage, self.max_storage)

def _set_storage(self, storage_value):
    self._storage = storage_value

def _get_spillage(self):
    new_storage = self._storage - self.release + self.inflows
    overflow = new_storage - self.max_storage
    return max(overflow, 0)

@on trait_change('storage')

def print_state(self):
    print 'Storage\tRelease\tInflows\tSpillage'
    str_format = '\t'.join(['{:.2f}'] * len(4))
    print str_format.format(self.storage, self.release, self.inflows,
                             self.spillage)
    print '-' * 79

if __name__ == '__main__':
    projectA = Reservoir(

name = 'Project A',
    max_storage = 30,
    max_release = 5,
    hydraulic_head = 60,
    efficiency = 0.8
)

state = ReservoirState(reservoir=projectA, storage=25)
state.release = 4
state.inflows = 0

The patterns supported by the `on trait_change` method and decorator are powerful. The reader should look at the docstring of `HasTraits.on trait_change` for the details.

### 18.3.7 Some more advanced traits

The following example demonstrates the usage of the `Enum` and `List` traits:

from traits.api import HasTraits, Str, Float, Range, Enum, List
from traitsuiapidoc import View, Item

class IrrigationArea(HasTraits):
    name = Str
    surface = FloatSurface [ha]')
    crop = Enum('Alfalfa', 'Wheat', 'Cotton')

class Reservoir(HasTraits):
    name = Str
    max_storage = Float(1e6, desc='Maximal storage [hm3]')
    max_release = Float(10, desc='Maximal release [m3/s]')
    head = Float(10, desc='Hydraulic head [m]')
    efficiency = Range(0, 1.)
    irrigated areas = List(IrrigationArea)

def energy_production(self, release):
    """Returns the energy production [Wh] for the given release [m3/s]"""
```python
''' power = 1000 * 9.81 * self head * release * self.efficiency
return power * 3600
'''

def __main__:
    traits_view = View(
        Item('name'),
        Item('max_storage'),
        Item('max_release'),
        Item('head'),
        Item('efficiency'),
        Item('irrigated areas'),
        resizable = True
    )

if __name__ == '__main__':
    upper_block = IrrigationArea(name='Section C', surface=2000, crop='Wheat')

    reservoir = Reservoir(
        name='Project A',
        max_storage=30,
        max_release=100.0,
        head=60,
        efficiency=0.8,
        irrigated areas=[upper_block]
    )

    release = 80
    print 'Releasing {} m3/s produces {} kWh'.format(
        release, reservoir.energy production(release)
    )
```

Trait listeners can be used to listen to changes in the content of the list to e.g. keep track of the total crop surface on linked to a given reservoir.

```python
from traits.api import HasTraits, Str, Float, Range, Enum, List, Property
from traitsui.api import View, Item

class IrrigationArea(HasTraits):
    name = Str
    surface = Float(description='Surface [ha]')
    crop = Enum('Alfalfa', 'Wheat', 'Cotton')

class Reservoir(HasTraits):
    name = Str
    max_storage = Float(1e6, description='Maximal storage [hm3]')
    max_release = Float(10, description='Maximal release [m3/s]')
    head = Float(10, description='Hydraulic head [m]')
    efficiency = Range(0, 1)

    irrigated areas = List(IrrigationArea)

total_cropsurface = Property depends_on='irrigated areas surface')

def get_total_cropsurface(self):
    return sum([iarea surface for iarea in self.irrigated areas])

def energy production(self, release):
    ''' Returns the energy production [Wh] for the given release [m3/s]
    '''
    power = 1000 * 9.81 * self head * release * self.efficiency
    return power * 3600
```
```python
traits_view = View(
    Item('name'),
    Item('max_storage'),
    Item('max_release'),
    Item('head'),
    Item('efficiency'),
    Item('irrigatedareas'),
    Item('total_cropsurface'),
    resizable = True
)
```

```python
if __name__ == '__main__':
    upper_block = IrrigationArea(name='Section C', surface=2000, crop='Wheat')

    reservoir = Reservoir(
        name='Project A',
        max_storage=30,
        max_release=100.0,
        head=60,
        efficiency=0.8,
        irrigated areas=[upper_block],
    )

    release = 80
    print 'Releasing {} m3/s produces {} kWh'.format(
        release, reservoir.energy_production(release)
    )
```

The next example shows how the Array trait can be used to feed a specialised TraitsUI Item, the ChacoPlotItem:

```python
import numpy as np

from traits.api import HasTraits, Array, Instance, Float, Property
from traits.api import DelegatesTo
from traitsui-api import View, Item, Group
from chaco.chaco_plot_editor import ChacoPlotItem

from reservoir import Reservoir

class ReservoirEvolution(HasTraits):
    reservoir = Instance(Reservoir)

    name = DelegatesTo('reservoir')

    inflows = Array(dtype= np.float64, shape=(None))
    releases = Array(dtype= np.float64, shape=(None))

    initial.stock = Float
    stock = Property depends_on='inflows, releases, initial.stock')

    month = Property depends_on='stock')

### Traits view

traits_view = View(
    Item('name'),
    Group(
        ChacoPlotItem('month', 'stock', show_label=False),
    ),
    width = 500,
    resizable = True
)
```
### Traits properties

def _get.stock(self):
    """
    fixme: should handle cases where we go over the max storage
    """
    return selfinitial�k + (self.inflows - self.releases).cumsum()

def _get_month(self):
    return np. range(self.stock.size)

if __name__ == '__main__':
    reservoir = Reservoir(
        name = 'Project A',
        max_storage = 30,
        max_release = 100.0,
        head = 60,
        efficiency = 0.8
    )

initial�k = 10.
inflows_ts = np. array([6., 6, 4, 4, 1, 2, 0, 0, 3, 1, 5, 3])
releases_ts = np. array([4., 5, 3, 5, 3, 5, 5, 3, 2, 1, 3, 3])

view = ReservoirEvolution(
    reservoir = reservoir,
    inflows = inflows_ts,
    releases = releases_ts
)
view.configure

<figure><img src="image_261.png" /></figure>

See also:

References

* ETS repositories

* Traits manual

* Traits UI manual

* Mailing list: enthought-dev@enthought.com
### 19.1.1 3D plotting functions

#### Points

<figure><img src="image_262.png" /></figure>

Points in 3D, represented with markers (or “glyphs”) and optionally different sizes.

```python
x, y, z, value = np.random.random((4, 40))
mlab points3d(x, y, z, value)
```

#### Lines

<figure><img src="image_263.png" /></figure>

A line connecting points in 3D, with optional thickness and varying color.

```python
mlab.clf()  # Clear the figure
t = np.linspace(0, 20, 200)
mlab.plot3d(np.sin(t), np.cos(t), 0.1*t, t)
```

#### Elevation surface

<figure><img src="image_264.png" /></figure>

A surface given by its elevation, coded as a 2D array
```python
mlab.clf()
x, y = np.mgrid[-10:10:100j, -10:10:100j]
r = np.sqrt(x**2 + y**2)
z = np.sin(r)/r
mlab.surf(z, warp_scale='auto')
```

## Arbitrary regular mesh

<figure><img src="image_265.png" /></figure>

A surface mesh given by x, y, z positions of its node points

```python
mlab.clf()
phi, theta = np.mgrid[0:np.pi:11j, 0:2*np.pi:11j]
x = np.sin(*phi) * np.cos(*theta)
y = np.sin District(*phi) * np.sin District(*theta)
z = np.cos District(*phi)
mlab.mesh(x, y, z)
mlab mesh(x, y, z, representation='wireframe', color=(0, 0, 0))
```

A surface is defined by points **connected** to form triangles or polygons. In `mayavi.mlab.surf()` and `mayavi.mlab mesh()`, the connectivity is implicitly given by the layout of the arrays. See also `mayavi.mlab.triangularMesh()`.

**Our data is often more than points and values: it needs some connectivity information**

## Volumetric data

<figure><img src="image_266.png" /></figure>

If your data is *dense* in 3D, it is more difficult to display. One option is to take iso-contours of the data.

```python
mlab.clf()
x, y, z = np.mgrid[-5:5:64j, -5:5:64j, -5:5:64j]
values = x*x*0.5 + y*y + z*z*2.0
mlab.contour3d(values)
```
<figure><img src="image_267.png" /></figure>

**This function works with a regular orthogonal grid:** the value array is a 3D array that gives the shape of the grid.

## 19.1.2 Figures and decorations

### Figure management

Here is a list of functions useful to control the current figure

<table><tr><td>Get the current figure:</td><td>mlab.gcf()</td></tr><tr><td>Clear the current figure:</td><td>mlab.clf()</td></tr><tr><td>Set the current figure:</td><td>mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0.5, 0.5, 0.5))</td></tr><tr><td>Save figure to image file:</td><td>mlab.savefig('foo.png', size=(300, 300))</td></tr><tr><td>Change the view:</td><td>mlab.view(azimuth=45, elevation=54, distance=1.)</td></tr></table>

### Changing plot properties

In general, many properties of the various objects on the figure can be changed. If these visualization are created via `mlab` functions, the easiest way to change them is to use the keyword arguments of these functions, as described in the docstrings.
**Example docstring: mlab.mesh**

Plots a surface using grid-spaced data supplied as 2D arrays.

**Function signatures:**

`mesh(x, y, z, ...)`

x, y, z are 2D arrays, all of the same shape, giving the positions of the vertices of the surface. The connectivity between these points is implied by the connectivity on the arrays.

For simple structures (such as orthogonal grids) prefer the surf function, as it will create more efficient data structures.

**Keyword arguments:**

**color** the color of the vtk object. Overrides the colormap, if any, when specified.
This is specified as a triplet of float ranging from 0 to 1, eg (1, 1, 1) for white.

**colormap** type of colormap to use.

**extent** [xmin, xmax, ymin, ymax, zmin, zmax] Default is the x, y, z arrays extents.
Use this to change the extent of the object created.

**figure** Figure to populate.

**line_width** The width of the lines, if any used. Must be a float. Default: 2.0

**mask** boolean mask array to suppress some data points.

**mask_points** If supplied, only one out of 'mask_points' data point is displayed.
This option is usefull to reduce the number of points displayed on large
datasets Must be an integer or None.

**mode** the mode of the glyphs. Must be '2darrow' or '2dcircle' or '2dcross' or '2ddash' or '2ddiamond' or '2dhookedArrow' or '2dsquare' or '2dthickArrow' or '2dthick-cross' or '2dtriangle' or '2dvertex' or 'arrow' or 'cone' or 'cube' or 'cylinder' or 'point' or 'sphere'. Default: sphere

**name** the name of the vtk object created.

**representation** the representation type used for the surface. Must be 'surface' or 'wireframe' or 'points' or 'mesh' or 'fancymesh'. Default: surface

**resolution** The resolution of the glyph created. For spheres, for instance, this is the number of divisions along theta and phi. Must be an integer. Default: 8

**scalsars** optional scalar data.

**scale_factor** scale factor of the glyphs used to represent the vertices, in
fancyMesh mode. Must be a float. Default: 0.05

**scale_mode** the scaling mode for the glyphs ('vector', 'scalar', or 'none').

**transparent** make the opacity of the actor depend on the scalar.

**tube radius** radius of the tubes used to represent the lines, in mesh mode. If
None, simple lines are used.

**tube_sides** number of sides of the tubes used to represent the lines. Must be an
integer. Default: 6

**vmax** vmax is used to scale the colormap If None, the max of the data will be used

**vmin** vmin is used to scale the colormap If None, the min of the data will be used
<figure><img src="image_268.png" /></figure>

### Example:

In [1]: import numpy as np

In [2]: r, theta = np.mgrid[0:10, -np.pi[1:10]]

In [3]: x = r * np.cos theta

In [4]: y = r * np.sin theta

In [5]: z = np.sin(r)/r

In [6]: from mayavi import mlab

In [7]: mlab.mesh(x, y, z, colormap='gist_earth', extent=[0, 1, 0, 1, 0, 1])

Out [7]: <mayavi.modules surface Surface object at 0xde6f08c>

In [8]: mlab mesh(x, y, z, extent=[0, 1, 0, 1, 0, 1],

...: representation='wireframe', line_width=1, color=(0.5, 0.5, 0.5))

Out [8]: <mayavi.modules surface Surface object at 0xd6a71c>

### Decorations

Different items can be added to the figure to carry extra information, such as a colorbar or a title.

In [9]: mlab.colorbar(Out[7], orientation='vertical')

Out [9]: <tvtk_classes(scalar bar actor.ScalarBarActor object at 0xd897f8c>

In [10]: mlab.title('polar mesh')

Out [10]: <enthought.mayavi.modules.text.Text object at 0xd8ed38c>

In [11]: mlab.outline(Out[7])

Out [11]: <enthought.mayavi.modules.outlineOutline object at 0xd21b6c>

In [12]: mlab.axes(Out[7])

Out [12]: <enthought.mayavi.modules.axes.Axes object at 0xd2e4bcc>
<figure><img src="image_269.png" /></figure>

<figure><img src="image_270.png" /></figure>

▲ **extent**: If we specified extents for a plotting object, `mlab.outline` and `mlab.axes` don't get them by default.

## 19.2 Interactive work

The quickest way to create beautiful visualization with Mayavi is probably to interactively tweak the various settings.

### 19.2.1 The "pipeline dialog"

Click on the 'Mayavi' button in the scene, and you can control properties of objects with dialogs.

<figure><img src="image_271.png" /></figure>

* Set the background of the figure in the Mayavi Scene node

* Set the colormap in the Colors and legends node

* Right click on the node to add modules or filters
### 19.2.2 The script recording button

To find out what code can be used to program these changes, click on the red button as you modify those properties, and it will generate the corresponding lines of code.

## 19.3 Slicing and dicing data: sources, modules and filters

### 19.3.1 An example: inspecting magnetic fields

Suppose we are simulating the magnetic field generated by Helmholtz coils. The `examples/compute_field.py` script does this computation and gives you a B array, that is (3 x n), where the first axis is the direction of the field (Bx, By, Bz), and the second axis the index number of the point. Arrays X, Y and Z give the positions of these data points.

**Excercise**

Visualize this field. Your goal is to make sure that the simulation code is correct.

<figure><img src="image_272.png" /></figure>

**Suggestions**

* If you compute the norm of the vector field, you can apply an isosurface to it.

* using `mayavi.mlab.quiver3d()` you can plot vectors. You can also use the 'masking' options (in the GUI) to make the plot a bit less dense.

### 19.3.2 Different views on data: sources and modules

As we see above, it may be desirable to look at the same data in different ways.

Mayavi visualization are created by loading the data in a *data source* and then displayed on the screen using *modules*.
This can be seen by looking at the "pipeline" view. By right-clicking on the nodes of the pipeline, you can add new modules.

### Quiz

Why is it not possible to add a `VectorCutPlane` to the vectors created by `mayavi.mlab.quiver3d()`?

## Different sources: scatters and fields

Data comes in different descriptions.

* A 3D block of regularly-spaced value is structured: it is easy to know how one measurement is related to another neighboring and how to continuously interpolate between these. We can call such data a **field**, borrowing from terminology used in physics, as it is continuously defined in space.

* A set of data points measured at random positions in a random order gives rise to much more difficult and ill-posed interpolation problems: the data structure itself does not tell us what are the neighbors of a data point. We call such data a **scatter**.

<figure><img src="image_273.png" /></figure>

Data sources corresponding to **scatters** can be created with `mayavi.mlab.pipeline()` or `mayavi.mlab.pipeline.vector_scatter()`; **field** data sources can be created with `mlab.pipeline()` or `mlab.pipelinevector_field()`.

### Exercise:

1. Create a contour (for instance of the magnetic field norm) by using one of those functions and adding the right *module* by clicking on the GUI dialog.

2. Create the right source to apply a 'vector_cutPlane' and reproduce the picture of the magnetic field shown previously.

Note that one of the difficulties is providing the data in the right form (number of arrays, shape) to the functions. This is often the case with real-life data.

### See also:

Sources are described in details in the *Mayavi manual*.

### Transforming data: filters

If you create a *vector field*, you may want to visualize the iso-contours of its magnitude. But the isosurface module can only be applied to scalar data, and not vector data. We can use a *filter*, `ExtractVectorNorm` to add this scalar value to the vector field.

Filters apply a transformation to data, and can be added between sources and modules
### Excercice

Using the GUI, add the ExtractVectorNorm filter to display iso-contours of the field magnitude.

## mlab.pipeline: the scripting layer

The mlab scripting layer builds pipelines for you. You can reproduce these pipelines programmatically with the mlab.pipeline interface: each step has a corresponding mlab.pipeline function (simply convert the name of the step to lower-case underscore-separated: ExtractVectorNorm gives extract_vector_norm). This function takes as an argument the node that it applies to, as well as optional parameters, and returns the new node.

For example, iso-contours of the magnitude are coded as:

mlab.pipeline.iso_surface(mlab.pipeline.extract_vector_norm(field),
    contours=[0.1*Bmax, 0.4*Bmax],
    opacity=0.5)

<figure><img src="image_274.png" /></figure>

### Excercice

Using the mlab.pipeline interface, generate a complete visualization, with iso-contours of the field magnitude, and a vector cut plane.

(click on the figure for a solution)

## 19.4 Animating the data

To make movies, or interactive application, you may want to change the data represented on a given visualization.

If you have built a visualization, using the mlab plotting functions, or the mlab.pipeline function, we can update the data by assigning new values to the mlab_source attributes

x, y, z = np.Ogrid[-5:5:100j, -5:5:100j, -5:5:100j]
scals = np.sin(x * y * z) / (x * y * z)

iso = mlab.contour3d(scals, transparent=True, contours=[0.5])
for i in range(1, 20):
    scalars = np.sin(i * x * y * z) / (x * y * z)
    iso.mlab_source.scals = scalars
**See also:**

More details in the Mayavi documentation

## Event loops

For the interaction with the user (for instance changing the view with the mouse), Mayavi needs some time to process these events. The for loop above prevents this. The Mayavi documentation details a workaround

# 19.5 Making interactive dialogs

It is very simple to make interactive dialogs with Mayavi using the Traits library (see the dedicated chapter *Traits: building interactive dialogs*).

## 19.5.1 A simple dialog

```python
from traits.api import HasTraits, Instance
from traitsui-api import View, Item, HGroup
from mayavi.core.ui-api import SceneEditor, MlabSceneModel

def curve(n Turns):
    "The function creating the x, y, z coordinates needed to plot"
    phi = np.linspace(0, 2*np.pi, 2000)
    return [np.cosenneurate(0, 2*np.pi, 2000)
            np.sinenneurate(0, 2*np.pi, 2000)
            0.5*np.sin(n turns*phi)]

class Visualization(HasTraits):
    "The class that contains the dialog"
    scene = Instance(MlabSceneModel,())

    def __init__(self):
        HasTraits.__init__(self)
        x, y, z = curve(n turns=2)
        # Populating our plot
        self.plot = self scene.mlab plot3d(x, y, z)

    # Describe the dialog
    view = Viewasser('scene', height=300, show_label=False,
                      editor=SceneEditor,
                      HGroup('n turns'), shareable=True)

    # Fire up the dialog
    Visualization().configure traits()
```

Let us read a bit the code above (examples/mlab_dialog.py).

First, the curve function is used to compute the coordinate of the curve we want to plot.

Second, the dialog is defined by an object inheriting from HasTraits, as it is done with *Traits*. The important point here is that a Mayavi scene is added as a specific *Traits* attribute (*Instance*). This is important for embedding it in the dialog. The view of this dialog is defined by the *view* attribute of the object. In the init of this object, we populate the 3D scene with a curve.

Finally, the `configure廳` method creates the dialog and starts the event loop.

**See also:**

There are a few things to be aware of when doing dialogs with Mayavi. Please read the Mayavi documentation
### 19.5.2 Making it interactive

We can combine the *Traits events handler* with the `mlab_source` to modify the visualization with the dialog.

We will enable the user to vary the `n Turns` parameter in the definition of the curve. For this, we need:

* to define an `n turns` attribute on our visualization object, so that it can appear in the dialog. We use a `Range` type.

* to wire modification of this attribute to a recomputation of the curve. For this, we use the `on traits change` decorator.

<figure><img src="image_275.png" /></figure>

```python
from traits.api import Range, on trait change

class Visualization(HasTraits):
    n turns = Range(0, 30, 11)
    scene = Instance(MlabSceneModel,())

    def __init__(self):
        HasTraits.__init__(self)
        x, y, z = curve(self.n turns)
        self plot = self scene.mlab plot 3d(x, y, z)

    @on trait change('n turns')
    def update plot(self):
        x, y, z = curve(self.n turns)
        self plot.mlab source.set(x=x, y=y, z=z)

    view = View( Item('scene', height=300, show_label=False,
                       editor=SceneEditor()),
                   HGroup('n turns'), resizable=True)

# Fire up the dialog
Visualization().configure traits()
```

Full code of the example: `examples/mlab_dialog.py`.

## 19.6 Putting it together
### Exercise

Using the code from the magnetic field simulation, create a dialog that enables to move the 2 coils: change their parameters.

Hint: to define a dialog entry for a vector of dimension 3

```
direction = Array(float, value=(0, 0, 1), cols=3, shape=(3,))
```

You can look at the *example_coil_application.py* to see a full-blown application for coil design in 270 lines of code.
scikit-learn: machine learning in Python

Authors: Fabian Pedregosa, Gael Varoquaux

<figure><img src="image_276.png" /></figure>

Prerequisites

- Numpy, Scipy

- IPython

- matplotlib

- scikit-learn (http://scikit-learn.org)

Chapters contents

• Loading an example dataset

- Learning and Predicting

• Classification

- k-Nearest neighbors classifier

- Support vector machines (SVMs) for classification

- Clustering: grouping observations together

- K-means clustering

- Dimension Reduction with Principal Component Analysis

- Putting it all together: face recognition

- Linear model: from regression to sparsity

- Sparse models

- Model selection: choosing estimators and their parameters

- Grid-search and cross-validated estimators
## 20.1 Loading an example dataset

<figure><img src="image_277.png" /></figure>

First we will load some data to play with. The data we will use is a very simple flower database known as the Iris dataset.

We have 150 observations of the iris flower specifying some measurements: sepal length, sepal width, petal length and petal width together with its subtype: *Iris setosa*, *Iris versicolor*, *Iris virginica*.

To load the dataset into a Python object:

```python
>>> from sklearn import datasets
```

```python
>>> iris = datasets.load_iris()
```

This data is stored in the `.data` member, which is a (n Samples, n Features) array.

```python
>>> iris.data.shape
```

(150, 4)

The class of each observation is stored in the `.target` attribute of the dataset. This is an integer 1D array of length nSamples:

```python
>>> iris.target.shape
```

(150,)

```python
>>> import numpy as np
```

```python
>>> np.unique(iris.target)
```

```
array([0, 1, 2])
```
An example of reshaping data: the digits dataset

<figure><img src="image_278.png" /></figure>

The digits dataset consists of 1797 images, where each one is an 8x8 pixel image representing a hand-written digit

```python
>>> digits = datasets.load_digits()
>>> digits/images.shape
(1797, 8, 8)
>>> import pylab as pl
>>> pl.imshowنتهimages[0], cmap=pl.cm .gray_r)
<matplotlib.image .AxesImage object at ...>
```

To use this dataset with the scikit, we transform each 8x8 image into a vector of length 64

```python
>>> data = digits images .reshape((digits images .shape [0], -1))
```

20.1.1 Learning and Predicting

Now that we've got some data, we would like to learn from it and predict on new one. In scikit-learn, we
learn from existing data by creating an estimator and calling its fit(X, Y) method.

```python
>>> from sklearn import svm
>>> clf = svm.LinearSVC()
>>> clf.fit(iris.data, iris.target) # learn from the data
LinearSVC(...)
```

Once we have learned from the data, we can use our model to predict the most likely outcome on unseen data:

```python
>>> clf Predict([ [5.0, 3.6, 1.3, 0.25]])
array([0])
```

We can access the parameters of the model via its attributes ending with an underscore:

```python
>>> clf.coef_
array([[0.]])
```

20.2 Classification

20.2.1 k-Nearest neighbors classifier

The simplest possible classifier is the nearest neighbor: given a new observation, take the label of the training samples closest to it in n-dimensional space, where n is the number of features in each sample.
<figure><img src="image_279.png" /></figure>

The k-nearest neighbors classifier internally uses an algorithm based on ball trees to represent the samples it is trained on.

**KNN (k-nearest neighbors) classification example:**

```python
>>> # Create and fit a nearest-neighbor classifier
>>> fromsklearnimportneighbors
>>> knn = neighbors.KNeighborsClassifier()
>>> knn.fit(iris.data, iris.target)
KNeighborsClassifier(...)
>>> knn Predict([ [0.1, 0.2, 0.3, 0.4] ])
array([0])
```

**Training set and testing set**

When experimenting with learning algorithms, it is important not to test the prediction of an estimator on the data used to fit the estimator. Indeed, with the kNN estimator, we would always get perfect prediction on the training set.

```python
>>> perm = np.random.permutation(iris.target.size)
>>> iris.data = iris.data[perm]
>>> iris.target = iris.target[perm]
>>> knn.fit(iris.data[:100], iris.target[:100])
KNeighborsClassifier(...)
>>> knn.score(iris.data[100:], iris.target[100:])
0.95999...
```

Bonus question: why did we use a random permutation?

## 20.2.2 Support vector machines (SVMs) for classification

### Linear Support Vector Machines

SVMs try to construct a hyperplane maximizing the margin between the two classes. It selects a subset of the input, called the support vectors, which are the observations closest to the separating hyperplane.
<figure><img src="image_280.png" /></figure>

```python
>>> from sklearn import svm
>>> svc = svm.SVC核='linear')
>>> svc.fit(iris.data, iris.target)
SVC(...)
```

There are several support vector machine implementations in `scikit-learn`. The most commonly used ones are `svm.SVC`, `svm.NuSVC` and `svm.LinearSVC`; "SVC" stands for Support Vector Classifier (there also exist SVMs for regression, which are called "SVR" in `scikit-learn`).

### Excercise

Train an `svm.SVC` on the digits dataset. Leave out the last 10% and test prediction performance on these observations.

## Using kernels

Classes are not always separable by a hyperplane, so it would be desirable to have a decision function that is not linear but that may be for instance polynomial or exponential:

<table><thead><tr><th>Linear kernel</th><th>Polynomial kernel</th><th>RBF kernel (Radial Basis Function)</th></tr></thead><tbody><tr><td></td><td></td><td></td></tr><tr><td>>>> svc = svm.SVC核='linear')svc = svm.SVC核='polynomial', degree=3)>>> # gamma: inverse of size of >>> # degree: polynomial degree>>> # radial kernel</td><td></td><td></td></tr></tbody></table>

### Exercise

Which of the kernels noted above has a better prediction performance on the digits dataset?

## 20.3 Clustering: grouping observations together

Given the iris dataset, if we knew that there were 3 types of iris, but did not have access to their labels, we could try **unsupervised learning**: we could **cluster** the observations into several groups by some criterion.
### 20.3.1 K-means clustering

The simplest clustering algorithm is k-means. This divides a set into *k* clusters, assigning each observation to a cluster so as to minimize the distance of that observation (in *n*-dimensional space) to the cluster's mean; the means are then recomputed. This operation is run iteratively until the clusters converge, for a maximum for `max_iter` rounds.

(An alternative implementation of k-means is available in SciPy's `cluster` package. The `scikit-learn` implementation differs from that by offering an object API and several additional features, including smart initialization.)

```python
>>> from sklearn import cluster, datasets
>>> iris = datasets.load_iris()
>>> k Means = cluster.KMeans(n clusters=3)
>>> k Means.fit(iris.data)
KMeans(...)
>>> print(k Means.labels_[::10])
[1 1 1 1 1 0 0 0 0 0 2 2 2 2]
>>> print(iris.target[:10])
[0 0 0 0 0 1 1 1 1 1 2 2 2 2]
```

<figure><img src="image_281.png" /></figure>

Ground truth

<figure><img src="image_282.png" /></figure>

K-means (3 clusters)

<figure><img src="image_283.png" /></figure>

K-means (8 clusters)

#### Application to Image Compression

Clustering can be seen as a way of choosing a small number of observations from the information. For instance, this can be used to posterize an image (conversion of a continuous gradation of tone to several regions of fewer tones):

```python
>>> from scipy import misc
>>> lena = misc.lena().astype(np.float32)
>>> X = lena.reshape((-1, 1))  # We need an (n_sample, nfeature) array
>>> kmeans = cluster.KMeans(n clusters=5)
>>> kmeans.fit(X)
KMeans(...)
>>> values = kmeans.clusterCenters_. squeeze()
>>> labels = kmeans.labels_
>>> lena_compressed = np choose (labels, values)
>>> lena_compressed.shape = lena.shape
```

<figure><img src="image_284.png" /></figure>

Raw image

<figure><img src="image_285.png" /></figure>

K-means quantization
## 20.4 Dimension Reduction with Principal Component Analysis

<figure><img src="image_286.png" /></figure>

The cloud of points spanned by the observations above is very flat in one direction, so that one feature can almost be exactly computed using the 2 other. PCA finds the directions in which the data is not *flat* and it can reduce the dimensionality of the data by projecting on a subspace.

<figure><img src="image_287.png" /></figure>

Depending on your version of scikit-learn PCA will be in module decomposition or pca.

```python
>>> fromsklearnimportdecomposition
>>> pca = decomposition.PCA(n_components=2)
>>> pca.fit(iris.data)
PCA(copy=True, n_components=2, whiten=False)
>>> X = pca transform(iris.data)
```

Now we can visualize the (transformed) iris dataset:

```python
>>> importpylab as pl
>>> pl(scatter(X[:,0], X[:,1], c=iris.target))
< matplotlib.collections ... Collection object at ... >
```

<figure><img src="image_288.png" /></figure>

PCA is not just useful for visualization of high dimensional datasets. It can also be used as a preprocessing step to help speed up supervised methods that are not efficient with high dimensions.
## 20.5 Putting it all together: face recognition

An example showcasing face recognition using Principal Component Analysis for dimension reduction and Support Vector Machines for classification.

<figure><img src="image_289.png" /></figure>

Stripped-down version of the face recognition example:

```python
import numpy as np
import pylab as pl
from sklearn import cross_val, datasets, decomposition, svm

# ..
# .. load data ..
lfw people = datasets.fetch_lfw people(min_faces_perPerson=70, resize=0.4)
perm = np.random.permutation(lfw people.target.size)
lfw people.data = lfw people.data[perm]
lfw people.target = lfw people.target[perm]
faces = np.reshape(lfw people.data, (lfw people.target.shape[0], -1))
train, test = iter(cross_val.StratifiedKFold(lfw people.target, k=4)).next()
X_train, X_test = faces[train], faces(test]
y_train, y_test = lfw people.target[train], lfw people.target test]

# ..
# .. dimension reduction ..
pca = decomposition.RandomizedPCA(n_components=150, whiten=True)
pca.fit(X_train)
X_train_pca = pca transform(X_train)
X_test_pca = pca transform(X_test)

# ..
# .. classification ..
clf = svm.SVC(C=5., gamma=0.001)
clf.fit(X_train_pca, y_train)
```
```python
# ..
# .. predict on new images ..
for i in range(10):
    print(lfw.people.target_names[clf Predict(X_test_pca[i])[0]])
    _ = pl.imshow(X_test[i].reshape(50, 37), cmap=pl.cm_gray)
    _ = raw_input()
```

## 20.6 Linear model: from regression to sparsity

### Diabetes dataset

The diabetes dataset consists of 10 physiological variables (age, sex, weight, blood pressure) measured on 442 patients, and an indication of disease progression after one year:

```python
>>> diabetes = datasets.load_diabetes()
>>> diabetes_X_train = diabetes.data[:-20]
>>> diabetes_X_test = diabetes.data[-20:]
>>> diabetes_y_train = diabetes.target[:-20]
>>> diabetes_y_test = diabetes.target[-20:]
```

The task at hand is to predict disease prediction from physiological variables.

### 20.6.1 Sparse models

To improve the conditioning of the problem (uninformative variables, mitigate the curse of dimensionality, as a feature selection preprocessing, etc.), it would be interesting to select only the informative features and set non-informative ones to 0. This penalization approach, called **Lasso**, can set some coefficients to zero. Such methods are called **sparse method**, and sparsity can be seen as an application of Occam's razor: prefer simpler models to complex ones.

```python
>>> from sklearn import linear_model
>>> regr = linear_model.Lasso(alpha=.3)
>>> regr.fit(diabetes_X_train, diabetes_y_train)
Lasso(...)
>>> regr.coef_ # very sparse coefficients
array([ 0.      , -0.      ,  497.34075682,  199.17441034,
       -0.      , -0.      , -118.89291545,    0.      ,
      430.9379595,    0.      ])
>>> regr.score(diabetes_X_test, diabetes_y_test)
0.5510835453...
```

being the score very similar to linear regression (Least Squares):

```python
>>> lin = linear_model.LinearRegression()
>>> lin.fit(diabetes_X_train, diabetes_y_train)
LinearRegression(...)
>>> lin.score(diabetes_X_test, diabetes_y_test)
0.5850753022...
```

### Different algorithms for a same problem

Different algorithms can be used to solve the same mathematical problem. For instance the Lasso object in the `sklearn` solves the lasso regression using a coordinate descent method, that is efficient on large datasets. However, the `sklearn` also provides the `LassoLARS` object, using the `LARS` which is very efficient for problems in which the weight vector estimated is very sparse, that is problems with very few observations.
## 20.7 Model selection: choosing estimators and their parameters

### 20.7.1 Grid-search and cross-validated estimators

#### Grid-search

The scikit-learn provides an object that, given data, computes the score during the fit of an estimator on a parameter grid and chooses the parameters to maximize the cross-validation score. This object takes an estimator during the construction and exposes an estimator API:

```python
>>> from sklearn import svm, grid_search
>>> gammas = np.logspace(-6, -1, 10)
>>> svc = svm.SVC()
>>> clf = grid_search GridSearchCV estimator = svc, param_grid = dict(gamma = gammas),
...                                                    n_jobs = -1)
>>> clf.fit digits .data[:1000], digits .target[:1000])
GridSearchCV(cv = None, ...)
>>> clf.best_score_
0.9...
>>> clf.best FaceEstimator .gamma
0.00059948425031894088
```

By default the GridSearchCV uses a 3-fold cross-validation. However, if it detects that a classifier is passed, rather than a regressor, it uses a stratified 3-fold.

#### Cross-validated estimators

Cross-validation to set a parameter can be done more efficiently on an algorithm-by-algorithm basis. This is why, for certain estimators, the scikit-learn exposes “CV” estimators, that set their parameter automatically by cross-validation:

```python
>>> from sklearn import linear_model, datasets
>>> lasso = linear_model .LassoCV()
>>> diabetes = datasets.load_diabetes()
>>> X_diabetes = diabetes .data
>>> y_diabetes = diabetes .target
>>> lasso.fit(X_diabetes, y_diabetes)
LassoCV(alphas = None, ...)
>>> # The estimator chose automatically its lambda:
>>> lasso .alpha_
0.012...
```

These estimators are called similarly to their counterparts, with 'CV' appended to their name.

#### Exercise

On the diabetes dataset, find the optimal regularization parameter alpha.
D

diff, 304, 307

differentiation, 304

dsolve, 307

E

equations

algebraic, 305

differential, 307

I

integration, 305

M

Matrix, 306

P

Python Enhancement Proposals

PEP 255, 145

PEP 3118, 182

PEP 3129, 155

PEP 318, 147, 155

PEP 342, 145

PEP 343, 155

PEP 380, 147

PEP 380#id13, 147

PEP 8, 150

S

solve, 305