and you may want to read directly the documentation on the wiki instead of the official documentation website. Note that anyone can create an account on the wiki and write better documentation; this is an easy way to contribute to an open-source project and improve the tools you are using!

* Scipy central http://central.scipy.org/ gives recipes on many common problems frequently encountered, such as fitting data points, solving ODE, etc.

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

* Matplotlib's website http:// matplotlib.org/ features a very nice **gallery** with a large number of plots, each of them shows both the source code and the resulting plot. This is very useful for learning by example. More standard documentation is also available.

Finally, two more “technical” possibilities are useful as well:

* In Ipython, the magical function `%psearch` search for objects matching patterns. This is useful if, for example, one does not know the exact name of a function.

```
In [3]: import numpy as np
In [4]: %psearch np.diag*
np.diag
np.diagflat
np.diagonal
```

* ` numpy lookfor ` looks for keywords inside the docstrings of specified modules.

```
In [45]: numpy lookfor('convolution')
Search results for 'convolution'
--------------------------------------

>>> numpy convolve
Returns the discrete, linear convolution of two one-dimensional
sequences.

>>> numpy bartlett
Return the Bartlett window.

>>> numpy correlate
Discrete, linear correlation of two 1-dimensional sequences.

In [46]: numpy lookfor('remove', module='os')
Search results for 'remove'
--------------------------------------

>>> os.remove
remove(path)
>>> os.removedirs
removedirs(path)
>>> os.rmdir
rmdir(path)
>>> os_unlink
unlink(path)
>>> osWalk
Directory tree generator.
```

* If everything listed above fails (and Google doesn't have the answer)... don't despair! Write to the mailing-list suited to your problem: you should have a quick answer if you describe your problem well. Experts on scientific python often give very enlightening explanations on the mailing-list.

– Numpy discussion ( numpy-discussion@scipy.org ): all about numpy arrays, manipulating them, indexation questions, etc.
- **SciPy Users List** (scipy-user@scipy.org): scientific computing with Python, high-level data processing, in particular with the scipy package.

- `matplotlib-users@lists.sourceforge.net` for plotting with `matplotlib`.
## Part II

# Advanced topics
This part of the *Scipy lecture notes* is dedicated to advanced usage. It strives to educate the proficient Python coder to be an expert and tackles various specific topics.
# Advanced Python Constructs

**Author** Zbigniew Jędrzejewski-Szmek

This section covers some features of the Python language which can be considered advanced — in the sense that not every language has them, and also in the sense that they are more useful in more complicated programs or libraries, but not in the sense of being particularly specialized, or particularly complicated.

It is important to underline that this chapter is purely about the language itself — about features supported through special syntax complemented by functionality of the Python stdlib, which could not be implemented through clever external modules.

The process of developing the Python programming language, its syntax, is very transparent; proposed changes are evaluated from various angles and discussed via *Python Enhancement Proposals* — *PEPs*. As a result, features described in this chapter were added after it was shown that they indeed solve real problems and that their use is as simple as possible.

## Chapter contents

* Iterators, generator expressions and generators

- Iterators

- Generator expressions

- Generators

- Bidirectional communication

- Chaining generators

* Decorators

- Replacing or tweaking the original object

- Decorators implemented as classes and as functions

- Copying the docstring and other attributes of the original function

- Examples in the standard library

- Deprecation of functions

- A while-loop removing decorator

- A plugin registration system

* Context managers

- Catching exceptions

- Using generators to define context managers

## 7.1 Iterators, generator expressions and generators

### 7.1.1 Iterators
## Simplicity

Duplication of effort is wasteful, and replacing the various home-grown approaches with a standard feature usually ends up making things more readable, and interoperable as well.

Guido van Rossum — Adding Optional Static Typing to Python

An iterator is an object adhering to the `iterator protocol` — basically this means that it has a `next` method, which, when called, returns the next item in the sequence, and when there's nothing to return, raises the `StopIteration` exception.

An iterator object allows to loop just once. It holds the state (position) of a single iteration, or from the other side, each loop over a sequence requires a single iterator object. This means that we can iterate over the same sequence more than once concurrently. Separating the iteration logic from the sequence allows us to have more than one way of iteration.

Calling the `__iter__` method on a container to create an iterator object is the most straightforward way to get hold of an iterator. The `iter` function does that for us, saving a few keystrokes.

```python
>>> nums = [1, 2, 3]      # note that ... varies: these are different objects
>>> iter(nums)
<...iterator object at ...>
>>> nums.__iter__()
<...iterator object at ...>
>>> nums.__reversed__()
<...reverseiterator object at ...>
>>> it = iter(nums)
>>> next(it)
1
>>> next(it)
2
>>> next(it)
3
>>> next(it)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
StopIteration
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
StopIteration
```

When used in a loop, `StopIteration` is swallowed and causes the loop to finish. But with explicit invocation, we can see that once the iterator is exhausted, accessing it raises an exception.

Using the `for..in` loop also uses the `__iter__` method. This allows us to transparently start the iteration over a sequence. But if we already have the iterator, we want to be able to use it in an `for` loop in the same way. In order to achieve this, iterators in addition to `next` are also required to have a method called `__iter__` which returns the iterator (`self`).

Support for iteration is pervasive in Python: all sequences and unordered containers in the standard library allow this. The concept is also stretched to other things: e.g. `file` objects support iteration over lines.

```python
>>> f = open('/etc/fstab')
>>> f is f.__iter__()
True
```

The `file` is an iterator itself and its `__iter__` method doesn't create a separate object: only a single thread of sequential access is allowed.
### 7.1.2 Generator expressions

A second way in which iterator objects are created is through **generator expressions**, the basis for **list comprehensions**. To increase clarity, a generator expression must always be enclosed in parentheses or an expression. If round parentheses are used, then a generator iterator is created. If rectangular parentheses are used, the process is short-circuited and we get a **list**.

```python
>>> (i for i in nums)
<generator object <genexpr> at 0x...>
>>> [i for i in nums]
[1, 2, 3]
>>> list(i for i in nums)
[1, 2, 3]
```

The list comprehension syntax also extends to **dictionary and set comprehensions**. A set is created when the generator expression is enclosed in curly braces. A dict is created when the generator expression contains "pairs" of the form `key:value`:

```python
>>> {i for i in range(3)}
set([0, 1, 2])
>>> {i:i**2 for i in range(3)}
{0: 0, 1: 1, 2: 4}
```

One *gotcha* should be mentioned: in old Python's the index variable (i) would leak, and in versions >= 3 this is fixed.

### 7.1.3 Generators

#### Generators

A generator is a function that produces a sequence of results instead of a single value.

David Beazley — A Curious Course on Coroutines and Concurrency

A third way to create iterator objects is to call a **generator function**. A **generator** is a function containing the keyword `yield`. It must be noted that the mere presence of this keyword completely changes the nature of the function: this `yield` statement doesn't have to be invoked, or even reachable, but causes the function to be marked as a generator. When a normal function is called, the instructions contained in the body start to be executed. When a generator is called, the execution stops before the first instruction in the body. An invocation of a generator function creates a generator object, adhering to the iterator protocol. As with normal function invocations, concurrent and recursive invocations are allowed.

When `next` is called, the function is executed until the first `yield`. Each encountered `yield` statement gives a value becomes the return value of `next`. After executing the `yield` statement, the execution of this function is suspended.

```python
>>> def f():
...     yield 1
...     yield 2
...
>>> f()
<generator object f at 0x...>
>>> gen = f()
>>> next(gen)
1
>>> next(gen)
2
>>> next(gen)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
StopIteration
```
Let's go over the life of the single invocation of the generator function.

```python
>>> def f():
...     print("-- start --")
...     yield 3
...     print("-- middle --")
...     yield 4
...     print("-- finished --")
>>> gen = f()
>>> next(gen)
-- start --
3
>>> next(gen)
-- middle --
4
>>> next(gen)
-- finished --
Traceback (most recent call last):
...
StopIteration
Traceback (most recent call last):
...
StopIteration
```

Contrary to a normal function, where executing `f()` would immediately cause the first `print` to be executed, `gen` is assigned without executing any statements in the function body. Only when `gen.next()` is invoked by `next`, the statements up to the first `yield` are executed. The second `next` prints `-- middle --` and execution halts on the second `yield`. The third `next` prints `-- finished --` and falls of the end of the function. Since no `yield` was reached, an exception is raised.

What happens with the function after a `yield`, when the control passes to the caller? The state of each generator is stored in the generator object. From the point of view of the generator function, it looks almost as if it was running in a separate thread, but this is just an illusion: execution is strictly single-threaded, but the interpreter keeps and restores the state in between the requests for the next value.

Why are generators useful? As noted in the parts about iterators, a generator function is just a different way to create an iterator object. Everything that can be done with `yield` statements, could also be done with `next` methods. Nevertheless, using a function and having the interpreter perform its magic to create an iterator has advantages. A function can be much shorter than the definition of a class with the required `next` and `__iter__` methods. What is more important, it is easier for the author of the generator to understand the state which is kept in local variables, as opposed to instance attributes, which have to be used to pass data between consecutive invocations of `next` on an iterator object.

A broader question is why are iterators useful? When an iterator is used to power a loop, the loop becomes very simple. The code to initialise the state, to decide if the loop is finished, and to find the next value is extracted into a separate place. This highlights the body of the loop — the interesting part. In addition, it is possible to reuse the iterator code in other places.

### 7.1.4 Bidirectional communication

Each `yield` statement causes a value to be passed to the caller. This is the reason for the introduction of generators by PEP 255 (implemented in Python 2.2). But communication in the reverse direction is also useful. One obvious way would be some external state, either a global variable or a shared mutable object. Direct communication is possible thanks to PEP 342 (implemented in 2.5). It is achieved by turning the previously boring `yield` statement into an expression. When the generator resumes execution after a `yield` statement, the caller can call a method on the generator object to either pass a value **into** the generator, which then is returned by the `yield` statement, or a different method to inject an exception into the generator.

The first of the new methods is `send(value)`, which is similar to `next()`, but passes value into the generator to be used for the value of the `yield` expression. In fact, `g.next()` and `g.send(None)` are equivalent.

The second of the new methods is `throw(type, value=None, traceback=None)` which is equivalent to:
```python
raise type, value, traceback
```

at the point of the `yield` statement.

Unlike `raise` (which immediately raises an exception from the current execution point), `throw()` first resumes the generator, and only then raises the exception. The word throw was picked because it is suggestive of putting the exception in another location, and is associated with exceptions in other languages.

What happens when an exception is raised inside the generator? It can be either raised explicitly or when executing some statements or it can be injected at the point of a `yield` statement by means of the `throw()` method. In either case, such an exception propagates in the standard manner: it can be intercepted by an `except` or `finally` clause, or otherwise it causes the execution of the generator function to be aborted and propagates in the caller.

For completeness' sake, it's worth mentioning that generator iterators also have a `close()` method, which can be used to force a generator that would otherwise be able to provide more values to finish immediately. It allows the generator `__del__` method to destroy objects holding the state of generator.

Let's define a generator which just prints what is passed in through `send` and `throw`.

```python
>>> import itertools
>>> def g():
...     print '--start--'
...     for i in itertools.count():
...         print '--yielding %i--' % i
...         try:
...             ans = yield i
...         except GeneratorExit:
...             print '--closing--'
...             raise
...         except Exception as e:
...             print '--yield raised %r--' % e)
...         else:
...             print '--yield returned %s--' % ans)
>>> it = g()
>>> next (it)
--start--
--yielding 0--
0
>>> it.send (11)
--yield returned 11--
--yielding 1--
1
>>> it.throw (IndexError)
--yield raisedIndexError()--
--yielding 2--
2
>>> it.close()
--closing--
```

### next or __next__?

In Python 2.x, the iterator method to retrieve the next value is called `next`. It is invoked implicitly through the global function `next`, which means that it should be called `__next__`. Just like the global function `iter` calls `__iter__`. This inconsistency is corrected in Python 3.x, where `it.next` becomes `it.__next__`. For other generator methods — `send` and `throw` — the situation is more complicated, because they are not called implicitly by the interpreter. Nevertheless, there's a proposed syntax extension to allow `continue` to take an argument which will be passed to `send` of the loop's iterator. If this extension is accepted, it's likely that `gen.send` will become `gen.__send__`. The last of generator methods, `close`, is pretty obviously named incorrectly, because it is already invoked implicitly.
### 7.1.5 Chaining generators

This is a preview of PEP 380 (not yet implemented, but accepted for Python 3.3).

Let's say we are writing a generator and we want to yield a number of values generated by a second generator, a **subgenerator**. If yielding of values is the only concern, this can be performed without much difficulty using a loop such as

```python
subgen = some_other_generator()
for v in subgen:
    yield v
```

However, if the subgenerator is to interact properly with the caller in the case of calls to `send()`, `throw()` and `close()`, things become considerably more difficult. The `yield` statement has to be guarded by a `try..except..finally` structure similar to the one defined in the previous section to "debug" the generator function. Such code is provided in PEP 380, here it suffices to say that new syntax to properly yield from a subgenerator is being introduced in Python 3.3:

```python
yield from some_other_generator()
```

This behaves like the explicit loop above, repeatedly yielding values from `some_other_generator` until it is exhausted, but also forwards `send`, `throw` and `close` to the subgenerator.

## 7.2 Decorators

### Summary

This amazing feature appeared in the language almost apologetically and with concern that it might not be that useful.

Bruce Eckel — An Introduction to Python Decorators

Since a function or a class are objects, they can be passed around. Since they are mutable objects, they can be modified. The act of altering a function or class object after it has been constructed but before is is bound to its name is called decorating.

There are two things hiding behind the name “decorator” — one is the function which does the work of decorating, i.e. performs the real work, and the other one is the expression adhering to the decorator syntax, i.e. an at-symbol and the name of the decorating function.

Function can be decorated by using the decorator syntax for functions:

```python
@decorator          # ②
def function():      # ①
    pass
```

* A function is defined in the standard way. ①

* An expression starting with `@` placed before the function definition is the decorator `②`. The part after `@` must be a simple expression, usually this is just the name of a function or class. This part is evaluated first, and after the function defined below is ready, the decorator is called with the newly defined function object as the single argument. The value returned by the decorator is attached to the original name of the function.

Decorators can be applied to functions and to classes. For classes the semantics are identical — the original class definition is used as an argument to call the decorator and whatever is returned is assigned under the original name.

Before the decorator syntax was implemented (PEP 318), it was possible to achieve the same effect by assigning the function or class object to a temporary variable and then invoking the decorator explicitly and then assigning the return value to the name of the function. This sounds like more typing, and it is, and also the
name of the decorated function doubling as a temporary variable must be used at least three times, which is
prone to errors. Nevertheless, the example above is equivalent to:

```python
def function():
    # ①
    pass

function = decorator(function) # ②
```

Decorators can be stacked — the order of application is bottom-to-top, or inside-out. The semantics are such that the originally defined function is used as an argument for the first decorator, whatever is returned by the first decorator is used as an argument for the second decorator, ..., and whatever is returned by the last decorator is attached under the name of the original function.

The decorator syntax was chosen for its readability. Since the decorator is specified before the header of the function, it is obvious that its is not a part of the function body and its clear that it can only operate on the whole function. Because the expression is prefixed with @ is stands out and is hard to miss ("in your face", according to the PEP :)). When more than one decorator is applied, each one is placed on a separate line in an easy to read way.

### 7.2.1 Replacing or tweaking the original object

Decorators can either return the same function or class object or they can return a completely different object. In the first case, the decorator can exploit the fact that function and class objects are mutable and add attributes, e.g. add a docstring to a class. A decorator might do something useful even without modifying the object, for example register the decorated class in a global registry. In the second case, virtually anything is possible: when something different is substituted for the original function or class, the new object can be completely different. Nevertheless, such behaviour is not the purpose of decorators: they are intended to tweak the decorated object, not do something unpredictable. Therefore, when a function is “decorated” by replacing it with a different function, the new function usually calls the original function, after doing some preparatory work. Likewise, when a class is “decorated” by replacing if with a new class, the new class is usually derived from the original class. When the purpose of the decorator is to do something “every time”, like to log every call to a decorated function, only the second type of decorators can be used. On the other hand, if the first type is sufficient, it is better to use it, because it is simpler.

### 7.2.2 Decorators implemented as classes and as functions

The only requirement on decorators is that they can be called with a single argument. This means that decorators can be implemented as normal functions, or as classes with a `__call__` method, or in theory, even as lambda functions.

Let's compare the function and class approaches. The decorator expression (the part after @) can be either just a name, or a call. The bare-name approach is nice (less to type, looks cleaner, etc.), but is only possible when no arguments are needed to customise the decorator. Decorators written as functions can be used in those two cases:

```python
>>> def simpledecorator(function):
...     print("doing decoration")
...     return function
>>> @simpledecorator
... def function():
...     print("inside function")
doing decoration
>>> function()
inside function

>>> def decorator_with_argument(arg):
...     print("defining the decorator")
...     def _decorator(function):
...         # in this inner function, arg is available too
...         print("doing decoration, %r" % arg)
...     return function
```
... return _decorator
>>> @decorator_with_argument
... def function():
...     print("inside function")
defining the decorator
doing decoration, 'abc'
>>> function()
inside function

The two trivial decorators above fall into the category of decorators which return the original function. If they
were to return a new function, an extra level of nestedness would be required. In the worst case, three levels of
nested functions.

>>> def replacingdecorator_with_args(arg):
...     print("defining the decorator")
...     def _decorator(function):
...         # in this inner function, arg is available too
...         print("doing decoration, %r" % arg)
...         def _wrapper(*args, **kwargs):
...             print("inside wrapper, %r %r" % (args, kwargs))
...             return function(*args, **kwargs)
...         return _wrapper
...     return _decorator
>>> @replacingdecorator_with_args("abc")
... def function(*args, **kwargs):
...     print("inside function, %r %r" % (args, kwargs))
...     return 14
defining the decorator
doing decoration, 'abc'
>>> function(11, 12)
inside wrapper, (11, 12) {}
inside function, (11, 12) {}
14

The `_wrapper` function is defined to accept all positional and keyword arguments. In general we cannot know
what arguments the decorated function is supposed to accept, so the wrapper function just passes everything
to the wrapped function. One unfortunate consequence is that the apparent argument list is misleading.

Compared to decorators defined as functions, complex decorators defined as classes are simpler. When an
object is created, the `__init__` method is only allowed to return `None`, and the type of the created object
cannot be changed. This means that when a decorator is defined as a class, it doesn't make much sense to use
the argument-less form: the final decorated object would just be an instance of the decorating class, returned
by the constructor call, which is not very useful. Therefore it's enough to discuss class-based decorators where
arguments are given in the decorator expression and the decorator `__init__` method is used for decorator
construction.

>>> class decorator_class(object):
...     def __init__(self, arg):
...         # this method is called in the decorator expression
...         print("in decorator init, %s" % arg)
...         self.arg = arg
...     def __call__(self, function):
...         # this method is called to do the job
...         print("in decorator call, %s" % self.arg)
...         return function
>>> deco_instance = decorator_class('foo')
in decorator init, foo
>>> @deco_instance
...     def function(*args, **kwargs):
...         print("in function, %s %s" % (args, kwargs))
in decorator call, foo
>>> function()
in function, () {}
Contrary to normal rules (PEP 8) decorators written as classes behave more like functions and therefore their name often starts with a lowercase letter.

In reality, it doesn't make much sense to create a new class just to have a decorator which returns the original function. Objects are supposed to hold state, and such decorators are more useful when the decorator returns a new object.

```python
>>> class replacingdecorator_class(object):
...     def __init__(self, arg):
...         # this method is called in the decorator expression
...         print("in decorator init, %s" % arg)
...         self.arg = arg
...     def __call__(self, function):
...         # this method is called to do the job
...         print("in decorator call, %s" % self.arg)
...         self.function = function
...         return self._wrapper
...     def _wrapper(self, *args, **kwargs):
...         print("in the wrapper, %s %s" % (args, kwargs))
...         return self.function(*args, **kwargs)
>>> deco_instance = replacingdecorator_class('foo')
in decorator init, foo
>>> @deco_instance
...     def function(*args, **kwargs):
...         print("in function, %s %s" % (args, kwargs))
in decorator call, foo
>>> function(11, 12)
in the wrapper, (11, 12) {}
in function, (11, 12) {}
```

A decorator like this can do pretty much anything, since it can modify the original function object and mangle the arguments, call the original function or not, and afterwards mangle the return value.

### 7.2.3 Copying the docstring and other attributes of the original function

When a new function is returned by the decorator to replace the original function, an unfortunate consequence is that the original function name, the original docstring, the original argument list are lost. Those attributes of the original function can partially be “transplanted” to the new function by setting `__doc__` (the docstring), `__module__` and `__name__` (the full name of the function), and `__annotations__` (extra information about arguments and the return value of the function available in Python 3). This can be done automatically by using `functools.update.wrapper`.
functools.update.wrapper.wrapper, wrapped)

"Update a wrapper function to look like the wrapped function."

```python
>>> import functools
>>> def replacingdecorator_with_args(arg):
...     print("defining the decorator")
...     def _decorator(function):
...         print("doing decoration, %r" % arg)
...         def _wrapper(*args, **kwargs):
...             print("inside wrapper, %r %r" % (args, kwargs))
...             return function(*args, **kwargs)
...         return functools.update.wrapper(_wrapper, function)
...     return _decorator
>>> @replacingdecorator_with_args("abc")
... def function():
...     "extensive documentation"
...     print("inside function")
...     return 14
defining the decorator
doing decoration, 'abc'
>>> function
<function function at 0x...>
>>> print(function.__doc__)
extensive documentation
```

One important thing is missing from the list of attributes which can be copied to the replacement function: the argument list. The default values for arguments can be modified through the `__defaults__`, `__kwdefaults__` attributes, but unfortunately the argument list itself cannot be set as an attribute. This means that `help(function)` will display a useless argument list which will be confusing for the user of the function. An effective but ugly way around this problem is to create the wrapper dynamically, using `eval`. This can be automated by using the external decorator module. It provides support for the decorator decorator, which takes a wrapper and turns it into a decorator which preserves the function signature.

To sum things up, decorators should always use `functools.update.wrapper` or some other means of copying function attributes.

### 7.2.4 Examples in the standard library

First, it should be mentioned that there's a number of useful decorators available in the standard library. There are three decorators which really form a part of the language:

* `classmethod` causes a method to become a "class method", which means that it can be invoked without creating an instance of the class. When a normal method is invoked, the interpreter inserts the instance object as the first positional parameter, `self`. When a class method is invoked, the class itself is given as the first parameter, often called `cls`.

Class methods are still accessible through the class' namespace, so they don't pollute the module's namespace. Class methods can be used to provide alternative constructors:

```python
class Array(object):
    def __init__(self, data):
        self.data = data

    @classmethod
    def fromfile(cls, file):
        data = numpy.load(file)
        return cls(data)
```

This is cleaner then using a multitude of flags to `__init__`.
* `staticmethod` is applied to methods to make them “static”, i.e. basically a normal function, but accessible through the class namespace. This can be useful when the function is only needed inside this class (its name would then be prefixed with `_`), or when we want the user to think of the method as connected to the class, despite an implementation which doesn’t require this.

* `property` is the pythonic answer to the problem of getters and setters. A method decorated with `property` becomes a getter which is automatically called on attribute access.

```python
>>> class A(object):
...     @property
...     def a(self):
...         "an important attribute"
...         return "a value"
...
>>> A.a
<property object at 0x...>
>>> A().a
'a value'
```

In this example, `A.a` is an read-only attribute. It is also documented: `help(A)` includes the docstring for attribute `a` taken from the getter method. Defining `a` as a `property` allows it to be calculated on the fly, and has the side effect of making it read-only, because no setter is defined.

To have a setter and a getter, two methods are required, obviously. Since Python 2.6 the following syntax is preferred:

```python
class Rectangle(object):
    def __init__(self, edge):
        self.edge = edge

    @property
    def area(self):
        """Computed area.
        """
        return self.edge**2

    @area.setter
    def area(self, area):
        self.edge = area ** 0.5
```

The way that this works, is that the `property` decorator replaces the getter method with a property object. This object in turn has three methods, `getter`, `setter`, and `deleter`, which can be used as decorators. Their job is to set the getter, setter and deleter of the property object (stored as attributes `fget`, `fset`, and `fdel`). The getter can be set like in the example above, when creating the object. When defining the setter, we already have the property object under `area`, and we add the setter to it by using the `setter` method. All this happens when we are creating the class.

Afterwards, when an instance of the class has been created, the property object is special. When the interpreter executes attribute access, assignment, or deletion, the job is delegated to the methods of the property object.

To make everything crystal clear, let's define a “debug” example:

```python
>>> class D(object):
...     @property
...     def a(self):
...         print("getting 1")
...         return 1
...
...     @a.setter
...     def a(self, value):
...         print("setting %r" % value)
...         @a.deleter
...         def a(self):
```
```python
...
    print("deleting")
>>> D.a
<property object at 0x...>
>>> D.a.fget
<function ...>
>>> D.a.fset
<function ...>
>>> D.a.fdel
<function ...>
>>> d = D()
# ... varies, this is not the same `a` function
>>> d.a
getting 1
1
>>> d.a = 2
setting 2
>>> del d.a
deleting
>>> d.a
getting 1
1
```

Properties are a bit of a stretch for the decorator syntax. One of the premises of the decorator syntax — that the name is not duplicated — is violated, but nothing better has been invented so far. It is just good style to use the same name for the getter, setter, and deleter methods.

Some newer examples include:

* `functools.lru_cache` memoizes an arbitrary function maintaining a limited cache of arguments: answer pairs (Python 3.2)

* `functools.total_ordering` is a class decorator which fills in missing ordering methods (`__lt__`, `__gt__`, `__le__`, ...) based on a single available one (Python 2.7).

### 7.2.5 Deprecation of functions

Let's say we want to print a deprecation warning on stderr on the first invocation of a function we don't like anymore. If we don't want to modify the function, we can use a decorator:

```python
class deprecated(object):
    """Print a deprecation warning once on first use of the function.

    >>> @deprecated()
    ... def f():
    ...     pass
    >>> f()
    f is deprecated
    """
    def _call_(self, func):
        self.func = func
        self.count = 0
        return self._wrapper
    def _wrapper(self, *args, **kwargs):
        self.count += 1
        if self.count == 1:
            print(self.func.__name__, 'is deprecated')
        return self.func(*args, **kwargs)
```

It can also be implemented as a function:

```python
def deprecated(func):
    """Print a deprecation warning once on first use of the function.

    >>> @deprecated
    # doctest: +SKIP
```
```python
... def f():
...     pass
>>> f()
# doctest: +SKIP
f is deprecated
"
"
count = [0]
def wrapper(*args, **kwargs):
    count[0] += 1
    if count[0] == 1:
        print func.__name__, 'is deprecated'
    return func(*args, **kwargs)
return wrapper
```

### 7.2.6 A while-loop removing decorator

Let's say we have function which returns a list of things, and this list created by running a loop. If we don't know how many objects will be needed, the standard way to do this is something like:

```python
def find():
    answers = []
    while True:
        ans = look_for_next_answer()
        if ans is None:
            break
        answers.append ans)
    return answers
```

This is fine, as long as the body of the loop is fairly compact. Once it becomes more complicated, as often happens in real code, this becomes pretty unreadable. We could simplify this by using `yield` statements, but then the user would have to explicitly call `list.find answers()`.

We can define a decorator which constructs the list for us:

```python
def vectorized(generator_func):
    def wrapper(*args, **kwargs):
        return list(generator_func(*args, **kwargs))
    return functools.update.wrapper.wrapper, generator_func)
```

Our function then becomes:

```python
@vectorized
def find answers():
    while True:
        ans = look_for_next_answer()
        if ans is None:
            break
    yield ans
```

### 7.2.7 A plugin registration system

This is a class decorator which doesn't modify the class, but just puts it in a global registry. It falls into the category of decorators returning the original object:

```python
class WordProcessor(object):
    PLUGINS = []
    def process(self, text):
        for plugin in self.PLUGINS:
            text = plugin().cleanup(text)
        return text
```
```python
@classmethod
def plugin(cls, plugin):
    cls.PLUGINS.append plugin)

@WordProcessor.plugin
class CleanMdashesExtension(object):
    def cleanup(self, text):
        return text.replace('&mdash; ', u'\N{em dash}')
```

Here we use a decorator to decentralise the registration of plugins. We call our decorator with a noun, instead of a verb, because we use it to declare that our class is a plugin for WordProcessor. Method plugin simply appends the class to the list of plugins.

A word about the plugin itself: it replaces HTML entity for em-dash with a real Unicode em-dash character. It exploits the `unicode literal notation` to insert a character by using its name in the unicode database ("EM DASH"). If the Unicode character was inserted directly, it would be impossible to distinguish it from an en-dash in the source of a program.

**See also:**

**More examples and reading**

* PEP 318 (function and method decorator syntax)

* PEP 3129 (class decorator syntax)

* http://wiki.python.org/moin/PythonDecoratorLibrary

* https://docs.python.org/dev/library/functools.html

* http://pypi.python.org/pypi/decorator

* Bruce Eckel

- Decorators I: Introduction to Python Decorators

- Python Decorators II: Decorator Arguments

- Python Decorators III: A Decorator-Based Build System

## 7.3 Context managers

A context manager is an object with `__enter__` and `__exit__` methods which can be used in the `with` statement:

```python
with manager as var:
    do_something(var)
```

is in the simplest case equivalent to

```python
var = manager.__enter__()
try:
    do_something(var)
finally:
    manager.__exit__()
```

In other words, the context manager protocol defined in PEP 343 permits the extraction of the boring part of a `try..except..finally` structure into a separate class leaving only the interesting `do_something` block.

1. The `__enter__` method is called first. It can return a value which will be assigned to `var`. The as-part is optional: if it isn't present, the value returned by `__enter__` is simply ignored.

2. The block of code underneath `with` is executed. Just like with `try` clauses, it can either execute successfully to the end, or it can `break`, `continue` or `return`, or it can throw an exception. Either way, after the block is finished, the `__exit__` method is called. If an exception was thrown, the information about the
exception is passed to `__exit__()`, which is described below in the next subsection. In the normal case, exceptions can be ignored, just like in a finally clause, and will be rethrown after `__exit__` is finished.

Let's say we want to make sure that a file is closed immediately after we are done writing to it:

```python
>>> class closing(object):
...     def __init__(self, obj):
...         self.obj = obj
...     def __enter__(self):
...         return self.obj
...     def __exit__(self, *args):
...         self.obj.close()
>>> with closing(open('/tmp/file', 'w')) as f:
...     f.write('the contents\n')
```

Here we have made sure that the `f.close()` is called when the `with` block is exited. Since closing files is such a common operation, the support for this is already present in the `file` class. It has an `__exit__` method which calls `close` and can be used as a context manager itself:

```python
>>> with open('/tmp/file', 'a') as f:
...     f.write('more contents\n')
```

The common use for `try...finally` is releasing resources. Various different cases are implemented similarly: in the `__enter__` phase the resource is acquired, in the `__exit__` phase it is released, and the exception, if thrown, is propagated. As with files, there's often a natural operation to perform after the object has been used and it is most convenient to have the support built in. With each release, Python provides support in more places:

* all file-like objects:

- `file` ➔ automatically closed

- `fileinput`, `tempfile` (py >= 3.2)

- `bz2.BZ2File`, `gzip.GzipFile`, `tarfile.TarFile`, `zipfile.ZipFile`

- `ftplib`, `nntplib` ➔ close connection (py >= 3.2 or 3.3)

* locks

- `multiprocessing.RLock` ➔ lock and unlock

- `multiprocessing.Semaphore`

- `memoryview` ➔ automatically release (py >= 3.2 and 2.7)

* `decimal.localcontext` ➔ modify precision of computations temporarily

* `_winreg.PyHKEY` ➔ open and close hive key

* `warnings catcher_warnings` ➔ kill warnings temporarily

* `contextlib.closing` ➔ the same as the example above, call `close`

* parallel programming

- `concurrent.futures.ThreadPoolExecutor` ➔ invoke in parallel then kill thread pool (py >= 3.2)

- `concurrent.futures.ProcessPoolExecutor` ➔ invoke in parallel then kill process pool (py >= 3.2)

- `nogil` ➔ solve the GIL problem temporarily (cython only: ()

### 7.3.1 Catching exceptions

When an exception is thrown in the `with-block`, it is passed as arguments to `__exit__`. Three arguments are used, the same as returned by `sys.exc_info()`: type, value, traceback. When no exception is thrown, `None` is used for all three arguments. The context manager can "swallow" the exception by returning a true value from
__exit__. Exceptions can be easily ignored, because if __exit__ doesn't use return and just falls of the end, None is returned, a false value, and therefore the exception is rethrown after __exit__ is finished.

The ability to catch exceptions opens interesting possibilities. A classic example comes from unit-tests — we want to make sure that some code throws the right kind of exception:

```python
class assert_raises(object):
    # based on pytest and unittest.testCase
    def __init__(self, type):
        self.type = type
    def __enter__(self):
        pass
    def __exit__(self, type, value, traceback):
        if type is None:
            raise AttributeError('exception expected')
        if isinstance(type, self.type):
            return True  # swallow the expected exception
        raise AttributeError('wrong exception type')

with assert_raises(KeyError):
    {}['foo']
```

### 7.3.2 Using generators to define context managers

When discussing *generators*, it was said that we prefer generators to iterators implemented as classes because they are shorter, sweeter, and the state is stored as local, not instance, variables. On the other hand, as described in *Bidirectional communication*, the flow of data between the generator and its caller can be bidirectional. This includes exceptions, which can be thrown into the generator. We would like to implement context managers as special generator functions. In fact, the generator protocol was designed to support this use case.

```python
@contextlib(contextmanager
def some_generator(<arguments>):
    <setup>
    try:
        yield <value>
    finally:
        <cleanup>
```

The `contextlib(contextmanager)` helper takes a generator and turns it into a context manager. The generator has to obey some rules which are enforced by the wrapper function — most importantly it must `yield` exactly once. The part before the `yield` is executed from `__enter__`, the block of code protected by the context manager is executed when the generator is suspended in `yield`, and the rest is executed in `__exit__`. If an exception is thrown, the interpreter hands it to the wrapper through `__exit__` arguments, and the wrapper function then throws it at the point of the `yield` statement. Through the use of generators, the context manager is shorter and simpler.

Let's rewrite the `closing` example as a generator:

```python
@contextlib(contextmanager
def closing(obj):
    try:
        yield obj
    finally:
        obj.close()
```

Let's rewrite the `assert_raises` example as a generator:

```python
@contextlib(contextmanager
def assert_raises(type):
    try:
        yield
    except type:
```
```python
return
except Exception as value:
    raise AttributeError('wrong exception type')
else:
    raise AttributeError('exception expected')
```

Here we use a decorator to turn generator functions into context managers!
Advanced Numpy

**Author: Pauli Virtanen**

Numpy is at the base of Python's scientific stack of tools. Its purpose is to implement efficient operations on many items in a block of memory. Understanding how it works in detail helps in making efficient use of its flexibility, taking useful shortcuts.

This section covers:

* Anatomy of Numpy arrays, and its consequences. Tips and tricks.

* Universal functions: what, why, and what to do if you want a new one.

* Integration with other tools: Numpy offers several ways to wrap any data in an ndarray, without unnecessary copies.

* Recently added features, and what's in them: PEP 3118 buffers, generalized ufuncs, ...

**Prerequisites**

* Numpy

* Cython

* Pillow (Python imaging library, used in a couple of examples)
## Chapter contents

* Life of `ndarray`

- It's...

- Block of memory

- Data types

- Indexing scheme: strides

- Findings in dissection

* Universal functions

- What they are?

- Exercise: building an `ufunc` from scratch

- Solution: building an `ufunc` from scratch

- Generalized `ufuncs`

* Interoperability features

- Sharing multidimensional, typed data

- The old buffer protocol

- The old buffer protocol

- Array interface protocol

* Array siblings: `chararray`, `maskedarray`, `matrix`

- `chararray`: vectorized string operations

- `masked_array` missing data

- `recarray`: purely convenience

- `matrix`: convenience?

* Summary

* Contributing to `Numpy/Scipy`

- Why

- Reporting bugs

- Contributing to documentation

- Contributing features

- How to help, in general

In this section, ` numpy ` will be imported as follows:

```python
>>> import numpy as np
```

# 8.1 Life of `ndarray`

## 8.1.1 It's...

`ndarray` =

block of memory + indexing scheme + data type descriptor

* raw data

* how to locate an element

* how to interpret an element
<figure><img src="image_119.png" /></figure>

```c
typedef struct PyArrayObject {
    PyObject_HEAD

    /* Block of memory */
    char *data;

    /* Data type descriptor */
    PyArrayDescr *descr;

    /* Indexing scheme */
    int nd;
    npy_intp *dimensions;
    npy_intp *strides;

    /* Other stuff */
    PyObject *base;
    int flags;
    PyObject *weakreflist;
} PyArrayObject;
```

### 8.1.2 Block of memory

```c
>>> x = np.array([1, 2, 3], dtype= np.int32)
>>> x.data
<... at ...>
>>> str(x.data)
'\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00'
```

Memory address of the data:

```python
>>> x.__array_interface__['data'] [0]
64803824
```

The whole `__array_interface__`:

```python
>>> x.__array_interface__
{'data': (35828928, False),
 'descr': [((), '<i4')],
 'shape': (4,),
 'strides': None,
 'typestr': '<i4',
 'version': 3}
```

Reminder: two ndarrays may share the same memory:

```python
>>> x = np.array([1, 2, 3, 4])
>>> y = x[:-1]
```
```python
>>> x[0] = 9
>>> y
array([9, 2, 3])
```

Memory does not need to be owned by an `ndarray`:

```python
>>> x = b'1234'  # The 'b' is for "bytes", necessary in Python 3
```

x is a string (in Python 3 a `bytes`), we can represent its data as an array of ints:

```python
>>> y = np.frombuffer(x, type(np.int8))
>>> y.data
<... at ...>
>>> y.base is x
True
>>> y.flags
C_CONTINGUOUS : True
F_CONTINGUOUS : True
OWNDATA : False
WRITEABLE : False
ALIGNED : True
UPDATEIFCOPY : False
```

The `owndata` and `writeable` flags indicate status of the memory block.

### 8.1.3 Data types

#### The descriptor

`dtype` describes a single item in the array:

<table><tr><td>type</td><td><strong>scalar type</strong> of the data, one of:<br>int8, int16, float64, et al. (fixed size)<br>str, unicode, void (flexible size)</td></tr><tr><td>itemsize</td><td><strong>size</strong> of the data block</td></tr><tr><td>byte-<br>order</td><td><strong>byte order</strong>: big-endian &gt; / little-endian &lt; / not applicable |</td></tr><tr><td>fields</td><td>sub-dtypes, if it's a <strong>structured data type</strong></td></tr><tr><td>shape</td><td>shape of the array, if it's a <strong>sub-array</strong></td></tr></table>

```python
>>> np.dtype(int).type
<type 'numpy.int64'>
>>> np.dtype(int).itemsize
8
>>> np.dtype(int).byteorder
'='
```

#### Example: reading `.wav` files

The `.wav` file header:
<table><tr><td>chunk_id</td><td>"RIFF"</td></tr><tr><td>chunk_size</td><td>4-byte unsigned little-endian integer</td></tr><tr><td>format</td><td>"WAVE"</td></tr><tr><td>fmt_id</td><td>"fmt "</td></tr><tr><td>fmt_size</td><td>4-byte unsigned little-endian integer</td></tr><tr><td>audio fmt</td><td>2-byte unsigned little-endian integer</td></tr><tr><td>num channels</td><td>2-byte unsigned little-endian integer</td></tr><tr><td>sample rate</td><td>4-byte unsigned little-endian integer</td></tr><tr><td>byte rate</td><td>4-byte unsigned little-endian integer</td></tr><tr><td>block align</td><td>2-byte unsigned little-endian integer</td></tr><tr><td>bits per sample</td><td>2-byte unsigned little-endian integer</td></tr><tr><td>data id</td><td>"data"</td></tr><tr><td>data size</td><td>4-byte unsigned little-endian integer</td></tr></table>

* 44-byte block of raw data (in the beginning of the file)

* ... followed by `data_size` bytes of actual sound data.

The `.wav` file header as a Numpy *structured* data type:

```python
>>> wav_header_dtype = npลง( [
...     ("chunk_id", (bytes, 4)), # flexible-sized scalar type, item size 4
...     ("chunk_size", "<u4"),    # little-endian unsigned 32-bit integer
...     ("format", "S4"),        # 4-byte string
...     ("fmt_id", "S4"),
...     ("fmt_size", "<u4"),
...     ("audio fmt", "<u2"),   #
...     ("num channels", "<u2"), # .. more of the same ...
...     ("sample rate", "<u4"), #
...     ("byte rate", "<u4"),
...     ("block align", "<u2"),
...     ("bits per sample", "<u2"),
...     ("data_id", ("S1", (2, 2))), # sub-array, just for fun!
...     ("data_size", "u4"),
...     #
...     # the sound data itself cannot be represented here:
...     # it does not have a fixed size
...   ])
```

**See also:**

wavreader.py

```python
>>> wav_header_dtype['format']
dtype('S4')
>>> wav_header_dtypefields
dict_proxy({'blockalign': (dtype('uint16'), 32), 'format': (dtype('S4'), 8), 'data_id': (dtype(( 'S1', (2, 2)))
>>> wav_header_dtypefields['format']
dtype('S4'), 8)
```

* The first element is the sub-dtype in the structured data, corresponding to the name `format`

* The second one is its offset (in bytes) from the beginning of the item
## Exercise

Mini-exercise, make a "sparse" dtype by using offsets, and only some of the fields:

```python
>>> wav_header_dtype = np.dtype(dict(
...     names=['format', 'sample_rate', 'data_id'],
...     offsets=[offset_1, offset_2, offset_3],  # counted from start of structure in bytes
...     formats=list of dtypes for each of the fields,
... ))
```

and use that to read the sample rate, and data_id (as sub-array).

```python
>>> f = open('data/test.wav', 'r')
>>> wav_header = np.fromfile(f, dtype=wav_header_dtype, count=1)
>>> f.close()
>>> print(wav_header)
[ ('RIFF', 17402L, 'WAVE', 'fmt ', 16L, 1, 1, 16000L, 32000L, 2, 16, [[ 'd', 'a'], ['t', 'a']], 17366L)]
>>> wav_header['sample_rate']
array([16000], dtype=uint32)
```

Let's try accessing the sub-array:

```python
>>> wav_header['data_id']
array([[[ 'd', 'a'],
      ['t', 'a']]]],
     dtype='|S1')
>>> wav_header.shape
(1,)
>>> wav_header['data_id'] .shape
(1, 2, 2)
```

When accessing sub-arrays, the dimensions get added to the end!

There are existing modules such as wavfile, audiolab, etc. for loading sound data...

## Casting and re-interpretation/views

### casting

* on assignment

* on array construction

* on arithmetic

* etc.

* and manually: `.astype(dtype)`

### data re-interpretation

* manually: `.view(dtype)`

### Casting

* Casting in arithmetic, in nutshell:

- only type (not value!) of operands matters

- largest "safe" type able to represent both is picked

- scalars can "lose" to arrays in some situations

* Casting in general copies data:
```python
>>> x = np.array([1, 2, 3, 4], dtype= np.float)
>>> x
array([1., 2., 3., 4.])
>>> y = xerve(np.int8)
>>> y
array([1, 2, 3, 4], ndtype=int8)
>>> y + 1
array([2, 3, 4, 5], ndtype=int8)
>>> y + 256
array([257, 258, 259, 260], ndtype=int16)
>>> y + 256.0
array([ 257.,  258.,  259.,  260.])
>>> y + np.array([256], ndtype= np.int32)
array([257, 258, 259, 260], ndtype=int32)
```

* Casting on setitem: ndtype of the array is not changed on item assignment:

```python
>>> y[:] = y + 1.5
>>> y
array([2, 3, 4, 5], ndtype=int8)
```

Exact rules: see ` numpy documentation`

## Re-interpretation / viewing

* Data block in memory (4 bytes)

<table><tr><td>0x01</td><td>||</td><td>0x02</td><td>||</td><td>0x03</td><td>||</td><td>0x04</td></tr></table>

- 4 of `uint8`, OR,

- 4 of `int8`, OR,

- 2 of `int16`, OR,

- 1 of `int32`, OR,

- 1 of `float32`, OR,

- ...

How to switch from one to another?

1. Switch the `dtype`:

```python
>>> x = np.array([1, 2, 3, 4], ndtype= npinx)
>>> x.dtype = "<i2"
>>> x
array([ 513, 1027], ndtype=int16)
>>> 0x0201, 0x0403
(513, 1027)
```

<table><tr><td>0x01</td><td>0x02</td><td>||</td><td>0x03</td><td>0x04</td></tr></table>

little-endian: least significant byte is on the left in memory

2. Create a new view:

```python
>>> y = x.view("<i4")
>>> y
array([67305985], ndtype=int32)
>>> 0x04030201
67305985
```

<table><tr><td>0x01</td><td>0x02</td><td>0x03</td><td>0x04</td></tr></table>
* `.view()` makes *views*, does not copy (or alter) the memory block

* only changes the `dtype` (and adjusts array shape):

```python
>>> x[1] = 5
>>> y
array([328193], ndtype=int32)
>>> y.base is x
True
```

## Mini-exercise: data re-interpretation

**See also:**

view-colors.py

You have `RGBA` data in an array:

```python
>>> x = np.zeros((10, 10, 4), ndtype(np.int8)
>>> x[: : , 0] = 1
>>> x[: : , 1] = 2
>>> x[: : , 2] = 3
>>> x[: : , 3] = 4
```

where the last three dimensions are the R, B, and G, and alpha channels.

How to make a (10, 10) structured array with field names 'r', 'g', 'b', 'a' without copying data?

```python
>>> y = ...
>>> assert (y['r'] == 1).all()
>>> assert (y['g'] == 2).all()
>>> assert (y['b'] == 3).all()
>>> assert (y['a'] == 4).all()
```

**Solution**

```python
>>> y = x.view((('r', 'i1'),
...             ('g', 'i1'),
...             ('b', 'i1'),
...             ('a', 'i1'))]
>>> )[: : , 0]
```
<figure><img src="image_120.png" /></figure>

Another array taking exactly 4 bytes of memory:

```python
>>> y = np.array([[1, 3], [2, 4]], dtype= nplass)
>>> x = y.copy()
>>> x
array([[1, 2],
      [3, 4]], ndtype= npint)
>>> y
array([[1, 2],
      [3, 4]], ndtype= npint)
>>> x.view(np.int16)
array([[ 513],
      [1027]], ndtype= npint)
>>> 0x0201, 0x0403
(513, 1027)
>>> y.view(np.int16)
array([[ 769, 1026]], ndtype= npint)
```

* What happened?

* ... we need to look into what x[0, 1] actually means

```python
>>> 0x0301, 0x0402
(769, 1026)
```

### 8.1.4 Indexing scheme: strides

#### Main point

##### The question:

```python
>>> x = np.array([[1, 2, 3],
...                [4, 5, 6],
...                [7, 8, 9]], ndtype= np.int8)
>>> str(x.data)
'\x01\x02\x03\x04\x05\x06\x07\x08\t'
```
At which byte in ``x.data`` does the item ``x[1, 2]`` begin?

##### The answer (in Numpy)

* strides: the number of bytes to jump to find the next element

* 1 stride per dimension

```python
>>> x.strides
(3, 1)
>>> byte_offset = 3*1 + 1*2  # to find x[1, 2]
>>> x flat[byte_offset]
6
>>> x[1, 2]
6
- simple, **flexible**
```

#### C and Fortran order

```python
>>> x = np.array([[1, 2, 3],
...                [4, 5, 6]], ndtype= np.int16, order='C')
>>> x.strides
(6, 2)
```
```python
>>> str(x.data)
'\x01\x00\x02\x00\x03\x00\x04\x00\x05\x00\x06\x00'
```

* Need to jump 6 bytes to find the next row

* Need to jump 2 bytes to find the next column

```python
>>> y = np.array(x, order='F')
>>> y.strides
(2, 4)
>>> str(y.data)
'\x01\x00\x04\x00\x02\x00\x05\x00\x03\x00\x06\x00'
```

* Need to jump 2 bytes to find the next row

* Need to jump 4 bytes to find the next column

* Similarly to higher dimensions:

- C: last dimensions vary fastest (= smaller strides)

- F: first dimensions vary fastest

```
shape = (d1, d2, ..., d_n)
```

```
strides = (s1, s2, ..., s_n)
```

```
s_j^C = d_{j+1}d_{j+2}...d_n \times \text{itemsize}
```

```
s_j^F = d_1d_2...d_{j-1} \times \text{itemsize}
```

Now we can understand the behavior of `.view()`:

```python
>>> y = np.array([[1, 3], [2, 4]], dtype(npмат.uint8)).transpose()
>>> x = y.copy()
```

Transposition does not affect the memory layout of the data, only strides

```python
>>> x.strides
(2, 1)
>>> y.strides
(1, 2)
```

```python
>>> str(x.data)
'\x01\x02\x03\x04'
>>> str(y.data)
'\x01\x03\x02\x04'
```

* the results are different when interpreted as 2 of `int16`

* `.copy()` creates new arrays in the C order (by default)

## Slicing with integers

* Everything can be represented by changing only `shape`, `strides`, and possibly adjusting the data pointer!

* Never makes copies of the data

```python
>>> x = np.array([1, 2, 3, 4, 5, 6], dtype(np.іnt32)
>>> y = x[::-1]
>>> y
array([6, 5, 4, 3, 2, 1], dtype=int32)
>>> y.strides
(-4,)
```
```python
>>> y = x[2:]
>>> y.__array_interface__['data'] [0] - x.__array_interface__['data'] [0]
8
>>> x = np(zeros((10, 10, 10), dtype(np.float))
>>> x.strides
(800, 80, 8)
>>> x[: :2, : :3, : :4].strides
(1600, 240, 32)
```

* Similarly, transposes never make copies (it just swaps strides):

```python
>>> x = np(zeros((10, 10, 10), dtype(np.float))
>>> x.strides
(800, 80, 8)
>>> x.T.strides
(8, 80, 800)
```

But: not all reshaping operations can be represented by playing with strides:

```python
>>> a = np.arange(6, dtype(np.int8)).reshape(3, 2)
>>> b = a.T
>>> b.strides
(1, 2)
```

So far, so good. However:

```python
>>> str(a.data)
'\x00\x01\x02\x03\x04\x05'
>>> b
array([[0, 2, 4],
      [1, 3, 5]], dtype=int8)
>>> c = b reshape(3*2)
>>> c
array([0, 2, 4, 1, 3, 5], dtype=int8)
```

Here, there is no way to represent the array `c` given one stride and the block of memory for `a`. Therefore, the `reshape` operation needs to make a copy here.

## Example: fake dimensions with strides

### Stride manipulation

```python
>>> from numpy.lib.stride_tricks import as_strided
>>> help(as_strided)
as_strided(x, shape=None, strides=None)
Make an ndarray from the given array with the given shape and strides
```

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

`as_strided` does **not** check that you stay inside the memory block bounds...

```python
>>> x = np.array([1, 2, 3, 4], dtype(np.int16)
>>> as_strided(x, strides=(2*2, ), shape=(2, ))
array([1, 3], dtype=int16)
>>> x[: :2]
array([1, 3], dtype=int16)
```

### See also:

`stride-fakedims.py`

### Exercise
```python
array([1, 2, 3, 4], dtype(np.int8))
-> array([[1, 2, 3, 4],
         [1, 2, 3, 4],
         [1, 2, 3, 4]], dtype(np.int8))
```

using only `as_strided`:

```
Hint: byte_offset = stride[0]*index[0] + stride[1]*index[1] + ...
```

### Spoiler

Stride can also be 0:

```python
>>> x = np.array([1, 2, 3, 4], dtype(np.int8))
>>> y = as_strided(x, strides=(0, 1), shape=(3, 4))
>>> y
array([[1, 2, 3, 4],
      [1, 2, 3, 4],
      [1, 2, 3, 4]], dtype=int8)
>>> y.base.base is x
True
```

## Broadcasting

* Doing something useful with it: outer product of [1, 2, 3, 4] and [5, 6, 7]

```python
>>> x = np.array([1, 2, 3, 4], dtype(np.int16)
>>> x2 = as_strided(x, strides=(0, 1*2), shape=(3, 4))
>>> x2
array([[1, 2, 3, 4],
      [1, 2, 3, 4],
      [1, 2, 3, 4]], dtype=int16)
```

```python
>>> y = np.array([5, 6, 7], dtype(np.int16)
>>> y2 = as_strided(y, strides=(1*2, 0), shape=(3, 4))
>>> y2
array([[5, 5, 5, 5],
      [6, 6, 6, 6],
      [7, 7, 7, 7]], dtype=int16)
```

```python
>>> x2 * y2
array([[ 5, 10, 15, 20],
      [ 6, 12, 18, 24],
      [ 7, 14, 21, 28]], dtype=int16)
```

... seems somehow familiar ...

```python
>>> x = np.array([1, 2, 3, 4], dtype(np.int16)
>>> y = np.array([5, 6, 7], dtype(np.int16)
>>> x[np.newaxis,:]*y[,np.newaxis]
array([[ 5, 10, 15, 20],
      [ 6, 12, 18, 24],
      [ 7, 14, 21, 28]], dtype=int16)
```

* Internally, array **broadcasting** is indeed implemented using 0-strides.

## More tricks: diagonals

### See also:
stride-diagonals.py

### Challenge

* Pick diagonal entries of the matrix: (assume C memory order):

```python
>>> x = np.array([[1, 2, 3],
...                [4, 5, 6],
...                [7, 8, 9]], dtype= np.int32)
>>> x_diag = as_strided(x, shape=(3,), strides=(???,))
```

* Pick the first super-diagonal entries [2, 6].

* And the sub-diagonals?

(Hint to the last two: slicing first moves the point where striding starts from.)

### Solution

Pick diagonals:

```python
>>> x_diag = as_strided(x, shape=(3, ), strides=((3+1)*x.itemsize, ))
>>> x_diag
array([1, 5, 9], type=int32)
```

Slice first, to adjust the data pointer:

```python
>>> as_strided(x[0, 1:], shape=(2, ), strides=((3+1)*x.itemsize, ))
array([2, 6], type=int32)
>>> as_strided(x[1:, 0], shape=(2, ), strides=((3+1)*x.itemsize, ))
array([4, 8], type=int32)
```

### Using np.diag

```python
>>> y = np.diag(x, k=1)
>>> y
array([2, 6], type=int32)
```

However,

```python
>>> y_flags.owndata
False
```

**Note** This behavior has changed: before numpy 1.9, np.diag would make a copy.

### See also:

stride-diagonals.py

### Challenge

Compute the tensor trace:

```python
>>> x = np.arange(5*5*5*5).reshape(5, 5, 5, 5)
>>> s = 0
>>> for i in range(5):
...     for j in range(5):
...         s += x[j, i, j, i]
```

by striding, and using `sum()` on the result.

```python
>>> y = as_strided(x, shape=(5, 5), strides=(TODO, TODO))
>>> s2 = ...
>>> assert s == s2
```

### Solution
```python
>>> y = as_strided(x, shape=(5, 5), strides=((5*5*5 + 5)*x.itemsize,
...                                         (5*5 + 1)*x.itemsize))
>>> s2 = y.sum()
```

## CPU cache effects

Memory layout can affect performance:

```python
In [1]: x = np.zeros((20000,))

In [2]: y = np.zeros((20000*67,))[::67]

In [3]: x.shape, y.shape
((20000,), (20000,))

In [4]: %timeit x.sum()
100000 loops, best of 3: 0.180 ms per loop

In [5]: %timeit y.sum()
100000 loops, best of 3: 2.34 ms per loop

In [6]: x.strides, y.strides
((8,), (536,))
```

### Smaller strides are faster?

x

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

y

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

### cache block size

* CPU pulls data from main memory to its cache in blocks

* If many array items consecutively operated on fit in a single block (small stride):

- ⇒ fewer transfers needed

- ⇒ faster

**See also:**

`numexpr` is designed to mitigate cache effects in array computing.

### Example: inplace operations (caveat emptor)

* Sometimes,

```python
>>> a -= b
```

is not the same as

```python
>>> a -= b.copy()
```
```python
>>> x = np.array([[1, 2], [3, 4]])
>>> x = x transpose()
>>> x
array([[ 0, -1],
      [ 4,  0]])
```

```python
>>> y = np.array([[1, 2], [3, 4]])
>>> y = y.T.copy()
>>> y
array([[ 0, -1],
      [ 1,  0]])
```

* x and `x transpose()` share data

* `x = x transpose()` modifies the data element-by-element...

* because x and `x transpose()` have different striding, modified data re-appears on the RHS

### 8.1.5 Findings in dissection

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

* *memory block*: may be shared, `.base`, `.data`

* *data type descriptor*: structured data, sub-arrays, byte order, casting, viewing, `.astype()`, `.view()`

* *strided indexing*: strides, C/F-order, slicing w/ integers, `as_strided`, broadcasting, stride tricks, diag, CPU cache coherence

## 8.2 Universal functions

### 8.2.1 What they are?

* Ufunc performs and elementwise operation on all elements of an array.

Examples:

`np.add`, `np.subtract`, `scipy.special.*, ...`

* Automatically support: broadcasting, casting, ...

* The author of an ufunc only has to supply the elementwise operation, Numpy takes care of the rest.

* The elementwise operation needs to be implemented in C (or, e.g., Cython)
## Parts of an Ufunc

1. Provided by user

```c
void ufunc_loop(void **args, int *dimensions, int *steps, void *data)
{
    /*
     * int8 output = elementwise_function(int8 input_1, int8 input_2)
     *
     * This function must compute the ufunc for many values at once,
     * in the way shown below.
     */
    char *input_1 = (char*)args[0];
    char *input_2 = (char*)args[1];
    char *output = (char*)args[2];
    int i;

    for (i = 0; i < dimensions[0]; ++i) {
        *output = elementwise_function(*input_1, *input_2);
        input_1 += steps[0];
        input_2 += steps[1];
        output += steps[2];
    }
}
```

2. The Numpy part, built by

```c
char types[3]

types[0] = NPY_BYTES   /* type of first input arg */
types[1] = NPY_BYTES   /* type of second input arg */
types[2] = NPY_BYTES   /* type of third input arg */

PyObject *python_ufunc = PyUFunc_FromFuncAndData(
    ufunc_loop,
    NULL,
    types,
    1, /* ntypes */
    2, /* numInputs */
    1, /* numOutputs */
    identity_element,
    name,
    docstring,
    unused)
```

* A ufunc can also support multiple different input-output type combinations.

## Making it easier

3. `ufunc_loop` is of very generic form, and Numpy provides pre-made ones

<table><tr><td>PyUFunc_f_f</td><td>float elementwise_func(float input_1)</td></tr><tr><td>PyUFunc_ff_f</td><td>float elementwise_func(float input_1, float input_2)</td></tr><tr><td>PyUFunc_d_d</td><td>double elementwise_func(double input_1)</td></tr><tr><td>PyUFunc_dd_d</td><td>double elementwise_func(double input_1, double input_2)</td></tr><tr><td>PyUFunc_D_D</td><td>elementwise_func(npy_cdouble *input, npy_cdouble* output)</td></tr><tr><td>PyUFunc_DD_D</td><td>elementwise_func(npy_cdouble *in1, npy_cdouble *in2, npy_cdouble* out)</td></tr></table>

* Only `elementwise_func` needs to be supplied

* ... except when your elementwise function is not in one of the above forms
### 8.2.2 Exercise: building an ufunc from scratch

The Mandelbrot fractal is defined by the iteration

$$z \leftarrow z^2 + c$$

where $c = x + iy$ is a complex number. This iteration is repeated – if $z$ stays finite no matter how long the iteration runs, $c$ belongs to the Mandelbrot set.

* Make ufunc called mandel(z0, c) that computes:

```python
z = z0
for k in range(iterations):
    z = z * z + c
```

say, 100 iterations or until $z$.real**2 + z.imag**2 > 1000$. Use it to determine which $c$ are in the Mandelbrot set.

* Our function is a simple one, so make use of the PyUFunc_* helpers.

* Write it in Cython

**See also:**

mandel.pyx, mandelplot.pyx

```c
#
# Fix the parts marked by TODO
#
#
#
# Compile this file by (Cython >= 0.12 required because of the complex vars)
#
#   cython mandel.pyx
#   python setup.py build_ext -i
#
# and try it out with, in this directory,
#
# >>> import mandel
# >>> mandel.mandel(0, 1 + 2j)
#
#
```

```c
# The elementwise function
# ------------------------
```

```c
cdef void mandel_single_point(double complex *z_in,
                              double complex *c_in,
                              double complex *z_out) nogil:

# The Mandelbrot iteration
#
#
# Some points of note:
#
# - It's *NOT* allowed to call any Python functions here.
#
#   The Ufunc loop runs with the Python Global Interpreter Lock released.
#   Hence, the ``nogil``.
#
# - And so all local variables must be declared with ``cdef``
#
# - Note also that this function receives *pointers* to the data
#
```
```c
cdef double complex z = z_in[0]
cdef double complex c = c_in[0]
cdef int k  # the integer we use in the for loop

#
# TODO: write the Mandelbrot iteration for one point here,
#       as you would write it in Python.

#
# Say, use 100 as the maximum number of iterations, and 1000
# as the cutoff for z.real**2 + z.imag**2.

#

TODO: mandelbrot iteration should go here

# Return the answer for this point
z_out[0] = z

# Boilerplate Cython definitions
#
# The litany below is particularly long, but you don't really need to
# read this part; it just pulls in stuff from the Numpy C headers.
#

cdef extern from "numpy/arrayobject.h":
    void import_array()
    ctypedef int npy_intp
    cdef enum NPY_types:
        NPY_DOUBLE
        NPY_CDOUBLE
        NPY_LONG

cdef extern from "numpy/ufuncobject.h":
    void import_ufunc()
    ctypedef void (*PyUFuncGenericFunction)(char**, npy_intp*, npy_intp*, void*)
    object PyUFunc_FromFuncAndData(PyUFuncGenericFunction* func, void** data,
                                  char* types, int ntypes, int nin, int nout,
                                  int identity, char* name, char* doc, int c)

# List of pre-defined loop functions

void PyUFunc_f_f_As_d_d(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_d_d(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_f_f(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_g_g(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_F_F_As_D_D(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_F_F(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_D_D(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_G_G(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_ff_f_As_dd_d(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_ff_f(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_dd_d(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_gg_g(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFuncFFE_F_As_DD_D(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_DD_D(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFuncFFE_F(char** args, npy_intp* dimensions, npy_intp* steps, void* func)
void PyUFunc_GG_G(char** args, npy_intp* dimensions, npy_intp* steps, void* func)

# Required module initialization
#
import_array()
```
```c
import_ufunc()

# The actual ufunc declaration
# --------------------------

cdef PyUFuncGenericFunction loop_func[1]
cdef char input_output_types[3]
cdef void *elementwise_funcs[1]

# 
# Reminder: some pre-made Ufunc loops:
# 
# ================= ====================================================
# ``PyUFunc_f_f`` ``float elementwise_func(float input_1)``
# ``PyUFunc_ff_f`` ``float elementwise_func(float input_1, float input_2)``
# ``PyUFunc_d_d`` ``double elementwise_func(double input_1)``
# ``PyUFunc_dd_d`` ``double elementwise_func(double input_1, double input_2)``
# ``PyUFunc_D_D`` ``elementwise_func(complex_double *input, complex_double* complex_double)``
# ``PyUFunc_DD_D`` ``elementwise_func(complex_double *in1, complex_double *in2, complex_double* out)``
# ===================
# 
# The full list is above.
# 
# 
# Type codes:
# 
# NPY Бог, NPY_BYTES, NPY_UBYTES, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UUID,
# NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG, NPY_FLOAT, NPY_DOUBLE,
# NPY_LONGDOUBLE, NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE, NPY_DATETIME,
# NPY_TIMEDELTA, NPY_OBJECT, NPY_STRING, NPY_UNICODE, NPY_VOID
# 
loop_func[0] = ... TODO: suitable PyUFunc_* ...
input_output_types[0] = ... TODO ...
... TODO: fill in rest of input_output_types ...

# This thing is passed as the ``data`` parameter for the generic
# PyUFunc_* loop, to let it know which function it should call.
elementwise_funcs[0] = <void*>mandel_single_point

# Construct the ufunc:
mandel = PyUFunc_FromFuncAndData(
    loop_func,
    elementwise_funcs,
    input_output_types,
    1, # number of supported input types
    TODO, # number of input args
    TODO, # number of output args
    0, # `identity` element, never mind this
    "mandel", # function name
    "mandel(z, c) -> computes z*z + c", # docstring
    0 # unused
)
```

Reminder: some pre-made Ufunc loops:
<table><tr><td>PyUfunc_f_f</td><td>float elementwise_func(float input_1)</td></tr><tr><td>PyUfunc_ff_f</td><td>float elementwise_func(float input_1, float input_2)</td></tr><tr><td>PyUfunc_d_d</td><td>double elementwise_func(double input_1)</td></tr><tr><td>PyUfunc_dd_d</td><td>double elementwise_func(double input_1, double input_2)</td></tr><tr><td>PyUfunc_D_D</td><td>elementwise_func(complex_double *input, complex_double* output)</td></tr><tr><td>PyUfunc_DD_D</td><td>elementwise_func(complex_double *in1, complex_double *in2, complex_double* out)</td></tr></table>

Type codes:

NPY Boolean, NPY BYTE, NPY_UBYTE, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UUID,
NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG, NPY_FLOAT, NPY_DOUBLE,
NPY_LONGDOUBLE, NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE, NPY_DATETIME,
NPY_TIMEDELTA, NPY_OBJECT, NPY_STRING, NPY_UNICODE, NPY_VOID

### 8.2.3 Solution: building an ufunc from scratch

# The elementwise function

# -----------------------

```c
cdef void mandel_single_point(double complex *z_in,
                              double complex *c_in,
                              double complex *z_out) nogil:

# The Mandelbrot iteration
#

# Some points of note:
#
# - It's *NOT* allowed to call any Python functions here.
#
# The Ufunc loop runs with the Python Global Interpreter Lock released.
# Hence, the ``nogil``.
#
# - And so all local variables must be declared with ``cdef``
#
# - Note also that this function receives *pointers* to the data;
#   the "traditional" solution to passing complex variables around
#
cdef double complex z = z_in[0]
cdef double complex c = c_in[0]
cdef int k  # the integer we use in the for loop

# Straightforward iteration

for k in range(100):
    z = z*z + c
    if z.real**2 + z.imag**2 > 1000:
        break

# Return the answer for this point
z_out[0] = z
```

# Boilerplate Cython definitions

#
# You don't really need to read this part, it just pulls in
# stuff from the Numpy C headers.

# ----------------------------------------------------------------------
```c
cdef extern from "numpy/arrayobject.h":
    void import_array()
    ctypedef int npy_intp
cdef enum NPY_types:
    NPY_CDOUBLE

cdef extern from "numpy/ufuncobject.h":
    void import_udf()
    ctypedef void (*PyUDFGenericFunction)(char**, npy_intp*, npy_intp*, void*)
    object PyUDFFromFuncAndData(PyUDFGenericFunction* func, void** data,
                              char* types, int ntypes, int nin, int nout,
                              int identity, char* name, char* doc, int c)
    void PyUDF_DD_D(char**, npy_intp*, npy_intp*, void*)
```

```c
# Required module initialization
# -----------------------------
```

```python
import_array()
import_udf()
```

```c
# The actual ufunc declaration
# -----------------------------
```

```c
cdef PyUDFGenericFunction loop_func[1]
cdef char input_output_types[3]
cdef void *elementwiseUDFs[1]

loop_func[0] = PyUDF_DD_D

input_output_types[0] = NPY_CDOUBLE
input_output_types[1] = NPY_CDOUBLE
input_output_types[2] = NPY_CDOUBLE

elementwiseUDFs[0] = <void*>mandel_single_point

mandel = PyUDF_FromFuncAndData(
    loop_func,
    elementwiseUDFs,
    input_output_types,
    1, # number of supported input types
    2, # number of input args
    1, # number of output args
    0, # `identity` element, never mind this
    "mandel", # function name
    "mandel(z, c) -> computes iterated z*z + c", # docstring
    0 # unused
)
```

```python
#"
Plot Mandelbrot
====================
```

Plot the Mandelbrot ensemble.

```python
#"
import numpy as np
import mandel
x = np.linspace(-1.7, 0.6, 1000)
y = np.linspace(-1.4, 1.4, 1000)
```
```python
c = x[None] + 1j*y[:3, None]
z = mandel.mandel(c, c)

import matplotlib.pyplot as plt
plt.imshow(abs(z)**2 < 1000, extent=[-1.7, 0.6, -1.4, 1.4])
plt gray()
plt show()
```

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

Most of the boilerplate could be automated by these Cython modules:
http://wiki.cython.org/MarkLodato/CreatingUfunc

## Several accepted input types

E.g. supporting both single- and double-precision versions

```c
cdef void mandel_single_point(double complex *z_in,
                              double complex *c_in,
                              double complex *z_out) nogil:

...

cdef void mandel_single_point_singleprec(float complex *z_in,
                                         float complex *c_in,
                                         float complex *z_out) nogil:

...

cdef PyUFuncGenericFunction loop_funcs[2]
cdef char input_output_types[3*2]
cdef void *elementwise_funcs[1*2]

loop_funcs[0] = PyUFunc_DD_D
input_output_types[0] = NPY_CDOUBLE
input_output_types[1] = NPY_CDOUBLE
input_output_types[2] = NPY_CDOUBLE
elementwise_funcs[0] = <void*>mandel_single_point

loop_funcs[1] = PyUFunc FF_F
input_output_types[3] = NPY_CFLOAT
input_output_types[4] = NPY_CFLOAT
input_output_types[5] = NPY_CFLOAT
elementwise_funcs[1] = <void*>mandel_single_point_singleprec

mandel = PyUFunc_FromFuncAndData(
    loop_func,
    elementwise_funcs,
    input_output_types,
```
2, # number of supported input types <------------------
2, # number of input args
1, # number of output args
0, # `identity` element, never mind this
"mandel", # function name
"mandel(z, c) -> computes iterated z*z + c", # docstring
0 # unused
)

### 8.2.4 Generalized ufunc

#### ufunc

output = elementwise_function(input)

Both output and input can be a single array element only.

#### generalized ufunc

output and input can be arrays with a fixed number of dimensions

For example, matrix trace (sum of diag elements):

input shape = (n, n)
output shape = ()      i.e. scalar

(n, n) -> ()

Matrix product:

input_1 shape = (m, n)
input_2 shape = (n, p)
output shape = (m, p)

(m, n), (n, p) -> (m, p)

* This is called the “signature” of the generalized ufunc

* The dimensions on which the g-ufunc acts, are “core dimensions”

#### Status in Numpy

* g-ufuncs are in Numpy already ...

* new ones can be created with PyUFunc_FromFuncAndDataAndSignature

* ... but we don't ship with public g-ufuncs, except for testing, ATM

```python
>>> import numpy.core.umath tests as ut
>>> ut.matrix_multiplier signature
'(m,n),(n,p)->(m,p)'
```

```python
>>> x = np.ones((10, 2, 4))
>>> y = np.ones((10, 4, 5))
>>> ut.matrix_multiplier(x, y).shape
(10, 2, 5)
```

* the last two dimensions became core dimensions, and are modified as per the signature

* otherwise, the g-ufunc operates “elementwise”

* matrix multiplication this way could be useful for operating on many small matrices at once
## Generalized ufunc loop

Matrix multiplication (m,n), (n,p) -> (m,p)

```c
void gufunc_loop(void **args, int *dimensions, int *steps, void *data)
{
    char *input_1 = (char*)args[0]; /* these are as previously */
    char *input_2 = (char*)args[1];
    char *output = (char*)args[2];

    int input_1_stride_m = steps[3]; /* strides for the core dimensions */
    int input_1_stride_n = steps[4]; /* are added after the non-core */
    int input_2_strides_n = steps[5]; /* steps */
    int input_2_strides_p = steps[6];
    int output_strides_n = steps[7];
    int output_strides_p = steps[8];

    int m = dimension[1]; /* core dimensions are added after */
    int n = dimension[2]; /* the main dimension; order as in */
    int p = dimension[3]; /* signature */

    int i;

    for (i = 0; i < dimensions[0]; ++i) {
        matmul_for_strideduner matrices(input_1, input_2, output,
                                          strides for each array...);

        input_1 += steps[0];
        input_2 += steps[1];
        output += steps[2];
    }
}
```

# 8.3 Interoperability features

## 8.3.1 Sharing multidimensional, typed data

Suppose you

1. Write a library that handles (multidimensional) binary data,

2. Want to make it easy to manipulate the data with Numpy, or whatever other library,

3. ... but would **not** like to have Numpy as a dependency.

Currently, 3 solutions:

1. the "old" buffer interface

2. the array interface

3. the "new" buffer interface (PEP 3118)

## 8.3.2 The old buffer protocol

* Only 1-D buffers

* No data type information

* C-level interface; PyBufferProcs tp_as_buffer in the type object

* But it's integrated into Python (e.g. strings support it)
Mini-exercise using `Pillow` (Python Imaging Library):

**See also:**

`pilbuffer.py`

```python
>>> from PIL import Image
>>> data = np.zeros((200, 200, 4), dtype(np.int8)
>>> data[: ,:] = [255, 0, 0, 255] # Red
>>> # In PIL, RGBA images consist of 32-bit integers whose bytes are [RR,GG,BB,AA]
>>> data = data.view(np.int32).squeeze()
>>> img = Image.frombuffer("RGBA", (200, 200), data, "raw", "RGBA", 0, 1)
>>> img.save('test.png')
```

**Q:**

Check what happens if `data` is now modified, and `img` saved again.

### 8.3.3 The old buffer protocol

```python
"""
From buffer
============

Show how to exchange data between numpy and a library that only knows
the buffer interface.
"""

import numpy as np
import Image

# Let's make a sample image, RGBA format
x = np.zeros((200, 200, 4), dtype(np.int8)

x[: : ,0] = 254 # red
x[: : ,3] = 255 # opaque

data = x.view(np.int32) # Check that you understand why this is OK!

img = Image.frombuffer("RGBA", (200, 200), data)
img.save('test.png')

#
# Modify the original data, and save again.
#
# It turns out that PIL, which knows next to nothing about Numpy,
# happily shares the same data.
#

x[: : ,1] = 254
img.save('test2.png')
```
<figure><img src="image_126.png" /></figure>

### 8.3.4 Array interface protocol

* Multidimensional buffers

* Data type information present

* Numpy-specific approach; slowly deprecated (but not going away)

* Not integrated in Python otherwise

**See also:**

Documentation: http://docs.scipy.org/doc/numpy/reference/arrays_interface.html

```python
>>> x = np.array([[1, 2], [3, 4]])
>>> x.__array_interface__
{'data': (171694552, False),          # memory address of data, is readonly?
 'descr': [('<i4')],                   # data type descriptor
 'typestr': '<i4',                    # same, in another form
 'strides': None,                    # strides; or None if in C-order
 'shape': (2, 2),
 'version': 3,
}
```

```python
::
```
## The mask

You can modify the mask by assigning:

```python
>>> mx[1] = np.ma masked
>>> mx
masked_array(data = [1 -- 3 -- 5],
               mask = [False True False True False],
               fill_value = 999999)
```

The mask is cleared on assignment:

```python
>>> mx[1] = 9
>>> mx
masked_array(data = [1 9 3 -- 5],
               mask = [False False False True False],
               fill_value = 999999)
```

The mask is also available directly:

```python
>>> mx mask
array([False, False, False, True, False], dtype=bool)
```

The masked entries can be filled with a given value to get an usual array back:

```python
>>> x2 = mx filled(-1)
>>> x2
array([ 1,  9,  3, -1,  5])
```

The mask can also be cleared:

```python
>>> mx mask = np.ma nomask
>>> mx
masked_array(data = [1 9 3 -99 5],
               mask = [False False False False False],
               fill_value = 999999)
```

## Domain-aware functions

The masked array package also contains domain-aware functions:

```python
>>> np.ma.log(np.array([1, 2, -1, -2, 3, -5]))
masked_array(data = [0.0 0.6931471805599453 -- -- 1.0986122886681098 --],
               mask = [False False True True False True],
               fill_value = 1e+20)
```

Streamlined and more seamless support for dealing with missing data in arrays is making its way into Numpy 1.7. Stay tuned!
## Example: Masked statistics

Canadian rangers were distracted when counting hares and lynxes in 1903-1910 and 1917-1918, and got the numbers are wrong. (Carrot farmers stayed alert, though.) Compute the mean populations over time, ignoring the invalid numbers.

```python
>>> data = np.loadtxt('data/populations.txt')
>>> populations = np.ma masked_array(data[:,1:])
>>> year = data[:,0]

>>> badYears = (((year >= 1903) & (year <= 1910))
...             | ((year >= 1917) & (year <= 1918)))
>>> # '&' means 'and' and '/' means 'or'
>>> populations[badYears, 0] = np.ma masked
>>> populations[badYears, 1] = np.ma masked

>>> populations.mean(axis=0)
masked_array(data = [40472.72727272727 18627.272727272728 42400.0],
              mask = [False False False],
              fill_value = 1e+20)

>>> populations.std(axis=0)
masked_array(data = [21087.656489006717 15625.799814240254 3322.5062255844787],
              mask = [False False False],
              fill_value = 1e+20)
```

Note that Matplotlib knows about masked arrays:

```python
>>> plt.plot(year, populations, 'o-')
[< matplotlib.lines Line2D object at ...>, ...]
```

<figure><img src="image_127.png" /></figure>
### 8.4.3 recarray: purely convenience

```python
>>> arr = np.array([['a', 1], ('b', 2)], dtype=[('x', 'S1'), ('y', int)])
>>> arr2 = arr.view(np.recarray)
>>> arr2.x
chararray(['a', 'b'],
           type='|S1')
>>> arr2.y
array([1, 2])
```

### 8.4.4 matrix: convenience?

* always 2-D

* `*` is the matrix product, not the elementwise one

```python
>>> np.matrix([[1, 0], [0, 1]]) * np.matrix([[1, 2], [3, 4]])
matrix([[1, 2],
        [3, 4]])
```

## 8.5 Summary

* Anatomy of the ndarray: data, type, strides.

* Universal functions: elementwise operations, how to make new ones

* Ndarray subclasses

* Various buffer interfaces for integration with other tools

* Recent additions: PEP 3118, generalized ufunc

## 8.6 Contributing to Numpy/Scipy

Get this tutorial: http://www.euroscipy.org/talk/882

### 8.6.1 Why

* "There's a bug?"

* "I don't understand what this is supposed to do?"

* "I have this fancy code. Would you like to have it?"

* "I'd like to help! What can I do?"

### 8.6.2 Reporting bugs

* Bug tracker (prefer **this**)

- http://projects.scipy.org/numpy

- http://projects.scipy.org/scipy

- Click the "Register" link to get an account

* Mailing lists (scipy.org/Mailing_Lists)

- If you're unsure
- No replies in a week or so? Just file a bug ticket.

## Good bug report

Title: numpy.random.permutations fails for non-integer arguments

I'm trying to generate random permutations, using numpy.random.permutations

When calling numpy.random.permutation with non-integer arguments
it fails with a cryptic error message::

```
>>> np.random.permutation(12)
array([ 6, 11,  4, 10,  2,  8,  1,  7,  9,  3,  0,  5])
>>> np.random.permutation(12.)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "mtrand.pyx", line 3311, in mtrand.RandomState.permutation
  File "mtrand.pyx", line 3254, in mtrand.RandomState shuf fle
TypeError: len() of unsized object
```

This also happens with long arguments, and so
`np.random.permutation(X.shape[0])` where X is an array fails on 64
bit windows (where shape is a tuple of longs).

It would be great if it could cast to integer or at least raise a
proper error for non-integer types.

I'm using Numpy 1.4.1, built from the official tarball, on Windows
64 with Visual studio 2008, on Python.org 64-bit Python.

0. What are you trying to do?

1. Small code snippet reproducing the bug (if possible)

* What actually happens

* What you'd expect

2. Platform (Windows / Linux / OSX, 32/64 bits, x86/PPC, ...)

3. Version of Numpy/Scipy

```
>>> print(np.__version__)
1...
```

Check that the following is what you expect

```
>>> print(np.__file__)
/...
```

In case you have old/broken Numpy installations lying around.

If unsure, try to remove existing Numpy installations, and reinstall...

### 8.6.3 Contributing to documentation

1. Documentation editor

* http://docs.scipy.org/doc/numpy

* Registration

- Register an account

- Subscribe to scipy-dev mailing list (subscribers-only)
- Problem with mailing lists: you get mail

* But: **you can turn mail delivery off**

* “change your subscription options”, at the bottom of
http://mail.scipy.org/mailman/listinfo/scipy-dev

- Send a mail @ scipy-dev mailing list; ask for activation:

To: scipy-dev@scipy.org

Hi,

I'd like to edit Numpy/Scipy docstrings. My account is XXXXX

Cheers,
N. N.

• Check the style guide:

- http://docs.scipy.org/doc/numpy/

- Don't be intimidated; to fix a small thing, just fix it

• Edit

2. Edit sources and send patches (as for bugs)

3. Complain on the mailing list

8.6.4 Contributing features

0. Ask on mailing list, if unsure where it should go

1. Write a patch, add an enhancement ticket on the bug tracklet

2. OR, create a Git branch implementing the feature + add enhancement ticket.

• Especially for big/invasive additions

• http://projects.scipy.org/numpy/wiki/GitMirror

• http://www.spheredev.org/wiki/Git_for_the_lazy

# Clone numpy repository
git clone --origin svn http://projects.scipy.org/git/numpy.git numpy

# Create a feature branch
git checkout -b name-of-my-feature-branch svn/trunk

<edit stuff>

git commit -a

• Create account on https://github.com (or anywhere)

• Create a new repository @ Github

• Push your work to github

git remote add github git@github:USERNAME/REPOSITORYNAME.git
git push github name-of-my-feature-branch
### 8.6.5 How to help, in general

* Bug fixes always welcome!

- What irks you most

- Browse the tracker

* Documentation work

- API docs: improvements to docstrings

* Know some Scipy module well?

- User guide

* Needs to be done eventually.

* Want to think? Come up with a Table of Contents
      http://scipy.org/Developer_Zone/UG_Toc

* Ask on communication channels:

- numpy-discussion list

- scipy-dev list
# Debugging code

**Author:** Gaël Varoquaux

This section explores tools to understand better your code base: debugging, to find and fix bugs.

It is not specific to the scientific Python community, but the strategies that we will employ are tailored to its needs.

## Prerequisites

* Numpy

* IPython

* nosetests (http://readthedocs.org/docs/nose/en/latest/)

* pyflakes (http://pypi.python.org/pypi/pyflakes)

* gdb for the C-debugging part.

## Chapter contents

* Avoiding bugs

- Coding best practices to avoid getting in trouble

- pyflakes: fast static analysis

* Debugging workflow

* Using the Python debugger

- Invoking the debugger

- Debugger commands and interaction

* Debugging segmentation faults using gdb

## 9.1 Avoiding bugs

### 9.1.1 Coding best practices to avoid getting in trouble

**Brian Kernighan**

*"Everyone knows that debugging is twice as hard as writing a program in the first place. So if you're as clever as you can be when you write it, how will you ever debug it?"*

* We all write buggy code. Accept it. Deal with it.

* Write your code with testing and debugging in mind.

* Keep It Simple, Stupid (KISS).
- What is the simplest thing that could possibly work?

* Don't Repeat Yourself (DRY).

- Every piece of knowledge must have a single, unambiguous, authoritative representation within a system.

- Constants, algorithms, etc...

* Try to limit interdependencies of your code. (Loose Coupling)

* Give your variables, functions and modules meaningful names (not mathematics names)

### 9.1.2 pyflakes: fast static analysis

They are several static analysis tools in Python; to name a few:

* `pylint`

* `pychecker`

* `pyflakes`

* `pep8`

* `flake8`

Here we focus on `pyflakes`, which is the simplest tool.

* **Fast, simple**

* Detects syntax errors, missing imports, typos on names.

Another good recommendation is the `flake8` tool which is a combination of `pyflakes` and `pep8`. Thus, in addition to the types of errors that `pyflakes` catches, `flake8` detects violations of the recommendation in `PEP8` style guide.

Integrating `pyflakes` (or `flake8`) in your editor or IDE is highly recommended, it **does yield productivity gains**.

#### Running pyflakes on the current edited file

You can bind a key to run `pyflakes` in the current buffer.

* **In kate** Menu: `settings -> configure kate`

- In plugins enable 'external tools'

- In external Tools, add `pyflakes`:

```
kdialog --title "pyflakes %filename" --msgbox "$(pyflakes %filename)"
```

* **In TextMate**

Menu: TextMate -> Preferences -> Advanced -> Shell variables, add a shell variable:

```
TM_PYCHECKER = /Library/Frameworks/Python.framework/Versions/Current/bin/pyflakes
```

Then `Ctrl-Shift-V` is bound to a `pyflakes` report

* **In vim** In your `vimrc` (binds `F5` to `pyflakes`):

```
autocmd FileType python let &mp = 'echo "*** running % ***" ; pyflakes %'
autocmd FileType tex,mp,rst,python imap <Esc>[15~ <C-O>:make!^M
autocmd FileType tex,mp,rst,python map <Esc>[15~ :make!^M
autocmd FileType tex,mp,rst,python set autowrite
```

* **In emacs** In your `emacs` (binds `F5` to `pyflakes`):
```python
(defun pyflakes-thisfile () (interactive)
  (compile (format "pyflakes %s" (buffer-file-name))
)

#define-minor-mode pyflakes-mode
"Toggle pyflakes mode.
With no argument, this command toggles the mode.
Non-null prefix argument turns on the mode.
Null prefix argument turns off the mode."
;; The initial value.
nil
;; The indicator for the mode line.
" Pyflakes"
;; The minor mode bindings.
'( ([f5] . pyflakes-thisfile) )
)

(add-hook 'python-mode-hook (lambda () (pyflakes-mode t)))
```

## A type-as-go spell-checker like integration

* In vim

- Use the `pyflakes.vim` plugin:

1. download the zip file from http://www.vim.org/scripts/script.php?script_id=2441

2. extract the files in ~/.vim/ftplugin/python

3. make sure your vimrc has `filetype` plugin indent on

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

- Alternatively: use the `syntastic` plugin. This can be configured to use `flake8` too and also handles on-the-fly checking for many other languages.

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

* In `emacs` Use the `flymake` mode with `pyflakes`, documented on http://www.plope.com/Members/chrism/flymake-mode: add the following to your `.emacs` file:

```elisp
(when (load "flymake" t)
  (defun flymake-pyflakes-init ()
    (let* ((temp-file (flymake-init-create-temp-buffer-copy
                        'flymake-create-temp-inplace)))
      (local-file (file-relative-name
                    temp-file
                    ))
```
```
(file-name-directory buffer-file-name)))
(list "pyflakes" (list local-file)))

(add-to-list 'flymake-allowed-file-name-masks
'("\\.py\\"" flymake-pyflakes-init)))

(add-hook 'find-file-hook 'flymake-find-file-hook)
```

## 9.2 Debugging workflow

If you do have a non trivial bug, this is when debugging strategies kick in. There is no silver bullet. Yet, strategies help:

**For debugging a given problem, the favorable situation is when the problem is isolated in a small number of lines of code, outside framework or application code, with short modify-run-fail cycles**

1. Make it fail reliably. Find a test case that makes the code fail every time.

2. Divide and Conquer. Once you have a failing test case, isolate the failing code.

* Which module.

* Which function.

* Which line of code.

=> isolate a small reproducible failure: a test case

3. Change one thing at a time and re-run the failing test case.

4. Use the debugger to understand what is going wrong.

5. Take notes and be patient. It may take a while.

Once you have gone through this process: isolated a tight piece of code reproducing the bug and fix the bug using this piece of code, add the corresponding code to your test suite.

## 9.3 Using the Python debugger

The python debugger, `pdb`: https://docs.python.org/library/pdb.html, allows you to inspect your code interactively.

Specifically it allows you to:

* View the source code.

* Walk up and down the call stack.

* Inspect values of variables.

* Modify values of variables.

* Set breakpoints.

### print

Yes, `print` statements do work as a debugging tool. However to inspect runtime, it is often more efficient to use the debugger.
### 9.3.1 Invoking the debugger

Ways to launch the debugger:

1. Postmortem, launch debugger after module errors.

2. Launch the module with the debugger.

3. Call the debugger inside the module

#### Postmortem

**Situation:** You're working in IPython and you get a traceback.

Here we debug the file `index_error.py`. When running it, an `IndexError` is raised. Type `%debug` and drop into the debugger.

```
In [1]: %run index_error.py
---------------------------------------------------------------------------
IndexError                                 Traceback (most recent call last)
  /home/varoquau/dev/scipy-lecture-notes/advanced/debugging/index_error.py in <module>(
    6
  7 if _name_ == '__main__':
---> 8      index_error()
        9

  /home/varoquau/dev/scipy-lecture-notes/advanced/debugging/index_error.py in index_error(
    3 def index_error():
      4      lst = list('foobar')
---> 5      print lst[len]
        6
  7 if _name_ == '__main__':
IndexError: list index out of range

In [2]: %debug
  /home/varoquau/dev/scipy-lecture-notes/advanced/debugging/index_error.py(5)index_error(
    4      lst = list('foobar')
---> 5      print lst[len lst])
        6

ipdb> list
  1 """Small snippet to raise an IndexError."""
  2
  3 def index_error():
      4      lst = list('foobar')
---> 5      print lst[len lst])
        6
  7 if _name_ == '__main__':
      8      index_error()
        9

ipdb> len lst)
  6
ipdb> print lst[len lst)-1]
  r
ipdb> quit
In [3]:
```
## Post-mortem debugging without IPython

In some situations you cannot use IPython, for instance to debug a script that wants to be called from the command line. In this case, you can call the script with `python -m pdb script.py`:

```python
$ python -m pdb index_error.py
> /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/index_error.py(1)<module>()
-> """Small snippet to raise an IndexError."""
(Pdb) continue
Traceback (most recent call last):
  File "/usr/lib/python2.6/pdb.py", line 1296, in main
    pdb._runscript(mainpyfile)
  File "/usr/lib/python2.6/pdb.py", line 1215, in _runscript
    self.run statement)
  File "/usr/lib/python2.6/bdb.py", line 372, in run
    exec cmd in globals, locals
  File "<string>", line 1, in <module>
  File "index_error.py", line 8, in <module>
    index_error()
  File "index_error.py", line 5, in index_error
    print lst[len(lst)]
IndexError: list index out of range
Uncaught exception. Entering post mortem debugging
Running 'cont' or 'step' will restart the program
> /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/index_error.py(5)index_error()
-> print lst[len(lst)]
(Pdb)
```

## Step-by-step execution

**Situation:** You believe a bug exists in a module but are not sure where.

For instance we are trying to debug `wiener_filtering.py`. Indeed the code runs, but the filtering does not work well.

* Run the script in IPython with the debugger using `%run -d wiener_filtering.py`:

```python
In [1]: %run -d wiener_filtering.py
*** Blank or comment
*** Blank or comment
*** Blank or comment
Breakpoint 1 at /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/wiener_filtering.py:4
NOTE: Enter 'c' at the ipdb> prompt to start your script.
> <string>(1)<module>()
```

* Set a break point at line 34 using b 34:

```python
ipdb> n
> /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/wiener_filtering.py(4)<module>()
3
1---> 4 import numpy as np
5 import scipy as sp

ipdb> b 34
Breakpoint 2 at /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/wiener_filtering.py:34
```

* Continue execution to next breakpoint with c (`cont(inue)`):

```python
ipdb> c
> /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/wiener_filtering.py(34)iterated_wiener()
33 """
```
```
2 --> 34 noisy img = noisy img
35 denoised img = local_mean(noisy img, size=size)
```

* Step into code with `n(ext)` and `s(tep)`: next jumps to the next statement in the current execution context, while step will go across execution contexts, i.e. enable exploring inside function calls:

```
ipdb> s
> /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/wiener_filtering.py(35)iterated_wiener()
2   34 noisy img = noisy img
---> 35 denoised img = local_mean(noisy img, size=size)
36   l_var = local_var(noisy img, size=size)
```

```
ipdb> n
> /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/wiener_filtering.py(36)iterated_wiener()
35   denoised img = local_mean(noisy img, size=size)
---> 36 l_var = local_var(noisy img, size=size)
37   for i in range(3):
```

* Step a few lines and explore the local variables:

```
ipdb> n
> /home/varoquau/dev/scipy-lecture-notes/advanced/optimizing/wiener_filtering.py(37)iterated_wiener()
36   l_var = local_var(noisy img, size=size)
---> 37 for i in range(3):
38     res = noisy img - denoised img
ipdb> print l_var
[[5868 5379 5316 ..., 5071 4799 5149]
[5013 363 437 ..., 346 262 4355]
[5379 410 344 ..., 392 604 3377]
...,
[ 435 362 308 ..., 275 198 1632]
[ 548 392 290 ..., 248 263 1653]
[ 466 789 736 ..., 1835 1725 1940]]
ipdb> print l_var.min()
0
```

Oh dear, nothing but integers, and 0 variation. Here is our bug, we are doing integer arithmetic.
## Raising exception on numerical errors

When we run the `wiener_filtering.py` file, the following warnings are raised:

```
In [2]: %run wiener_filtering.py
wiener_filtering.py:40: RuntimeWarning: divide by zero encountered in divide
    noise_level = (1 - noise/l_var)
```

We can turn these warnings in exception, which enables us to do post-mortem debugging on them, and find our problem more quickly:

```
In [3]: np.seterr(all='raise')
Out[3]: {'divide': 'print', 'invalid': 'print', 'over': 'print', 'under': 'ignore'}
In [4]: %run wiener_filtering.py
---------------------------------------------------------------------------
FloatingPointError                      Traceback (most recent call last)
  /home/esc/anaconda/lib/python2.7/site-packages/IPython/utilS/py3compat.pyc in execfile(fname, *where)
    176                  else:
    177                  filename = fname
 --> 178                  __builtins__.execfile(filename, *where)

/home/esc/physique-cuso-python-2013/scipy-lecture-notes/advanced/debugging/wiener_filtering.py in <module>(
    55 pl.matshow(noisy_lena[cut], cmap=pl.cm.Gray)
    56
 --> 57 denoised_lena = iterated_wiener(noisy_lena)
    58 pl.matshow(denoised_lena[cut], cmap=pl.cm.Gray)
    59

/home/esc/physique-cuso-python-2013/scipy-lecture-notes/advanced/debugging/wiener_filtering.py in <module>(
    38          res = noisy img - denoised img
    39          noise = (res**2).sum()/res.size
 --> 40          noise_level = (1 - noise/l_var )
    41          noise_level[noise_level<0] = 0
    42          denoised img += noise_level*res

FloatingPointError: divide by zero encountered in divide
```

## Other ways of starting a debugger

* **Raising an exception as a poor man break point**

If you find it tedious to note the line number to set a break point, you can simply raise an exception at the point that you want to inspect and use IPython's `%debug`. Note that in this case you cannot step or continue the execution.

* **Debugging test failures using nosetests**

You can run `nose tests --pdb` to drop in post-mortem debugging on exceptions, and `nose tests --pdb-failure` to inspect test failures using the debugger.

In addition, you can use the IPython interface for the debugger in nose by installing the nose plugin `ipdbplugin`. You can then pass `--ipdb` and `--ipdb-failure` options to `nose tests`.

* **Calling the debugger explicitly**

Insert the following line where you want to drop in the debugger:

```
import pdb; pdb.set_trace()
```

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

When running `nose tests`, the output is captured, and thus it seems that the debugger does not work. Simply run the `nose tests` with the `-s` flag.
## Graphical debuggers and alternatives

* For stepping through code and inspecting variables, you might find it more convenient to use a graphical debugger such as `winpdb`.

* Alternatively, `pudb` is a good semi-graphical debugger with a text user interface in the console.

* Also, the `pydbgr` project is probably worth looking at.

### 9.3.2 Debugger commands and interaction

<table><tr><td>l(list)</td><td>Lists the code at the current position</td></tr><tr><td>u(p)</td><td>Walk up the call stack</td></tr><tr><td>d(own)</td><td>Walk down the call stack</td></tr><tr><td>n(ext)</td><td>Execute the next line (does not go down in new functions)</td></tr><tr><td>s(tep)</td><td>Execute the next statement (goes down in new functions)</td></tr><tr><td>bt</td><td>Print the call stack</td></tr><tr><td>a</td><td>Print the local variables</td></tr><tr><td>!command</td><td>Execute the given Python command (by opposition to `pdb` commands)</td></tr></table>

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

**Debugger commands are not Python code**

You cannot name the variables the way you want. For instance, if in you cannot override the variables in the current frame with the same name: **use different names than your local variable when typing code in the debugger.**

#### Getting help when in the debugger

Type `h` or `help` to access the interactive help:

```
ipdb> help
```

Documented commands (type help <topic>):

========================================================

EOF      bt          cont       enable     jump     pdef      r          tbreak     w
a        c           continue   exit       l         pdoc      restart    u          whatis
alias    cl          d           h          list      pinfo     return     unalias   where
args     clear       debug      help       n         pp        run        unt
b        commands    disable     ignore     next      q         s          until
break    condition   down       j          p         quit      step       up

Miscellaneous help topics:

========================

exec  pdb

Undocumented commands:

========================

retval rv

## 9.4 Debugging segmentation faults using gdb

If you have a segmentation fault, you cannot debug it with `pdb`, as it crashes the Python interpreter before it can drop in the debugger. Similarly, if you have a bug in C code embedded in Python, `pdb` is useless. For this we turn to the `gnu` debugger, `gdb`, available on Linux.

Before we start with `gdb`, let us add a few Python-specific tools to it. For this we add a few macros to our `~/.gdbinit`. The optimal choice of macro depends on your Python version and your `gdb` version. I have added a simplified version in `gdbinit`, but feel free to read *DebuggingWithGdb*.
To debug with gdb the Python script `segfault.py`, we can run the script in gdb as follows

```
$ gdb python
...
(gdb) run segfault.py
Starting program: /usr/bin/python segfault.py
[Thread debugging using libthread_db enabled]

Program received signal SIGSEGV, Segmentation fault.
_strided_byte_copy (dst=0x8537478 "\360\343G", outstrides=4, src=
0x86c0690 <Address 0x86c0690 out of bounds>, instrides=32, N=3,
elsize=4)
    at numpy/core/src/multiarray/actors.c:365
365         _FAST_MOVE(Int32);
(gdb)
```

We get a segfault, and gdb captures it for post-mortem debugging in the C level stack (not the Python call stack). We can debug the C call stack using gdb's commands:

```
(gdb) up
#1  0x004af4f5 in _copy_from same shape (dest=<value optimized out>,
    src=<value optimized out>, myfunc=0x496780 <_strided_byte_copy>,
    swap=0)
at numpy/core/src/multiarray/actors.c:748
748         myfunc(dit->dataptr, dest->strides(maxaxis],
```

As you can see, right now, we are in the C code of numpy. We would like to know what is the Python code that triggers this segfault, so we go up the stack until we hit the Python execution loop:

```
(gdb) up
#8  0x080ddd23 in call_function (f=
    Frame 0x85371ec, for file /home/varoquau/usr/lib/python2.6/site-packages/numpy/core/arrayprint.py, line 156
    at ../Python/ceval.c:3750
3750      ../Python/ceval.c: No such file or directory.
          in ../Python/ceval.c

(gdb) up
#9  PyEval_EvalFrameEx (f=
    Frame 0x85371ec, for file /home/varoquau/usr/lib/python2.6/site-packages/numpy/core/arrayprint.py, line 156
    at ../Python/ceval.c:2412
2412      in ../Python/ceval.c
(gdb)
```

Once we are in the Python execution loop, we can use our special Python helper function. For instance we can find the corresponding Python code:

```
(gdb) pyframe
/home/varoquau/usr/lib/python2.6/site-packages/numpy/core/arrayprint.py (158): _leading_trailing
(gdb)
```

This is numpy code, we need to go up until we find code that we have written:

```
(gdb) up
...
(gdb) up
#34 0x080dc97a in PyEval_EvalFrameEx (f=
    Frame 0x82f064c, for file segfault.py, line 11, in print_big_array (small_array=<numpy.ndarray at remote 0x
    1630      ../Python/ceval.c: No such file or directory.
          in ../Python/ceval.c
(gdb) pyframe
segfault.py (12): print_big_array
```

The corresponding code is:
```python
def make_big_array(small_array):
    big_array = stride_tricks.as_strided(small_array,
                                     shape=(2e6, 2e6), strides=(32, 32))
    return big_array

def print_big_array(small_array):
    big_array = make_big_array(small_array)
```

Thus the segfault happens when printing `big_array[-10:]`. The reason is simply that `big_array` has been allocated with its end outside the program memory.

For a list of Python-specific commands defined in the `gdbinit`, read the source of this file.

## Wrap up exercise

The following script is well documented and hopefully legible. It seeks to answer a problem of actual interest for numerical computing, but it does not work... Can you debug it?

**Python source code:** `to_debug.py`
Optimizing code

Donald Knuth

"Premature optimization is the root of all evil"

Author: Gaël Varoquaux

This chapter deals with strategies to make Python code go faster.

Prerequisites

* line_profiler

Chapters contents

* Optimization workflow

* Profiling Python code

- Timeit

- Profiler

- Line-profiler

* Making code go faster

- Algorithmic optimization

* Example of the SVD

* Writing faster numerical code

- Additional Links

10.1 Optimization workflow

1. Make it work: write the code in a simple **legible** ways.

2. Make it work reliably: write automated test cases, make really sure that your algorithm is right and that if you break it, the tests will capture the breakage.

3. Optimize the code by profiling simple use-cases to find the bottlenecks and speeding up these bottlenecks, finding a better algorithm or implementation. Keep in mind that a trade off should be found between profiling on a realistic example and the simplicity and speed of execution of the code. For efficient work, it is best to work with profiling runs lasting around 10s.

10.2 Profiling Python code
No optimization without measuring!

* Measure: profiling, timing

* You'll have surprises: the fastest code is not always what you think

10.2.1 Timeit

In IPython, use `timeit` (https://docs.python.org/library/timeit.html) to time elementary operations:

In [1]: import numpy as np

In [2]: a = np.arange(1000)

In [3]: %timeit a ** 2
100000 loops, best of 3: 5.73 us per loop

In [4]: %timeit a ** 2.1
1000 loops, best of 3: 154 us per loop

In [5]: %timeit a * a
100000 loops, best of 3: 5.56 us per loop

Use this to guide your choice between strategies.

For long running calls, using `%time` instead of `%timeit`; it is less precise but faster

10.2.2 Profiler

Useful when you have a large program to profile, for example the following file:

# For this example to run, you also need the 'ica.py' file

import numpy as np
from scipy import linalg

from ica import fastica

def test():
    data = np.random.random((5000, 100))
    u, s, v = linalg.svd(data)
    pca = np.dot(u[, :10].T, data)
    results = fastica(pca.T, whiten=False)

if __name__ == '__main__':
    test()

This is a combination of two unsupervised learning techniques, principal component analysis (PCA) and in-
dependent component analysis (ICA). PCA is a technique for dimensionality reduction, i.e. an algorithm to
explain the observed variance in your data using less dimensions. ICA is a source seperation technique, for ex-
ample to unmix multiple signals that have been recorded through multiple sensors. Doing a PCA first and then
an ICA can be useful if you have more sensors than signals. For more information see: the FastICA example
from scikits-learn.

To run it, you also need to download the ica module. In IPython we can time the script:

In [1]: %run -t demo.py

IPython CPU timings (estimated):
```
User : 14.3929 s.
System: 0.256016 s.
```

and profile it:

```
In [2]: %run -p demo.py

916 function calls in 14.551 CPU seconds
```

Ordered by: internal time

```
ncalls   tottime  percall  cumtime  percall filename:eno (function)
      1   14.457   14.457   14.479   14.479 decomp.py:849 (svd)
      1      0.054      0.054      0.054      0.054 {method 'random_sample' of 'mtrand.RandomState' objects}
      1      0.017      0.017      0.021      0.021 function_base.py:645 (asarray_chkfinite)
     54      0.011      0.000      0.011      0.000 {numpy.core._dotblas.dot}
      2      0.005      0.002      0.005      0.002 {method 'any' of 'numpy.ndarray' objects}
      6      0.001      0.000      0.001      0.000 ica.py:195 (gprime)
      6      0.001      0.000      0.001      0.000 ica.py:192 (g)
     14      0.001      0.000      0.001      0.000 {numpy.lINALG.LAPACK lite.dsyevd}
     19      0.001      0.000      0.001      0.000 twodim_base.py:204 (diag)
      1      0.001      0.001      0.008      0.008 ica.py:69 (_ica_par)
      1      0.001      0.001   14.551   14.551 {execfile}
     107      0.000      0.000      0.001      0.000 defmatrix.py:239 (_array_finalize_)
      7      0.000      0.000      0.004      0.001 ica.py:58 (_symdecorrelation)
      7      0.000      0.000      0.002      0.000 lINALG.py:841 (eigh)
     172      0.000      0.000      0.000      0.000 {isinstance}
      1      0.000      0.000   14.551   14.551 demo.py:1 (<module>)
      29      0.000      0.000      0.000      0.000 numeric.py:180 (asarray)
     35      0.000      0.000      0.000      0.000 defmatrix.py:193 (_new_)
     35      0.000      0.000      0.001      0.000 defmatrix.py:43 (asmatrix)
     21      0.000      0.000      0.001      0.000 defmatrix.py:287 (_mul_)
     41      0.000      0.000      0.000      0.000 {numpy.core.multiparray.zeros}
     28      0.000      0.000      0.000      0.000 {method 'transpose' of 'numpy.ndarray' objects}
      1      0.000      0.000      0.008      0.008 ica.py:97 (fastica)
...
```

Clearly the svd (in decomp.py) is what takes most of our time, a.k.a. the bottleneck. We have to find a way to make this step go faster, or to avoid this step (algorithmic optimization). Spending time on the rest of the code is useless.

### Profiling outside of IPython, running "cProfile"

Similar profiling can be done outside of IPython, simply calling the built-in Python profilers `cProfile` and `profile`.

```
$ python -m cProfile -o demo.prof demo.py
```

Using the `-o` switch will output the profiler results to the file `demo.prof` to view with an external tool. This can be useful if you wish to process the profiler output with a visualization tool.

### 10.2.3 Line-profiler

The profiler tells us which function takes most of the time, but not where it is called.

For this, we use the `line_profiler`: in the source file, we decorate a few functions that we want to inspect with `@profile` (no need to import it)

```
@profile
def test():
    data = np.random.random((5000, 100))
```
```python
u, s, v = lINALG.svd(data)
pca = np.dot(u[:, :10], data)
results = FASTICA(pca.T, whiten=False)
```

Then we run the script using the `kernprof.py` program, with switches `-l`, `--line-by-line` and `-v`, `--view` to use the line-by-line profiler and view the results in addition to saving them:

```bash
$ kernprof.py -l -v demo.py

Wrote profile results to demo.py.lprof
Timer unit: 1e-06 s

File: demo.py
Function: test at line 5
Total time: 14.2793 s

Line #   Hits   Time  Per Hit  % Time  Line Contents
========== ========= = ======== ========= = ======== =========
5                   @profile
6                   def test():
7           1      19015  19015.0   0.1      data = np.random.random((5000, 100))
8           1      14242163  14242163.0  99.7      u, s, v = lINALG.svd(data)
9           1      10282   10282.0   0.1      PCA = np.dot(u[:10, :], data)
10          1      7799   7799.0   0.1      RESULTS = FASTICA(pca.T, whiten=False)
```

The SVD is taking all the time. We need to optimise this line.

## 10.3 Making code go faster

Once we have identified the bottlenecks, we need to make the corresponding code go faster.

### 10.3.1 Algorithmic optimization

The first thing to look for is algorithmic optimization: are there ways to compute less, or better?

For a high-level view of the problem, a good understanding of the maths behind the algorithm helps. However, it is not uncommon to find simple changes, like **moving computation or memory allocation outside a for loop**, that bring in big gains.

#### Example of the SVD

In both examples above, the SVD - *Singular Value Decomposition* - is what takes most of the time. Indeed, the computational cost of this algorithm is roughly $n^3$ in the size of the input matrix.

However, in both of these examples, we are not using all the output of the SVD, but only the first few rows of its first return argument. If we use the svd implementation of `scipy`, we can ask for an incomplete version of the SVD. Note that implementations of linear algebra in `scipy` are richer than those in `numpy` and should be preferred.

```python
In [3]: %timeit np.lINALG.svd(data)
1 loops, best of 3: 14.5 s per loop

In [4]: from scipy import lINALG

In [5]: %timeit lINALG.svd(data)
1 loops, best of 3: 14.2 s per loop

In [6]: %timeit lINALG.svd(data, fullFalse)
1 loops, best of 3: 295 ms per loop
```
In [7]: %timeit np lilalg.svd(data, full佔mrices=False)
1 loops, best of 3: 293 ms per loop

We can then use this insight to optimize the previous code:

def test():
    data = np.random.random((5000, 100))
    u, s, v = lilalg.svd(data, full佔mrices=False)
    pca = np.dot(u[:, :10].T, data)
    results = fastica(pca.T, whiten=False)

In [1]: import demo

In [2]: %timeit demo.
demo.fastica demo.np demo.prof.pdf demo.py demo.pyc
demo.lilalg demo.prof demo.prof.png demo.py.lprof demo.test

In [2]: %timeit demo.test()
ica.py:65: RuntimeWarning: invalid value encountered in sqrt
W = (u * np(1.0/np.sq
1 loops, best of 3: 17.5 s per loop

In [3]: import demo_opt

In [4]: %timeit demo_opt.test()
1 loops, best of 3: 208 ms per loop

Real incomplete SVDs, e.g. computing only the first 10 eigenvectors, can be computed with arpack, available
in scipy sparse lalg eigsh.

Computational linear algebra

For certain algorithms, many of the bottlenecks will be linear algebra computations. In this case, using
the right function to solve the right problem is key. For instance, an eigenvalue problem with a symmetric
matrix is easier to solve than with a general matrix. Also, most often, you can avoid inverting a matrix and
use a less costly (and more numerically stable) operation.

Know your computational linear algebra. When in doubt, explore `scipy.lINALg`, and use `%timeit` to try out different alternatives on your data.

10.4 Writing faster numerical code

A complete discussion on advanced use of numpy is found in chapter *Advanced Numpy*, or in the article *The NumPy array: a structure for efficient numerical computation* by van der Walt et al. Here we discuss only some commonly encountered tricks to make code faster.

* **Vectorizing for loops**

Find tricks to avoid for loops using numpy arrays. For this, masks and indices arrays can be useful.

* Broadcasting

Use *broadcasting* to do operations on arrays as small as possible before combining them.

* In place operations

In [1]: a = np Article of 1e7)

In [2]: %timeit global a ; a = 0*a
10 loops, best of 3: 111 ms per loop
```
In [3]: %timeit global a; a *= 0
10 loops, best of 3: 48.4 ms per loop
```

**note:** we need `global a` in the `timeit` so that it work, as it is assigning to `a`, and thus considers it as a local variable.

* **Be easy on the memory: use views, and not copies**

Copying big arrays is as costly as making simple numerical operations on them:

```
In [1]: a = np(zeros(1e7))

In [2]: %timeit a.copy()
10 loops, best of 3: 124 ms per loop

In [3]: %timeit a + 1
10 loops, best of 3: 112 ms per loop
```

* **Beware of cache effects**

Memory access is cheaper when it is grouped: accessing a big array in a continuous way is much faster than random access. This implies amongst other things that **smaller strides are faster** (see *CPU cache effects*):

```
In [1]: c = np(zeros((1e4, 1e4), order='C'))

In [2]: %timeit c.sum(axis=0)
1 loops, best of 3: 3.89 s per loop

In [3]: %timeit c.sum(axis=1)
1 loops, best of 3: 188 ms per loop

In [4]: c.strides
Out[4]: (80000, 8)
```

This is the reason why Fortran ordering or C ordering may make a big difference on operations:

```
In [5]: a = np.randomrand(20, 2**18)

In [6]: b = np.randomrand(20, 2**18)

In [7]: %timeit np.dot(b, a.T)
1 loops, best of 3: 194 ms per loop

In [8]: c = np.ascontiguousarray(a.T)

In [9]: %timeit np.dot(b, c)
10 loops, best of 3: 84.2 ms per loop
```

Note that copying the data to work around this effect may not be worth it:

```
In [10]: %timeit c = np.ascontiguousarray(a.T)
10 loops, best of 3: 106 ms per loop
```

Using `numexpr` can be useful to automatically optimize code for such effects.

* **Use compiled code**

The last resort, once you are sure that all the high-level optimizations have been explored, is to transfer the hot spots, i.e. the few lines or functions in which most of the time is spent, to compiled code. For compiled code, the preferred option is to use *Cython*: it is easy to transform exiting Python code in compiled code, and with a good use of the * numpy support* yields efficient code on `numpy` arrays, for instance by unrolling loops.

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

For all the above: profile and time your choices. Don't base your optimization on theoretical considerations.
### 10.4.1 Additional Links

* If you need to profile memory usage, you could try the `memory_profiler`

* If you need to profile down into C extensions, you could try using `gperftools` from Python with `yep`.

* If you would like to track performance of your code across time, i.e. as you make new commits to your repository, you could try: `vbench`

* If you need some interactive visualization why not try `RunSnakeRun`
# Sparse Matrices in SciPy

Author: Robert Cimrman

## 11.1 Introduction

(dense) matrix is:

* mathematical object

* data structure for storing a 2D array of values

important features:

* memory allocated once for all items

- usually a contiguous chunk, think NumPy ndarray

* fast access to individual items (*)

### 11.1.1 Why Sparse Matrices?

* the memory, that grows like $n**2$

* small example (double precision matrix):

```python
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(0, 1e6, 10)
>>> plt plot(x, 8.0 * (x**2) / 1e6, lw=5)
< matplotlib lines Line2D object at ... >
>>> plt xlabel('size n')
< matplotlib text Text object at ... >
>>> plt ylabel('memory [MB]')
< matplotlib text Text object at ... >
```

### 11.1.2 Sparse Matrices vs. Sparse Matrix Storage Schemes

* sparse matrix is a matrix, which is almost empty

* storing all the zeros is wasteful -> store only nonzero items

* think compression

* pros: huge memory savings
* cons: depends on actual storage scheme, (*) usually does not hold

### 11.1.3 Typical Applications

* **solution of partial differential equations (PDEs)**

- the *finite element method*

- mechanical engineering, electrotechnics, physics, ...

* **graph theory**

- nonzero at (i, j) means that node i is connected to node j

* ...

### 11.1.4 Prerequisites

recent versions of

* numpy

* scipy

* matplotlib (optional)

* ipython (the enhancements come handy)

### 11.1.5 Sparsity Structure Visualization

* spy() from matplotlib

* example plots:

<figure><img src="image_133.png" /></figure>
<figure><img src="image_134.png" /></figure>

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

## 11.2 Storage Schemes

* seven sparse matrix types in `scipy sparse`:

1. `csc_matrix`: Compressed Sparse Column format

2. `csr_matrix`: Compressed Sparse Row format

3. `bsr_matrix`: Block Sparse Row format

4. `lil_matrix`: List of Lists format

5. `dok_matrix`: Dictionary of Keys format

6. `coo_matrix`: COOrdinate format (aka IJV, triplet format)

7. `dia_matrix`: DIAgonal format

* each suitable for some tasks

* many employ sparsetools C++ module by Nathan Bell

* assume the following is imported:

```python
>>> import numpy as np
>>> import scipy sparse as sps
>>> import matplotlib.pyplot as plt
```

* warning for NumPy users:

- the multiplication with '*' is the *matrix multiplication* (dot product)

- not part of NumPy!
* passing a sparse matrix object to NumPy functions expecting `ndarray/matrix` does not work

### 11.2.1 Common Methods

* all `scipy sparse` classes are subclasses of `spmatrix`

- default implementation of arithmetic operations

* always converts to CSR

* subclasses override for efficiency

- shape, data type set/get

- nonzero indices

- format conversion, interaction with NumPy (`toarray()`, `todense()`)

- ...

* attributes:

- `mtx.A` - same as `mtx.toarray()`

- `mtx.T` - transpose (same as `mtx transpose()`) - fixed length -> waste space a bit when far from main diagonal

- `mtx.H` - Hermitian (conjugate) transpose

- `mtx.real` - real part of complex matrix

- `mtx.imag` - imaginary part of complex matrix

- `mtx.size` - the number of nonzeros (same as `self.getnnz()`) - subclass of `_data_matrix` (sparse matrix classes with data attribute)

- `mtx.shape` - the number of rows and columns (tuple)

* data usually stored in NumPy arrays

### 11.2.2 Sparse Matrix Classes

#### Diagonal Format (DIA)

* very simple scheme

* diagonals in dense NumPy array of shape (`n_diag`, `length`)

- fixed length -> waste space a bit when far from main diagonal

- subclass of `_data_matrix` (sparse matrix classes with data attribute)

* offset for each diagonal

- 0 is the main diagonal

- negative offset = below

- positive offset = above

* fast matrix * vector (`sparsetools`)

* fast and easy item-wise operations

- manipulate data array directly (fast NumPy machinery)

* constructor accepts:

- dense matrix (array)

- sparse matrix
- shape tuple (create empty matrix)

- (data, offsets) tuple

* no slicing, no individual item access

* use:

- rather specialized

- solving PDEs by finite differences

- with an iterative solver

Examples

* create some DIA matrices:

```python
>>> data = np.array([[1, 2, 3, 4]]).repeat(3, axis=0)
>>> data
array([[1, 2, 3, 4],
        [1, 2, 3, 4],
        [1, 2, 3, 4]])
>>> offsets = np.array([0, -1, 2])
>>> mtx = sparse.dia_matrix((data, offsets), shape=(4, 4))
>>> mtx
classical matrix of type '<... ' numpy.int64'>
with 9 stored elements (3 diagonals) in DIAagonal format>
>>> mtx.todense()
matrix([[1, 0, 3, 0],
        [1, 2, 0, 4],
        [0, 2, 3, 0],
        [0, 0, 3, 4]])
>>> data = np.arange(12).reshape((3, 4)) + 1
>>> data
array([[ 1,  2,  3,  4],
        [ 5,  6,  7,  8],
        [ 9, 10, 11, 12]])
>>> mtx = sparse.dia_matrix((data, offsets), shape=(4, 4))
>>> mtx.data
array([[ 1,  2,  3,  4],
        [ 5,  6,  7,  8],
        [ 9, 10, 11, 12]])...
>>> mtx Offsets
array([ 0, -1,  2], dtype=int32)
>>> print(mtx)
(0, 0)      1
(1, 1)      2
(2, 2)      3
(3, 3)      4
(1, 0)      5
(2, 1)      6
(3, 2)      7
(0, 2)     11
(1, 3)     12
>>> mtx.todense()
matrix([[ 1,  0, 11,  0],
        [ 5,  2,  0, 12],
        [ 0,  6,  3,  0],
        [ 0,  0,  7,  4]])
```

* explanation with a scheme:
```
offset: row

    2:  9
    1: --10------
    0:  1 . 11 .
-1:  5  2 . 12
-2: .  6  3 .
-3: .  .  7  4
-------8
```

* matrix-vector multiplication

```python
>>> vec = np.ones((4, ))
>>> vec
array([ 1.,  1.,  1.,  1.])
>>> mtx * vec
array([ 12.,  19.,  9., 11.])
>>> mtx.toarray() * vec
array([[ 1.,  0., 11.,  0.],
      [ 5.,  2.,  0., 12.],
      [ 0.,  6.,  3.,  0.],
      [ 0.,  0.,  7.,  4.]])
```

## List of Lists Format (LIL)

* row-based linked list

- each row is a Python list (sorted) of column indices of non-zero elements

- rows stored in a NumPy array (dtype(npobject))

- non-zero values stored analogously

* efficient for constructing sparse matrices incrementally

* constructor accepts:

- dense matrix (array)

- sparse matrix

- shape tuple (create empty matrix)

* flexible slicing, changing sparsity structure is efficient

* slow arithmetics, slow column slicing due to being row-based

* use:

- when sparsity pattern is not known apriori or changes

- example: reading a sparse matrix from a text file

## Examples

* create an empty LIL matrix:

```python
>>> mtx = sparse.lil_matrix((4, 5))
```

* prepare random data:

```python
>>> from numpy.random import rand
>>> data = np_round(rand(2, 3))
>>> data
array([[ 1.,  1.,  1.],
      [ 1.,  0.,  1.]])
```
* assign the data using fancy indexing:

```python
>>> mtx[:2, [1, 2, 3]] = data
>>> mtx
<4x5 sparse matrix of type '<... 'numpy.float64'>'
  with 5 stored elements in LInked List format>
>>> print(mtx)
(0, 1)  1.0
(0, 2)  1.0
(0, 3)  1.0
(1, 1)  1.0
(1, 3)  1.0
>>> mtx.todense()
matrix([[ 0.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])
>>> mtx.toarray()
array([[ 0.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])
```

* more slicing and indexing:

```python
>>> mtx = sparse.lil_matrix([[0, 1, 2, 0], [3, 0, 1, 0], [1, 0, 0, 1]])
>>> mtx.todense()
matrix([[0, 1, 2, 0],
        [3, 0, 1, 0],
        [1, 0, 0, 1]]...)
>>> print(mtx)
(0, 1)  1
(0, 2)  2
(1, 0)  3
(1, 2)  1
(2, 0)  1
(2, 3)  1
>>> mtx[:2, :]
<2x4 sparse matrix of type '<... 'numpy.int64'>'
  with 4 stored elements in LInked List format>
>>> mtx[:2, :].todense()
matrix([[0, 1, 2, 0],
        [3, 0, 1, 0]]...)
>>> mtx[1:2, [0,2]].todense()
matrix([[3, 1]]...)
>>> mtx.todense()
matrix([[0, 1, 2, 0],
        [3, 0, 1, 0],
        [1, 0, 0, 1]]...)
```

## Dictionary of Keys Format (DOK)

* subclass of Python `dict`

- keys are `(row, column)` index tuples (no duplicate entries allowed)

- values are corresponding non-zero values

* efficient for constructing sparse matrices incrementally

* constructor accepts:

- dense matrix (array)

- sparse matrix
- shape tuple (create empty matrix)

* efficient O(1) access to individual elements

* flexible slicing, changing sparsity structure is efficient

* can be efficiently converted to a `coo_matrix` once constructed

* slow arithmetics (for loops with `dict iteritems()`)

### use:

- when sparsity pattern is not known apriori or changes

## Examples

* create a DOK matrix element by element:

```python
>>> mtx = sparse.dok_matrix((5, 5), type= np.float64)
>>> mtx
<5x5 sparse matrix of type '<... 'numpy.float64'>'
  with 0 stored elements in Dictionary Of Keys format>
>>> for ir in range(5):
...     for ic in range(5):
...         mtx[ir, ic] = 1.0 * (ir != ic)
>>> mtx
<5x5 sparse matrix of type '<... 'numpy.float64'>'
  with 20 stored elements in Dictionary Of Keys format>
>>> mtx.todense()
matrix([[ 0.,  1.,  1.,  1.,  1.],
        [ 1.,  0.,  1.,  1.,  1.],
        [ 1.,  1.,  0.,  1.,  1.],
        [ 1.,  1.,  1.,  0.,  1.],
        [ 1.,  1.,  1.,  1.,  0.]])
```

* slicing and indexing:

```python
>>> mtx[1, 1]
0.0
>>> mtx[1, 1:3]
<1x2 sparse matrix of type '<... 'numpy.float64'>'
  with 1 stored elements in Dictionary Of Keys format>
>>> mtx[1, 1:3].todense()
matrix([[ 0.,  1.]])
>>> mtx[[2,1], 1:3].todense()
matrix([[ 1.,  0.],
        [ 0.,  1.]])
```

## Coordinate Format (COO)

* also known as the 'ijv' or 'triplet' format

- three NumPy arrays: `row`, `col`, `data`

- `data[i]` is value at `(row[i], col[i])` position

- permits duplicate entries

- subclass of `_data_matrix` (sparse matrix classes with `data` attribute)

* fast format for constructing sparse matrices

### constructor accepts:

- dense matrix (array)
- sparse matrix

- shape tuple (create empty matrix)

- `(data, ij)` tuple

* very fast conversion to and from CSR/CSC formats

* fast matrix \* vector (sparsetools)

* **fast and easy item-wise operations**

- manipulate data array directly (fast NumPy machinery)

* no slicing, no arithmetics (directly)

* **use:**

- facilitates fast conversion among sparse formats

- when converting to other format (usually CSR or CSC), duplicate entries are summed together

* facilitates efficient construction of finite element matrices

## Examples

* create empty COO matrix:

```python
>>> mtx = sparse.coo_matrix((3, 4), dtype= np.int8)
>>> mtx.todense()
matrix([[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]]), dtype=int8)
```

* create using `(data, ij)` tuple:

```python
>>> row = np.array([0, 3, 1, 0])
>>> col = np.array([0, 3, 1, 2])
>>> data = np.array([4, 5, 7, 9])
>>> mtx = sparse.coo_matrix((data, (row, col)), shape=(4, 4))
>>> mtx
<4x4 sparse matrix of type '<... 'numpy.int64'>'
  with 4 stored elements in COOordinate format>
>>> mtx.todense()
matrix([[4, 0, 9, 0],
        [0, 7, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 5]])
```

* duplicates entries are summed together:

```python
>>> row = np.array([0, 0, 1, 3, 1, 0, 0])
>>> col = np.array([0, 2, 1, 3, 1, 0, 0])
>>> data = np.array([1, 1, 1, 1, 1, 1, 1])
>>> mtx = sparse.coo_matrix((data, (row, col)), shape=(4, 4))
>>> mtx.todense()
matrix([[3, 0, 1, 0],
        [0, 2, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 1]])
```

* no slicing...:

```python
>>> mtx[2, 3]
Traceback (most recent call last):
  ...
TypeError: 'coo_matrix' object ...
```
```
Traceback (most recent call last):
  ...
TypeError: 'coo_matrix' object ...
```

## Compressed Sparse Row Format (CSR)

* row oriented

- three NumPy arrays: `indices`, `indptr`, `data`

* `indices` is array of column indices

* `data` is array of corresponding nonzero values

* `indptr` points to row starts in `indices` and `data`

* length is `n_row + 1`, last item = number of values = length of both `indices` and `data`

* nonzero values of the i-th row are `data[indptr[i]:indptr[i+1]]` with column indices `indices[indptr[i]:indptr[i+1]]`

* item `(i, j)` can be accessed as `data[indptr[i]+k]`, where `k` is position of `j` in `indices[indptr[i]:indptr[i+1]]`

- subclass of `_cs_matrix` (common CSR/CSC functionality)

* subclass of `_data_matrix` (sparse matrix classes with `data` attribute)

* fast matrix vector products and other arithmetics (sparsetools)

* constructor accepts:

- dense matrix (array)

- sparse matrix

- shape tuple (create empty matrix)

- `(data, ij)` tuple

- `(data, indices, indptr)` tuple

* efficient row slicing, row-oriented operations

* slow column slicing, expensive changes to the sparsity structure

* use:

- actual computations (most linear solvers support this format)

## Examples

* create empty CSR matrix:

```python
>>> mtx = sparse.csr_matrix((3, 4), dtype= np.int8)
>>> mtx.todense()
matrix([[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]], ndtype=int8)
```

* create using `(data, ij)` tuple:

```python
>>> row = np.array([0, 0, 1, 2, 2, 2])
>>> col = np.array([0, 2, 2, 0, 1, 2])
>>> data = np.array([1, 2, 3, 4, 5, 6])
>>> mtx = sparse.csr_matrix((data, (row, col)), shape=(3, 3))
>>> mtx
```
```python
<3x3 sparse matrix of type '<... 'numpy.int64'>'
    with 6 stored elements in Compressed Sparse Row format>
>>> mtx.todense()
matrix([[1, 0, 2],
        [0, 0, 3],
        [4, 5, 6]])...
>>> mtx.data
array([1, 2, 3, 4, 5, 6])...
>>> mtx.indexes
array([0, 2, 2, 0, 1, 2], dtype=int32)
>>> mtx.indptr
array([0, 2, 3, 6], type=int32)
```

* create using (`data`, `indices`, `indptr`) tuple:

```python
>>> data = np.array([1, 2, 3, 4, 5, 6])
>>> indices = np.array([0, 2, 2, 0, 1, 2])
>>> indptr = np.array([0, 2, 3, 6])
>>> mtx = sparse.csr_matrix((data, indices, indptr), shape=(3, 3))
>>> mtx.todense()
matrix([[1, 0, 2],
        [0, 0, 3],
        [4, 5, 6]])
```

## Compressed Sparse Column Format (CSC)

* column oriented

- three NumPy arrays: `indices`, `indptr`, `data`

* `indices` is array of row indices

* `data` is array of corresponding nonzero values

* `indptr` points to column starts in `indices` and `data`

* length is `n_col` + 1, last item = number of values = length of both `indices` and `data`

* nonzero values of the i-th column are `data[indptr[i]:indptr[i+1]]` with row indices `indices[indptr[i]:indptr[i+1]]`

* item `(i, j)` can be accessed as `data[indptr[j]+k]`, where `k` is position of `i` in `indices[indptr[j]:indptr[j+1]]`

- subclass of `_cs_matrix` (common CSR/CSC functionality)

* subclass of `_data_matrix` (sparse matrix classes with `data` attribute)

* fast matrix vector products and other arithmetics (sparsetools)

* constructor accepts:

- dense matrix (array)

- sparse matrix

- shape tuple (create empty matrix)

- `(data, ij)` tuple

- `(data, indices, indptr)` tuple

* efficient column slicing, column-oriented operations

* slow row slicing, expensive changes to the sparsity structure

* use:

- actual computations (most linear solvers support this format)
## Examples

* create empty CSC matrix:

```python
>>> mtx = sparse.csc_matrix((3, 4), dtype= np.int8)
>>> mtx.todense()
matrix([[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]]), dtype= int8)
```

* create using (data, ij) tuple:

```python
>>> row = np.array([0, 0, 1, 2, 2, 2])
>>> col = np.array([0, 2, 2, 0, 1, 2])
>>> data = np.array([1, 2, 3, 4, 5, 6])
>>> mtx = sparse.csc_matrix((data, (row, col)), shape=(3, 3))
>>> mtx
<3x3 sparse matrix of type '<... 'numpy.int64'>'
  with 6 stored elements in Compressed Sparse Column format>
>>> mtx.todense()
matrix([[1, 0, 2],
        [0, 0, 3],
        [4, 5, 6]])...
>>> mtx.data
array([1, 4, 5, 2, 3, 6])...
>>> mtx.indexes
array([0, 2, 2, 0, 1, 2], dtype=int32)
>>> mtx.indptr
array([0, 2, 3, 6], dtype=int32)
```

* create using (data, indices, indptr) tuple:

```python
>>> data = np.array([1, 4, 5, 2, 3, 6])
>>> indices = np.array([0, 2, 2, 0, 1, 2])
>>> indptr = np.array([0, 2, 3, 6])
>>> mtx = sparse.csc_matrix((data, indices, indptr), shape=(3, 3))
>>> mtx.todense()
matrix([[1, 0, 2],
        [0, 0, 3],
        [4, 5, 6]])
```

## Block Compressed Row Format (BSR)

* basically a CSR with dense sub-matrices of fixed shape instead of scalar items

- block size (R, C) must evenly divide the shape of the matrix (M, N)

- **three NumPy arrays: `indices`, `indptr`, `data`**

* `indices` is array of column indices for each block

* `data` is array of corresponding nonzero values of shape (nnz, R, C)

* ...

- **subclass of `_cs_matrix` (common CSR/CSC functionality)**

* subclass of `_data_matrix` (sparse matrix classes with `data` attribute)

* fast matrix vector products and other arithmetics (sparsetools)

* **constructor accepts:**

- dense matrix (array)

- sparse matrix
- shape tuple (create empty matrix)

- (data, ij) tuple

- (data, indices, indptr) tuple

* many arithmetic operations considerably more efficient than CSR for sparse matrices with dense sub-matrices

* use:

- like CSR

- vector-valued finite element discretizations

Examples

* create empty BSR matrix with (1, 1) block size (like CSR...):

```python
>>> mtx = sparse.bsr_matrix((3, 4), dtype= np.int8)
>>> mtx
<3x4 sparse matrix of type '<... 'numpy.int8'>
  with 0 stored elements (blocksize = 1x1) in Block Sparse Row format>
>>> mtx.todense()
matrix([[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]], dtype=int8)
```

* create empty BSR matrix with (3, 2) block size:

```python
>>> mtx = sparse.bsr_matrix((3, 4), blocksize=(3, 2), dtype= np.int8)
>>> mtx
<3x4 sparse matrix of type '<... 'numpy.int8'>
  with 0 stored elements (blocksize = 3x2) in Block Sparse Row format>
>>> mtx.todense()
matrix([[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]], dtype=int8)
```

- a bug?

* create using (data, ij) tuple with (1, 1) block size (like CSR...):

```python
>>> row = np.array([0, 0, 1, 2, 2, 2])
>>> col = np.array([0, 2, 2, 0, 1, 2])
>>> data = np.array([1, 2, 3, 4, 5, 6])
>>> mtx = sparse.bsr_matrix((data, (row, col)), shape=(3, 3))
>>> mtx
<3x3 sparse matrix of type '<... 'numpy.int64'>
  with 6 stored elements (blocksize = 1x1) in Block Sparse Row format>
>>> mtx.todense()
matrix([[1, 0, 2],
        [0, 0, 3],
        [4, 5, 6]])...
>>> mtx.data
array([[1]]
      [[2]]
      [[3]]
      [[4]]
      [[5]]
```
```
[[6]]...)
>>> mtx.indexes
array([0, 2, 2, 0, 1, 2], dtype=int32)
>>> mtx.indptr
array([0, 2, 3, 6], type=int32)
```

* create using (`data`, `indices`, `indptr`) tuple with `(2, 2)` block size:

```
>>> indptr = np.array([0, 2, 3, 6])
>>> indices = np.array([0, 2, 2, 0, 1, 2])
>>> data = np.array([1, 2, 3, 4, 5, 6]).repeat(4).reshape(6, 2, 2)
>>> mtx = sparse.bsr_matrix((data, indices, indptr), shape=(6, 6))
>>> mtx.todense()
matrix([[1, 1, 0, 0, 2, 2],
        [1, 1, 0, 0, 2, 2],
        [0, 0, 0, 0, 3, 3],
        [0, 0, 0, 0, 3, 3],
        [4, 4, 5, 5, 6, 6],
        [4, 4, 5, 5, 6, 6]])
>>> data
array([[1, 1],
        [1, 1]]
        [[2, 2]]
        [[2, 2]]
        [[3, 3]]
        [[3, 3]]
        [[4, 4]]
        [[4, 4]]
        [[5, 5]]
        [[5, 5]]
        [[6, 6]]
        [[6, 6]])
```

### 11.2.3 Summary

<table><caption>Table 11.1: Summary of storage schemes.</caption><thead><tr><th>format</th><th>matrix * vector</th><th>get item</th><th>fancy get</th><th>set item</th><th>fancy set</th><th>solvers</th><th>note</th></tr></thead><tbody><tr><td>DIA</td><td>sparse-tools</td><td>.</td><td>.</td><td>.</td><td>.</td><td>iterative</td><td>has data array, specialized</td></tr><tr><td>LIL</td><td>via CSR</td><td>yes</td><td>yes</td><td>yes</td><td>yes</td><td>iterative</td><td>arithmetic via CSR, incremental construction</td></tr><tr><td>DOK</td><td>python</td><td>yes</td><td>one axis only</td><td>yes</td><td>yes</td><td>iterative</td><td>O(1) item access, incremental construction</td></tr><tr><td>COO</td><td>sparse-tools</td><td>.</td><td>.</td><td>.</td><td>.</td><td>iterative</td><td>has data array, facilitates fast conversion</td></tr><tr><td>CSR</td><td>sparse-tools</td><td>yes</td><td>yes</td><td>slow</td><td>.</td><td>any</td><td>has data array, fast row-wise ops</td></tr><tr><td>CSC</td><td>sparse-tools</td><td>yes</td><td>yes</td><td>slow</td><td>.</td><td>any</td><td>has data array, fast column-wise ops</td></tr><tr><td>BSR</td><td>sparse-tools</td><td>.</td><td>.</td><td>.</td><td>.</td><td>specialized</td><td>has data array, specialized</td></tr></tbody></table>
## 11.3 Linear System Solvers

* sparse matrix/eigenvalue problem solvers live in `scipy sparse linalg`

* the submodules:

- `dsolve`: direct factorization methods for solving linear systems

- `isolve`: iterative methods for solving linear systems

- `eigen`: sparse eigenvalue problem solvers

* all solvers are accessible from:

```python
>>> import scipy sparse linalg as spla
>>> spla.__all__
['LinearOperator', 'Tester', 'arpack', 'aslinearoperator', 'bicg',
'bicgstab', 'cg', 'cgs', 'csc_matrix', 'csr_matrix', 'dsolve',
'eigen', 'eigen_symmetric', 'factorized', 'gmres', 'interface',
'isolve', 'iterative', 'lgmres', 'linsolve', 'lobpcg', 'lsqr',
'minres', 'np', 'qmr', 'speigs', 'spilu', 'splu', 'spsolve', 'svd',
'test', 'umfpack', 'use Solver', 'utils', 'warnings']
```

### 11.3.1 Sparse Direct Solvers

* default solver: SuperLU 4.0

- included in SciPy

- real and complex systems

- both single and double precision

* optional: umfpack

- real and complex systems

- double precision only

- recommended for performance

- wrappers now live in `scikits.umfpack`

- check-out the new `scikits.suitesparse` by Nathaniel Smith

#### Examples

* import the whole module, and see its docstring:

```python
>>> from scipy sparse linalg import dsolve
>>> help(dsolve)
```

* both `superlu` and `umfpack` can be used (if the latter is installed) as follows:

- prepare a linear system:

```python
>>> import numpy as np
>>> from scipy import sparse
>>> mtx = sparse.spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]], [0, 1], 5, 5)
>>> mtx.todense()
matrix([[ 1,  5,  0,  0,  0],
        [ 0,  2,  8,  0,  0],
        [ 0,  0,  3,  9,  0],
        [ 0,  0,  0,  4, 10],
        [ 0,  0,  0,  0,  5]])
>>> rhs = np.array([1, 2, 3, 4, 5], dtype= np.float32)
```
- solve as single precision real:

```python
>>> mtx1 = mtx.astype(np.float32)
>>> x = dsolve.spsolve(mtx1, rhs, use_umfpack=False)
>>> print(x)
[ 106. -21.    5.5  -1.5    1. ]
>>> print("Error: %s" % (mtx1 * x - rhs))
Error: [ 0. 0. 0. 0. 0.]
```

- solve as double precision real:

```python
>>> mtx2 = mtx.astype(np.float64)
>>> x = dsolve.spsolve(mtx2, rhs, use_umfpack=True)
>>> print(x)
[ 106. -21.    5.5  -1.5    1. ]
>>> print("Error: %s" % (mtx2 * x - rhs))
Error: [ 0. 0. 0. 0. 0.]
```

- solve as single precision complex:

```python
>>> mtx1 = mtx.astype(np.complex64)
>>> x = dsolve.spsolve(mtx1, rhs, use_umfpack=False)
>>> print(x)
[ 106.0+0.j  -21.0+0.j   5.5+0.j  -1.5+0.j   1.0+0.j]
>>> print("Error: %s" % (mtx1 * x - rhs))
Error: [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
```

- solve as double precision complex:

```python
>>> mtx2 = mtx.astype(np.complex128)
>>> x = dsolve.spsolve(mtx2, rhs, use_umfpack=True)
>>> print(x)
[ 106.0+0.j  -21.0+0.j   5.5+0.j  -1.5+0.j   1.0+0.j]
>>> print("Error: %s" % (mtx2 * x - rhs))
Error: [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
```

```python
#"
Solve a linear system
========================

Construct a 1000x1000 lil_matrix and add some values to it, convert it
to CSR format and solve A x = b for x: and solve a linear system with a
direct solver.
#"
import numpy as np
import scipy sparse as sps
from matplotlib import pyplot as plt
from scipy sparse linalg.dsolve import linsolve

rand = np.randomrand

mtx = sps lil_matrix((1000, 1000), type = np float64)
mtx[0, :100] = rand(100)
mtx[1, 100:200] = mtx[0, :100]
mtx.setdiag(rand(1000))

plt clf()
plt spy(mtx, marker = '.', markersize = 2)
plt show()

mtx = mtx.tocsr()
rhs = rand(1000)

x = linsolve spsolve(mtx, rhs)
```
```python
print('rezidual: %r' % np.lINALG.norm(mtx * x - rhs))
```

* examples/direct_solve.py

### 11.3.2 Iterative Solvers

* the `isolve` module contains the following solvers:

- bicg (BIConjugate Gradient)

- bicgstab (BIConjugate Gradient STABilized)

- cg (Conjugate Gradient) - symmetric positive definite matrices only

- cgs (Conjugate Gradient Squared)

- gmres (Generalized Minimal RESidual)

- minres (MINimum RESidual)

- qmr (Quasi-Minimal Residual)

#### Common Parameters

* mandatory:

A [{sparse matrix, dense matrix, LinearOperator}] The N-by-N matrix of the linear system.

b [{array, matrix}] Right hand side of the linear system. Has shape (N,) or (N,1).

* optional:

x0 [{array, matrix}] Starting guess for the solution.

tol [float] Relative tolerance to achieve before terminating.

maxiter [integer] Maximum number of iterations. Iteration will stop after maxiter steps even if the specified tolerance has not been achieved.

M [{sparse matrix, dense matrix, LinearOperator}] Preconditioner for A. The preconditioner should approximate the inverse of A. Effective preconditioning dramatically improves the rate of convergence, which implies that fewer iterations are needed to reach a given error tolerance.

callback [function] User-supplied function to call after each iteration. It is called as callback(xk), where xk is the current solution vector.

#### LinearOperator Class

```python
from scipy sparse linalg interface import LinearOperator
```

* common interface for performing matrix vector products

* useful abstraction that enables using dense and sparse matrices within the solvers, as well as *matrix-free* solutions

* has shape and `matvec()` (+ some optional parameters)

* example:

```python
>>> import numpy as np
>>> from scipyett.sparse.linalg import LinearOperator
>>> def mv(v):
...     return np.array([2*v[0], 3*v[1]])
...
>>> A = LinearOperator((2, 2), matvec=mv)
```
```python
>>> A
<2x2 LinearOperator with unspecified dtype>
>>> A.matvec(np.ones(2))
array([ 2., 3.])
>>> A * np.ones(2)
array([ 2., 3.])
```

## A Few Notes on Preconditioning

* problem specific

* often hard to develop

* if not sure, try ILU

- available in `dsolve` as `spilu()`

## 11.3.3 Eigenvalue Problem Solvers

### The `eigen` module

* `arpack` * a collection of Fortran77 subroutines designed to solve large scale eigenvalue problems

* `lobpcg` (Locally Optimal Block Preconditioned Conjugate Gradient Method) * works very well in combination with `PyAMG` * example by Nathan Bell:

```python
"""
Compute eigenvectors and eigenvalues using a preconditioned eigensolver
====================================================================

In this example Smoothed Aggregation (SA) is used to precondition
the LOBPCG eigensolver on a two-dimensional Poisson problem with
Dirichlet boundary conditions.
"""

import spipy
from spipy.sparse.linalg import lobpcg

from pyamg import smoothed_aggregation_solver
from pyamg.gallery import poisson

N = 100
K = 9
A = poisson((N,N), format='CSR')

# create the AMG hierarchy
ml = smoothed_aggregation_solver(A)

# initial approximation to the K eigenvectors
X = spipyrand(A.shape[0], K)

# preconditioner based on ml
M = ml.aspreconditioner()

# compute eigenvalues and eigenvectors with LOBPCG
W,V = lobpcg(A, X, M=M, tol=1e-8, largest=False)

# plot the eigenvectors
import pylab
```
```python
pylab.figure(figsize=(9,9))

for i in range(K):
    LABLab.subplot(3, 3, i+1)
    LABLab.title('Eigenvector %d' % i)
    LABLab.pcolor(V[:i].reshape(N,N))
    LABLab.axis('equal')
    LABLab.axis('off')
pylab.show()
```

- examples/pyamg_with_lobpcg.py

* example by Nils Wagner:

- examples/lobpcg_sakurai.py

• output:

```bash
$ python examples/lobpcg_sakurai.py
Results by LOBPCG for n=2500

[ 0.06250083  0.06250028  0.06250007]

Exact eigenvalues

[ 0.06250005  0.0625002  0.06250044]

Elapsed time 7.01
```

<figure><img src="image_136.png" /></figure>
## 11.4 Other Interesting Packages

* **PyAMG**

- algebraic multigrid solvers

- http://code.google.com/p/pyamg

* **Pysparse**

- own sparse matrix classes

- matrix and eigenvalue problem solvers

- http://pysparse.sourceforge.net/
# Image manipulation and processing using Numpy and Scipy

**Authors:** Emmanuelle Gouillart, Gaël Varoquaux

This section addresses basic image manipulation and processing using the core scientific modules NumPy and SciPy. Some of the operations covered by this tutorial may be useful for other kinds of multidimensional array processing than image processing. In particular, the submodule `scipy.ndimage` provides functions operating on n-dimensional NumPy arrays.

**See also:**

For more advanced image processing and image-specific routines, see the tutorial *Scikit-image: image processing*, dedicated to the `skimage` module.

## Image = 2-D numerical array

(or 3-D: CT, MRI, 2D + time; 4-D, ...)

Here, `image == Numpy array np.array`

**Tools used in this tutorial:**

* ` numpy`: basic array manipulation

* `scipy: scipy.ndimage` submodule dedicated to image processing (n-dimensional images). See the documentation:

>>> from scipy import ndimage

**Common tasks in image processing:**

* Input/Output, displaying images

* Basic manipulations: cropping, flipping, rotating, ...

* Image filtering: denoising, sharpening

* Image segmentation: labeling pixels corresponding to different objects

* Classification

* Feature extraction

* Registration

* ...
## Chapters contents

* Opening and writing to image files

* Displaying images

* Basic manipulations

- Statistical information

- Geometrical transformations

* Image filtering

- Blurring/smoothing

- Sharpening

- Denoising

- Mathematical morphology

* Feature extraction

- Edge detection

- Segmentation

* Measuring objects properties: `ndimage.measurements`

## 12.1 Opening and writing to image files

Writing an array to a file:

```
Displaying Lena
==================

Small example to plot lena.
"
"

from scipy import misc
l = misc.lena()
misc.imsave('lena.png', l) # uses the Image module (PIL)

import matplotlib.pyplot as plt
plt.imshow(l)
plt.show()
```

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

Creating a numpy array from an image file:
```python
>>> from scipy import misc
>>> lena = misc.lena()
>>> misc.imsave('lena.png', lena)  # First we need to create the PNG file
>>> lena = misc.imread('lena.png')
>>> type(lena)
<... 'numpy.ndarray'>
>>> lena.shape, lena.dtype
((512, 512), type('uint8'))
```

dtype is `uint8` for 8-bit images (0-255)

Opening raw files (camera, 3-D images)

```python
>>> lena.tofile('lena/raw')  # Create raw file
>>> lena_from_raw = np.fromfile('lena/raw', type(npuint8)
>>> lena_from_raw.shape
(262144,.)
>>> lena_from_raw.shape = (512, 512)
```

Need to know the shape and `dtype` of the image (how to separate data bytes).

For large data, use `np.mmap` for memory mapping:

```python
>>> lena_mmap = np.mmap('lena/raw', type(npuint8), shape=(512, 512))
```

(data are read from the file, and not loaded into memory)

Working on a list of image files

```python
>>> for i in range(10):
...     im = np.random.random_integers(0, 255, 10000).reshape((100, 100))
...     misc.imsave('random_%02d.png' % i, im)
>>> from glob import glob
>>> filelist = glob('random*.png')
>>> filelist.sort()
```

## 12.2 Displaying images

Use ` matplotlib` and `imshow` to display an image inside a ` matplotlib` figure:

```python
>>> l = misc.lena()
>>> import matplotlib.pyplot as plt
>>> plt.imshow(l, cmap= plt.cm .gray)
< matplotlib .image .AxesImage object at 0x... >
```

Increase contrast by setting `min` and `max` values:

```python
>>> plt.imshow(l, cmap= plt.cm .gray, vmin=30, vmax=200)
< matplotlib .image .AxesImage object at 0x... >
>>> # Remove axes and ticks
>>> plt .axis('off')
(-0.5, 511.5, 511.5, -0.5)
```

Draw contour lines:

```python
>>> plt.contour(l, [60, 211])
< matplotlib .contour .QuadContourSet ... >
```

For fine inspection of intensity variations, use `interpolation='nearest'`:
<figure><img src="image_138.png" /></figure>

```python
>>> plt.imshow(l[200:220, 200:220], cmap= plt.cm .gray)
< matplotlib .image .AxesImage object at 0x... >
>>> plt.imshow(l[200:220, 200:220], cmap= plt.cm .gray, interpolation='nearest')
< matplotlib .image .AxesImage object at 0x... >
```

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

### 3-D visualization: Mayavi

See *3D plotting with Mayavi*.

* Image plane widgets

* Isosurfaces

* ...

## 12.3 Basic manipulations

Images are arrays: use the whole ` numpy ` machinery.
<figure><img src="image_140.png" /></figure>

<table><tr><td>0</td><td>1</td><td>2</td></tr><tr><td>3</td><td>4</td><td>5</td></tr><tr><td>6</td><td>7</td><td>8</td></tr></table>

```python
>>> lena = misc.lena()
>>> lena[0, 40]
166
>>> # Slicing
>>> lena[10:13, 20:23]
array([[158, 156, 157],
      [157, 155, 155],
      [157, 157, 158]])
>>> lena[100:120] = 255
>>> 
>>> lx, ly = lena.shape
>>> X, Y = np.o grid [0:lx, 0:ly]
>>> mask = (X - lx / 2) ** 2 + (Y - ly / 2) ** 2 > lx * ly / 4
>>> # Masks
>>> lena [mask] = 0
>>> # Fancy indexing
>>> lena[ range (400), range (400)] = 255
```

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

### 12.3.1 Statistical information

```python
>>> lena = misc.lena()
>>> lena .mean()
124.04678344726562
>>> lena .max() , lena .min()
(245, 25)
```
np.histogram

### Exercise

* Open as an array the scikit-image logo (http://scikit-image.org/_static/img/logo.png), or an image that you have on your computer.

* Crop a meaningful part of the image, for example the python circle in the logo.

* Display the image array using ` matplotlib`. Change the interpolation method and zoom to see the difference.

* Transform your image to greyscale

* Increase the contrast of the image by changing its minimum and maximum values. **Optional**: use `scipy.stats.scoreatpercentile` (read the docstring!) to saturate 5% of the darkest pixels and 5% of the lightest pixels.

* Save the array to two different file formats (png, jpg, tiff)

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

scikits-image
image processing in python

### 12.3.2 Geometrical transformations

```python
>>> lena = misc.lena()
>>> lx, ly = lena.shape
>>> # Cropping
>>> crop_lena = lena[lx / 4: - lx / 4, ly / 4: - ly / 4]
>>> # up <-> down flip
>>> flip_ud_lena = np flipped(lena)
>>> # rotation
>>> rotate_lena = ndimage Rotate(lena, 45)
>>> rotate_lena_noreshape = ndimage Rotate(lena, 45, reshape=False)
```

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

## 12.4 Image filtering

**Local filters:** replace the value of pixels by a function of the values of neighboring pixels.

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

<figure><img src="image_144.png" /></figure>
### 12.4.1 Blurring/smoothing

**Gaussian filter** from `scipy.ndimage`:

```python
>>> from scipy import misc
```

```python
>>> lena = misc.lena()
```

```python
>>> blurred_lena = ndimage.gaussian_filter(lena, sigma=3)
```

```python
>>> very_blurred = ndimage.gaussian_filter(lena, sigma=5)
```

**Uniform filter**

```python
>>> local_mean = ndimage uniform_filter(lena, size=11)
```

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

### 12.4.2 Sharpening

Sharpen a blurred image:

```python
>>> from scipy import misc
```

```python
>>> lena = misc.lena()
```

```python
>>> blurred_l = ndimage.gaussian_filter(lena, 3)
```

increase the weight of edges by adding an approximation of the Laplacian:

```python
>>> filter.blurred_l = ndimage.gaussian_filter(blurred_l, 1)
```

```python
>>> alpha = 30
```

```python
>>> sharpened = blurred_l + alpha * (blurred_l - filter.blurred_l)
```

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

### 12.4.3 Denoising

Noisy lena:

```python
>>> from scipy import misc
```

```python
>>> l = misc.lena()
```

```python
>>> l = l[230:310, 210:350]
```

```python
>>> noisy = l + 0.4 * l.std() * np.random.random(l.shape)
```
A **Gaussian filter** smoothes the noise out... and the edges as well:

```python
>>> gauss_denoised = ndimage.gaussian_filter(noisy, 2)
```

Most local linear isotropic filters blur the image (`ndimage uniform_filter`)

A **median filter** preserves better the edges:

```python
>>> med_denoised = ndimage.median_filter(noisy, 3)
```

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

Median filter: better result for straight boundaries (**low curvature**):

```python
>>> im = np(zeros((20, 20)))
```

```python
>>> im[5:-5, 5:-5] = 1
```

```python
>>> im = ndimage.distance_transform_bf(im)
```

```python
>>> im_noise = im + 0.2 * np.randomrandn(*im.shape)
```

```python
>>> im_med = ndimage.median_filter(im_noise, 3)
```

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

Other rank filter: `ndimage.maximmfilter`, `ndimage.percentile_filter`

Other local non-linear filters: Wiener (`scipy.signal.wiener`), etc.

### Non-local filters

#### Exercise: denoising

* Create a binary image (of 0s and 1s) with several objects (circles, ellipses, squares, or random shapes).

* Add some noise (e.g., 20% of noise)

* Try two different denoising methods for denoising the image: gaussian filtering and median filtering.

* Compare the histograms of the two different denoised images. Which one is the closest to the histogram of the original (noise-free) image?

#### See also:

More denoising filters are available in `skimage.denoising`, see the *Scikit-image: image processing* tutorial.

### 12.4.4 Mathematical morphology

See *wikipedia* for a definition of 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.

#### Structuring element:
```python
>>> el = ndimage(generate_binary_structure(2, 1)
>>> el
array([[False, True, False],
      [ True, True, True],
      [False, True, False]], dtype=bool)
>>> el.astype(np.int)
array([[0, 1, 0],
      [1, 1, 1],
      [0, 1, 0]])
```

<figure><img src="image_149.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, 0, 1, 1, 1, 0, 0],
      [0, 0, 1, 1, 1, 0, 0],
      [0, 0, 1, 1, 1, 0, 0],
      [0, 0, 1, 1, 1, 0, 0],
      [0, 0, 0, 0, 0, 0, 0]])
>>> ndimage-binary_erosion(a).astype(a.dtype)
array([[0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 1, 0, 0, 0],
      [0, 0, 0, 1, 0, 0, 0],
      [0, 0, 0, 1, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0]])
>>> #Erosion removes objects smaller than the structure
>>> ndimage-binary_erosion(a, structure(np.ones((5,5))).astype(a.dtype)
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, 0, 0, 0, 0, 0, 0]])
```
Erosion

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

Dilation

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

Opening

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

Closing

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

**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.]])
>>> ndimage binary dilation(a).astype(a da type)
array([[ 0.,  0.,  0.,  0.,  0.],
      [ 0.,  0.,  1.,  0.,  0.],
      [ 0.,  1.,  1.,  1.,  0.],
      [ 0.,  0.,  1.,  0.,  0.],
      [ 0.,  0.,  0.,  0.,  0.]])
```

Also works for grey-valued images:

```python
>>> np.random.seed(2)
>>> im = np zeros((64, 64))
>>> x, y = (63* np random random((2, 8))).astype(np int)
>>> im[x, y] = np arange(8)
>>> bigger_points = ndimage grey dilation(im, size=(5, 5), structure = np ones((5, 5)))
>>> square = np zeros((16, 16))
>>> square[4:-4, 4:-4] = 1
>>> dist = ndimage distance transform_bf square
>>> dilate_dist = ndimage grey dilation(dist, size=(3, 3), \
...     structure = np ones((3, 3)))
```

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

**Opening:** erosion + dilation:

```python
>>> a = np zeros((5,5), da type = 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]])
>>> # Opening removes small objects
>>> ndimage binary_opening(a, structure = np ones((3,3))).astype(np int)
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, 0]]
>>> # Opening can also smooth corners
>>> ndimage-binary_opening(a).astype(np.int)
array([[0, 0, 0, 0, 0],
      [0, 0, 1, 0, 0],
      [0, 1, 1, 1, 0],
      [0, 0, 1, 0, 0],
      [0, 0, 0, 0, 0]])
```

**Application:** remove noise:

```
>>> square = np.zeros((32, 32))
>>> square[10:-10, 10:-10] = 1
>>> np.random.seed(2)
>>> x, y = (32* np.random.random((2, 20))).astype(np.int)
>>> square[x, y] = 1
>>> open_square = ndimage-binary_opening.square)
>>> eroded_square = ndimage-binary_erosion.square)
>>> reconstruction = ndimage-binary Propagation(eroded_square, mask=square)
```

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

**Closing:** dilation + erosion

Many other mathematical morphology operations: hit and miss transform, tophat, etc.

# 12.5 Feature extraction

## 12.5.1 Edge detection

Synthetic data:

```
>>> im = np.zeros((256, 256))
>>> im[64:-64, 64:-64] = 1
>>> im = ndimage Rotate(im, 15, mode='constant')
>>> im = ndimage.gaussian_filter(im, 8)
```

Use a **gradient operator (Sobel)** to find high intensity variations:

```
>>> sx = ndimage.sobel(im, axis=0, mode='constant')
>>> sy = ndimage.sobel(im, axis=1, mode='constant')
>>> sob = np.hypot(sx, sy)
```

<figure><img src="image_156.png" /></figure>
12.5.2 Segmentation

• Histogram-based segmentation (no spatial information)

>>> n = 10
>>> l = 256
>>> im = np(zeros((1, 1)))
>>> np.random.seed(1)
>>> points = l*np.random.random((2, n**2))
>>> im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
>>> im = ndimage.gaussian_filter(im, sigma=1/(4.*n))
>>> mask = (im > im.mean()).astype(np.float)
>>> mask += 0.1 * im
>>> img = mask + 0.2*np.random.rndn(*mask.shape)
>>> hist, bin_edges = np.histogram(img, bins=60)
>>> binCenters = 0.5*(bin_edges[:-1] + bin_edges[1:])
>>> binaryImg = img > 0.5

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

Use mathematical morphology to clean up the result:

>>> # Remove small white regions
>>> openImg = ndimage-binary_opening(binaryImg)
>>> # Remove small black hole
>>> closeImg = ndimage-binary_closing(openImg)

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

Exercise

Check that reconstruction operations (erosion + propagation) produce a better result than open-
ing/closing:

>>> erodedImg = ndimage-binary_erosion(binaryImg)
>>> reconstructImg = ndimage-binary Propagation(erodedImg, mask=binaryImg)
>>> tmp = np(logical_not(reconstructImg)
>>> eroded_tmp = ndimage-binary_erosion(tmp)
>>> reconstruct_final = np(logical_not(ndimage-binary-propagation(eroded_tmp, mask(tmp)))
>>> np abs mask - closeImg).mean()
0.00727836...
>>> np abs mask - reconstruct_final).mean()
0.00059502...
## Exercise

Check how a first denoising step (e.g. with a median filter) modifies the histogram, and check that the resulting histogram-based segmentation is more accurate.

**See also:**

More advanced segmentation algorithms are found in the `scikit-image`: see *Scikit-image: image processing*.

**See also:**

Other Scientific Packages provide algorithms that can be useful for image processing. In this example, we use the spectral clustering function of the `scikit-learn` in order to segment glued objects.

```python
>>> from sklearn feature_extraction import image
>>> from sklearn.cluster import spectral_clustering
>>> l = 100
>>> x, y = np indices((1, 1))
>>> center1 = (28, 24)
>>> center2 = (40, 50)
>>> center3 = (67, 58)
>>> center4 = (24, 70)
>>> radius1, radius2, radius3, radius4 = 16, 14, 15, 14
>>> circle1 = (x - center1[0])*2 + (y - center1[1])*2 < radius1*2
>>> circle2 = (x - center2[0])*2 + (y - center2[1])*2 < radius2*2
>>> circle3 = (x - center3[0])*2 + (y - center3[1])*2 < radius3*2
>>> circle4 = (x - center4[0])*2 + (y - center4[1])*2 < radius4*2
>>> # 4 circles
>>> img = circle1 + circle2 + circle3 + circle4
>>> mask = img .astype(bool)
>>> img = img .astype(float)
>>> img += 1 + 0.2* np random .randn(*img .shape)
>>> # Convert the image into a graph with the value of the gradient on
>>> # the edges.
>>> graph = image .img_to_graph(img, mask =mask)
>>> # Take a decreasing function of the gradient: we take it weakly
>>> # dependant from the gradient the segmentation is close to a voronoi
>>> graph.data = np exp(-graph.data/graph.data .std())
>>> labels = spectral_clustering(graph, n clusters=4, eigen solver='arpack')
>>> label_im = -np ones mask .shape)
>>> label_im[mask] = labels
```
<figure><img src="image_159.png" /></figure>

## 12.6 Measuring objects properties: ndimage.measurements

Synthetic data:

```python
>>> n = 10
>>> 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 = ndimage.gaussian_filter(im, sigma=1/(4.*n))
>>> mask = im > im.mean()
```

* **Analysis of connected components**

Label connected components: ndimage.label:

```python
>>> label_im, nb labels = ndimage.label mask
>>> nb labels # how many regions?
16
>>> plt.imshow label im
< matplotlib image AxesImage object at 0x... >
```

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

Compute size, mean_value, etc. of each region:

```python
>>> sizes = ndimage.sum mask, label im, range nb labels + 1)
>>> mean vals = ndimage.sum(im, label im, range(1, nb labels + 1))
```

Clean up small connected components:

```python
>>> mask_size = sizes < 1000
>>> remove pixel = mask_size label im
```
```python
>>> removePixel.shape
(256, 256)
>>> label_im[removePixel] = 0
>>> plt.imshow(label_im)
< matplotlib.image.AxesImage object at 0x... >
```

Now reassign labels with np.searchsorted:

```python
>>> labels = np.unique(label_im)
>>> label_im = np.searchsorted(labels, label_im)
```

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

Find region of interest enclosing object:

```python
>>> slice_x, slice_y = ndimage.find_objects(label_im==4)[0]
>>> roi = im[slice_x, slice_y]
>>> plt.imshow(roi)
< matplotlib.image.AxesImage object at 0x... >
```

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

Other spatial measures: ndimage.center_of, ndimage.maxim...

Can be used outside the limited scope of segmentation applications.

Example: block mean:

```python
>>> from scipy import misc
>>> l = misc.lena()
```
```python
>>> sx, sy = l.shape
>>> X, Y = np.o grid[0:sx, 0:sy]
>>> regions = (sy // 6) * (X // 4) + (Y // 6)  # note that we use broadcasting
>>> block_mean = ndimage .mean(1, labels=regions, index =np .arange(1,
...       regions .max() +1))
>>> block_mean .shape = (sx // 4, sy // 6)
```

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

When regions are regular blocks, it is more efficient to use stride tricks (Example: fake dimensions with strides).

Non-regularly-spaced blocks: radial mean:

```python
>>> sx, sy = l.shape
>>> X, Y = np .o grid[0:sx, 0:sy]
>>> r = np .hypot(X - sx /2, Y - sy /2)
>>> rbin = (20* r/r .max()).astype(np .int)
>>> radial_mean = ndimage .mean(1, labels=rbin, index =np .arange(1, rbin .max() +1))
```

• Other measures

Correlation function, Fourier/wavelet spectrum, etc.

One example with mathematical morphology: granulometry

```python
>>> def disk_structure(n):
...     struct = np .zeros((2 * n + 1, 2 * n + 1))
...     x, y = np .indices((2 * n + 1, 2 * n + 1))
...     mask = (x - n)**2 + (y - n)**2 <= n**2
...     struct [mask] = 1
...     return struct .astype(np .bool)
...
>>> def granulometry(data, sizes = None ):
...     s = max(data .shape)
...     if sizes == None :
...         sizes = range(1, s/2, 2)
...     granulo = [ndimage .binary_opening(data, \
...             structure=disk_structure(n)).sum() for n in sizes]
...     return granulo
...
>>>
```
<figure><img src="image_164.png" /></figure>

```python
>>> np.random.seed(1)
>>> n = 10
>>> 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 = ndimage.gaussian_filter(im, sigma=1/(4.*n))
>>> mask = im > im.mean()
>>> granulo = granulometry mask, sizes = np.arange(2, 19, 4))
```

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

**See also:**

More on image-processing:

* The chapter on *Scikit-image*
* Other, more powerful and complete modules: OpenCV (Python bindings), CellProfiler, ITK with Python bindings
# Mathematical optimization: finding minima of functions

**Authors:** Gaël Varoquaux

Mathematical optimization deals with the problem of finding numerically minimums (or maximums or zeros) of a function. In this context, the function is called *cost function*, or *objective function*, or *energy*.

Here, we are interested in using `scipy optimize` for black-box optimization: we do not rely on the mathematical expression of the function that we are optimizing. Note that this expression can often be used for more efficient, non black-box, optimization.

## Prerequisites

* Numpy, Scipy

* IPython

* matplotlib

## References

Mathematical optimization is very ... mathematical. If you want performance, it really pays to read the books:

* *Convex Optimization* by Boyd and Vandenberghe (pdf available free online).

* *Numerical Optimization*, by Nocedal and Wright. Detailed reference on gradient descent methods.

* *Practical Methods of Optimization* by Fletcher: good at hand-waving explanations.
## Chapters contents

* Knowing your problem

- Convex versus non-convex optimization

- Smooth and non-smooth problems

- Noisy versus exact cost functions

- Constraints

* A review of the different optimizers

- Getting started: 1D optimization

- Gradient based methods

- Newton and quasi-newton methods

- Gradient-less methods

- Global optimizers

* Practical guide to optimization with scipy

- Choosing a method

- Making your optimizer faster

- Computing gradients

- Synthetic exercices

* Special case: non-linear least-squares

- Minimizing the norm of a vector function

- Curve fitting

* Optimization with constraints

- Box bounds

- General constraints

## 13.1 Knowing your problem

Not all optimization problems are equal. Knowing your problem enables you to choose the right tool.

### Dimensionality of the problem

The scale of an optimization problem is pretty much set by the *dimensionality of the problem*, i.e. the number of scalar variables on which the search is performed.
### 13.1.1 Convex versus non-convex optimization

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

Optimizing convex functions is easy. Optimizing non-convex functions can be very hard.

It can be proven that for a convex function a local minimum is also a global minimum. Then, in some sense, the minimum is unique.

### 13.1.2 Smooth and non-smooth problems

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

Optimizing smooth functions is easier (true in the context of *black-box* optimization, otherwise *Linear Programming* is an example of methods which deal very efficiently with piece-wise linear functions).
### 13.1.3 Noisy versus exact cost functions

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

#### Noisy gradients

Many optimization methods rely on gradients of the objective function. If the gradient function is not given, they are computed numerically, which induces errors. In such situation, even if the objective function is not noisy,

### 13.1.4 Constraints

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

## 13.2 A review of the different optimizers

### 13.2.1 Getting started: 1D optimization

Use `scipyannyaHallidomine`. `scipyannyaHallidomine` to minimize 1D functions. It combines a bracketing strategy with a parabolic approximation.
**Brent's method on a quadratic function:** it converges in 3 iterations, as the quadratic approximation is then exact.

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

**Brent's method on a non-convex function:** note that the fact that the optimizer avoided the local minimum is a matter of luck.

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

```python
>>> from scipy import optimize
>>> def f(x):
...     return -np.exp(-(x - .7)**2)
>>> x_min = optimize.brent(f)  # It actually converges in 9 iterations!
>>> x_min
0.699999999...
>>> x_min - .7
-2.1605...e-10
```

Brent's method can also be used for optimization constrained to an interval using `scipy/optimize.fminbound()`

In `scipy.0.11`, `scipy.optimize.minimize_scalar()` gives a generic interface to 1D scalar minimization

### 13.2.2 Gradient based methods

#### Some intuitions about gradient descent

Here we focus on **intuitions**, not code. Code will follow.

**Gradient descent** basically consists in taking small steps in the direction of the gradient, that is the direction of the **steepest descent**.

<table><caption>Table 13.1: Fixed step gradient descent</caption><tr><td><strong>A well-conditionned quadratic function.</strong></td><td><img src="image_0_0_0.png" alt="A well-conditionned quadratic function. The function has a local minimum at the point (0,0). The x-axis is on the left and the y-axis is on the top. The function has a local maximum at (0,0). The x-axis is on the left and the y-axis is on the top. The function is continuous and has a single minimum (the local minimum at (0,0)). The gradient at (0,0) is (0,1). The value of the function at (0,0) is 0."></td><td><img src="image_0_0_1.png" alt="A graph of a function. The function has a local maximum at (0,0). The x-axis is on the left and the y-axis is on the top. The function is continuous and has a single minimum (the local maximum at (0,0)). The gradient at (0,0) is (0,1). The value of the function at (0,0) is 0."></td></tr><tr><td><strong>An ill-conditionned quadratic function.</strong><br>The core problem of gradient-methods on ill-conditioned problems is that the gradient tends not to point in the direction of the minimum.</td><td><img src="image_0_1_0.png" alt="An ill-conditionned quadratic function. The function has a local maximum at the point (0,0). The x-axis is on the left and the y-axis is on the top. The function has a local minimum at (0,0). The x-axis is on the left and the y-axis is on the top. The function is continuous and has a single maximum (the local maximum at (0,0)). The gradient at (0,0) is (0,0). The value of the function at (0,0) is 0."></td><td><img src="image_0_1_1.png" alt="An ill-conditionned quadratic function. The function has a local maximum at the point (0,0). The x-axis is on the left and the y-axis is on the top. The function has a local minimum at (0,0). The x-axis is on the left and the y-axis is on the top. The function is continuous and has a single minimum (the local minimum at (0,0)). The gradient at (0,0) is (0,0). The value of the function at (0,0) is 0."></td></tr></table>

We can see that very anisotropic (ill-conditionned) functions are harder to optimize.
## Take home message: conditioning number and preconditioning

If you know natural scaling for your variables, prescale them so that they behave similarly. This is related to preconditioning.

Also, it clearly can be advantageous to take bigger steps. This is done in gradient descent code using a line search.

<table><caption>Table 13.2: Adaptive step gradient descent</caption><tr><td>A well-conditionned quadratic function.</td><td></td><td></td></tr><tr><td>An ill-conditionned quadratic function.</td><td></td><td></td></tr><tr><td>An ill-conditionned non-quadratic function.</td><td></td><td></td></tr><tr><td>An ill-conditionned very non-quadratic function.</td><td></td><td></td></tr></table>

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

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

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

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

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

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

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

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

The more a function looks like a quadratic function (elliptic iso-curves), the easier it is to optimize.

## Conjugate gradient descent

The gradient descent algorithms above are toys not to be used on real problems.

As can be seen from the above experiments, one of the problems of the simple gradient descent algorithms, is that it tends to oscillate across a valley, each time following the direction of the gradient, that makes it cross
the valley. The conjugate gradient solves this problem by adding a *friction* term: each step depends on the two last values of the gradient and sharp turns are reduced.

<table><caption>Table 13.3: Conjugate gradient descent</caption><tr><td>An ill-conditionned non-quadratic function.</td><td></td><td></td></tr><tr><td>An ill-conditionned very non-quadratic function.</td><td></td><td></td></tr></table>

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

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

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

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

Methods based on conjugate gradient are named with 'cg' in scipy. The simple conjugate gradient method to minimize a function is `scipyfmin_cg(f, [2, 2])`

```python
>>> def f(x): # The rosenbrock function
...     return .5*(1 - x[0]**2 + (x[1] - x[0]**2)**2)
>>> optimize.fmin_cg(f, [2, 2])
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 13
Function evaluations: 120
Gradient evaluations: 30
array([ 0.99998968,  0.99997855])
```

These methods need the gradient of the function. They can compute it, but will perform better if you can pass them the gradient:

```python
>>> def fprime(x):
...     return np.array((-2.*.5*(1 - x[0])) - 4*x[0)*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.fmin_cg(f, [2, 2], fprime=fprime)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 13
Function evaluations: 30
Gradient evaluations: 30
array([ 0.99999199,  0.99998336])
```

Note that the function has only been evaluated 30 times, compared to 120 without the gradient.

### 13.2.3 Newton and quasi-newton methods

#### Newton methods: using the Hessian (2nd differential)

Newton methods use a local quadratic approximation to compute the jump direction. For this purpose, they rely on the 2 first derivative of the function: the *gradient* and the *Hessian*.
**An ill-conditionned quadratic function:**

Note that, as the quadratic approximation is exact, the Newton method is blazing fast

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

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

**An ill-conditionned non-quadratic function:**

Here we are optimizing a Gaussian, which is always below its quadratic approximation. As a result, the Newton method overshoots and leads to oscillations.

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

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

**An ill-conditionned very non-quadratic function:**

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

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

In scipy, the Newton method for optimization is implemented in `scipy-optimize.fmin_ncg()` (cg here refers to that fact that an inner operation, the inversion of the Hessian, is performed by conjugate gradient). `scipy-optimize.fmin_tnc()` can be use for constraint problems, although it is less versatile:

```python
>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0]**2) + (x[1] - x[0]**2)**2
>>> def fprime(x):
...     return np.array((-2.*.5*(1 - x[0]) - 4*x[0)*(x[1] - x[0]**2),
...                      2*(x[1] - x[0]**2)))
>>> optimize.fmin_ncg(f, [2, 2], fprime=fprime)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 9
Function evaluations: 11
Gradient evaluations: 51
Hessian evaluations: 0
array([ 1., 1.])
```

Note that compared to a conjugate gradient (above), Newton's method has required less function evaluations, but more gradient evaluations, as it uses it to approximate the Hessian. Let's compute the Hessian and pass it to the algorithm:

```python
>>> def hessian(x): # Computed with sympy
...     return np.array(((1 - 4*x[1] + 12*x[0]**2, -4*x[0]), (-4*x[0], 2)))
>>> optimize.fmin_ncg(f, [2, 2], fprime=fprime, fhess=hessian)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 9
Function evaluations: 11
Gradient evaluations: 19
Hessian evaluations: 9
array([ 1., 1.])
```

At very high-dimension, the inversion of the Hessian can be costly and unstable (large scale > 250).

Newton optimizers should not be confused with Newton's root finding method, based on the same principles, `scipy-optimize.newton()`.
## Quasi-Newton methods: approximating the Hessian on the fly

**BFGS:** BFGS (Broyden-Fletcher-Goldfarb-Shanno algorithm) refines at each step an approximation of the Hessian.

**An ill-conditionned quadratic function:**

On an exactly quadratic function, BFGS is not as fast as Newton's method, but still very fast.

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

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

**An ill-conditionned non-quadratic function:**

Here BFGS does better than Newton, as its empirical estimate of the curvature is better than that given by the Hessian.

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

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

**An ill-conditionned very non-quadratic function:**

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

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

```python
>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0]**2) + (x[1] - x[0]**2)**2
>>> def fprime(x):
...     return np.array((-2.*.5*(1 - x[0]) - 4*x[0)*(x[1] - x[0]**2),
...                      2*(x[1] - x[0]**2)))
>>> optimize.fmin_bfgs(f, [2, 2], fprime=fprime)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 16
Function evaluations: 24
Gradient evaluations: 24
array([ 1.00000017,  1.00000026])
```

**L-BFGS:** Limited-memory BFGS sits between BFGS and conjugate gradient: in very high dimensions (> 250) the Hessian matrix is too costly to compute and invert. L-BFGS keeps a low-rank version. In addition, the sciPY version, `scipy optimize.fmin_l_bfgs_b()`, includes box bounds:

```python
>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0]**2) + (x[1] - x[0]**2)**2
>>> def fprime(x):
...     return np.array((-2.*.5*(1 - x[0]) - 4*x[0)*(x[1] - x[0]**2),
...                      2*(x[1] - x[0]**2)))
>>> optimize.fmin_l_bfgs_b(f, [2, 2], fprime=fprime)
(array([ 1.00000005,  1.00000009]), 1.4417677473011859e-15, {...})
```

If you do not specify the gradient to the L-BFGS solver, you need to add `approxGrad=1`

### 13.2.4 Gradient-less methods

#### A shooting method: the Powell algorithm

Almost a gradient approach
**An ill-conditionned quadratic function:**
Powell's method isn't too sensitive to local ill-conditionning in low dimensions

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

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

**An ill-conditionned very non-quadratic function:**

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

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

## Simplex method: the Nelder-Mead

The Nelder-Mead algorithms is a generalization of dichotomy approaches to high-dimensional spaces. The algorithm works by refining a simplex, the generalization of intervals and triangles to high-dimensional spaces, to bracket the minimum.

**Strong points:** it is robust to noise, as it does not rely on computing gradients. Thus it can work on functions that are not locally smooth such as experimental data points, as long as they display a large-scale bell-shape behavior. However it is slower than gradient-based methods on smooth, non-noisy functions.

**An ill-conditionned
non-quadratic function:**

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

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

**An ill-conditionned very
non-quadratic function:**

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

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

In `scipyteroptimize.fmin()` implements the Nelder-Mead approach:

```python
>>> def f(x):  # The rosenbrock function
...     return .5*(1 - x[0])*2 + (x[1] - x[0]**2)**2
>>> optimize.fmin(f, [2, 2])
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 46
Function evaluations: 91
array([ 0.99998568,  0.99996682])
```
### 13.2.5 Global optimizers

If your problem does not admit a unique local minimum (which can be hard to test unless the function is convex), and you do not have prior information to initialize the optimization close to the solution, you may need a global optimizer.

#### Brute force: a grid search

`scipyannya.brute()` evaluates the function on a given grid of parameters and returns the parameters corresponding to the minimum value. The parameters are specified with ranges given to `numpy.mgrid`. By default, 20 steps are taken in each direction:

```python
>>> def f(x):  # The rosenbrock function
...     return .5*(1 - x[0]**2 + (x[1] - x[0]**2)**2)
>>> optimize.brute(f, ((-1, 2), (-1, 2)))
array([ 1.00001462,  1.00001547])
```

## 13.3 Practical guide to optimization with scipy

### 13.3.1 Choosing a method

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

**Without knowledge of the gradient**

* In general, prefer `BFGS` (`scipyannya.fmin_bfgs()`) or `L-BFGS` (`scipyannya.fmin_l_bfgs_b()`), even if you have to approximate numerically gradients

* On well-conditioned problems, Powell (`scipyannya.fmin_powell()`) and Nelder-Mead (`scipyannya.fmin()`), both gradient-free methods, work well in high dimension, but they collapse for ill-conditioned problems.

**With knowledge of the gradient**

* `BFGS` (`scipyannya.fmin_bfgs()`) or `L-BFGS` (`scipyannya.fmin_l_bfgs_b()`).

* Computational overhead of `BFGS` is larger than that of `L-BFGS`, itself larger than that of conjugate gradient. On the other side, `BFGS` usually needs less function evaluations than `CG`. Thus conjugate gradient method is better than `BFGS` at optimizing computationally cheap functions.

**With the Hessian**

* If you can compute the Hessian, prefer the Newton method (`scipyannya.fmin_ncg()`).
### If you have noisy measurements

* Use Nelder-Mead (scipy optimize.fmin()) or Powell (scipy optimize.fmin_powell()).

### 13.3.2 Making your optimizer faster

* Choose the right method (see above), do compute analytically the gradient and Hessian, if you can.

* Use preconditionning when possible.

* Choose your initialization points wisely. For instance, if you are running many similar optimizations, warm-restart one with the results of another.

* Relax the tolerance if you don't need precision

### 13.3.3 Computing gradients

Computing gradients, and even more Hessians, is very tedious but worth the effort. Symbolic computation with Sympy may come in handy.

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

A very common source of optimization not converging well is human error in the computation of the gradient. You can use `scipy optimize.checkGrad()` to check that your gradient is correct. It returns the norm of the different between the gradient given, and a gradient computed numerically:

```python
>>> optimize.checkGrad(f, fprime, [2, 2])
2.384185791015625e-07
```

See also `scipy optimize.approx_fprime()` to find your errors.

### 13.3.4 Synthetic exercices

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

#### Exercise: A simple (?) quadratic function

Optimize the following function, using K[0] as a starting point:

```python
np.random.seed(0)
K = np.random.normal(size=(100, 100))

def f(x):
    return np.sum((np.dot(K, x - 1))**2) + np.sum(x**2)**2
```

Time your approach. Find the fastest approach. Why is BFGS not working well?
Exercice: A locally flat minimum

Consider the function $exp(-1/(.1*x**2 + y**2))$. This function admits a minimum in (0, 0). Starting from an initialization at (1, 1), try to get within 1e-8 of this minimum point.

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

13.4 Special case: non-linear least-squares

13.4.1 Minimizing the norm of a vector function

Least square problems, minimizing the norm of a vector function, have a specific structure that can be used in
the Levenberg-Marquardt algorithm implemented in `scipy-optimize.leastsq()`.

Lets try to minimize the norm of the following vectorial function:

```python
>>> def f(x):
...     return np.arctan(x) - np.arctan(np.linspace(0, 1, len(x)))

>>> x0 = np(zeros(10)
>>> optimize.leastsq(f, x0)
(array([ 0.000000,  0.1111111,  0.2222222,  0.3333333,  0.4444444,
       0.5555556,  0.6666667,  0.7777778,  0.8888889,   1.0 ]), 2)
```

This took 67 function evaluations (check it with 'full_output=1'). What if we compute the norm ourselves and
use a good generic optimizer (BFGS):

```python
>>> def g(x):
...     return np.sum(f(x)**2)

>>> optimize.fmin_bfgs(g, x0)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 11
Function evaluations: 144
Gradient evaluations: 12
array([ -7.4...-09,   1.1...e-01,   2.2...e-01,
    3.3...e-01,   4.4...e-01,   5.5...e-01,
    6.6...e-01,   7.7...e-01,   8.8...e-01,
    1.0...e+00])
```

BFGS needs more function calls, and gives a less precise result.

`leastsq` is interesting compared to BFGS only if the dimensionality of the output vector is large, and larger
than the number of parameters to optimize.

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

If the function is linear, this is a linear-algebra problem, and should be solved with `scipy.linalg.lstsq()`.
### 13.4.2 Curve fitting

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

Least square problems occur often when fitting a non-linear to data. While it is possible to construct our optimization problem ourselves, scipy provides a helper function for this purpose: `scipy-optimize-curve-fit()`:

```python
>>> def f(t, omega, phi):
...     return np.cos(omega * t + phi)

>>> x = np.linspace(0, 3, 50)
>>> y = f(x, 1.5, 1) + .1 * np.random.normal(size=50)

>>> optimize.curve_fit(f, x, y)
(array([ 1.51854577,  0.92665541]),
  array([[ 0.00037994, -0.00056796],
        [-0.00056796,  0.00123978]]))
```

**Exercise**

Do the same with `omega = 3`. What is the difficulty?

## 13.5 Optimization with constraints

### 13.5.1 Box bounds

Box bounds correspond to limiting each of the individual parameters of the optimization. Note that some problems that are not originally written as box bounds can be rewritten as such via change of variables.

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

* `scipy-optimize-fminbound()` for 1D-optimization

* `scipy-optimize-fmin_l_bfgs_b()` a *quasi-Newton* method with bound constraints:
```python
>>> def f(x):
...     return np OV 3) **2 + (x[1] - 2) **2
>>> optimize.fmin_l_bfgs_b(f, np.array([0, 0]), approxGrad=1, bounds=(-1.5, 1.5), (-1.5, 1.5))
(array([ 1.5,  1.5]), 1.5811388300841898, {...})
```

### 13.5.2 General constraints

Equality and inequality constraints specified as functions: f(x) = 0 and g(x) < 0.

* `scipy optimize.fmin_slsqp()` Sequential least square programming: equality and inequality con-

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

straints:

```python
>>> def f(x):
...     return np OV 3) **2 + (x[1] - 2) **2
>>> def constraint(x):
...     return np.atleast_1d(1.5 - np.sum(np.abs(x)))
>>> optimize.fmin_slsqp(f, np.array([0, 0]), ieqcons=[constraint, ])
Optimization terminated successfully. (Exit mode 0)
Current function value: 2.47487373504
Iterations: 5
Function evaluations: 20
Gradient evaluations: 5
array([ 1.25004696,  0.24995304])
```

* `scipy optimize.fmin_cobyla()` Constraints optimization by linear approximation: inequality con-straints only:

```python
>>> optimize.fmin_cobyla(f, np.array([0, 0]), cons=constraint)
array([ 1.25009622,  0.24990378])
```

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

The above problem is known as the Lasso problem in statistics, and there exists very efficient solvers for it (for instance in `scikit-learn`). In general do not use generic solvers when specific ones exist.

#### Lagrange multipliers

If you are ready to do a bit of math, many constrained optimization problems can be converted to non-constrained optimization problems using a mathematical trick known as *Lagrange multipliers*.
## Interfacing with C

**Author:** Valentin Haenel

This chapter contains an *introduction* to the many different routes for making your native code (primarily C/C++) available from Python, a process commonly referred to *wrapping*. The goal of this chapter is to give you a flavour of what technologies exist and what their respective merits and shortcomings are, so that you can select the appropriate one for your specific needs. In any case, once you do start wrapping, you almost certainly will want to consult the respective documentation for your selected technique.

### Chapters contents

* Introduction

* Python-C-Api

* Ctypes

* SWIG

* Cython

* Summary

* Further Reading and References

* Exercises

### 14.1 Introduction

This chapter covers the following techniques:

* Python-C-Api

* Ctypes

* SWIG (Simplified Wrapper and Interface Generator)

* Cython

These four techniques are perhaps the most well known ones, of which Cython is probably the most advanced one and the one you should consider using first. The others are also important, if you want to understand the wrapping problem from different angles. Having said that, there are other alternatives out there, but having understood the basics of the ones above, you will be in a position to evaluate the technique of your choice to see if it fits your needs.

The following criteria may be useful when evaluating a technology:

* Are additional libraries required?

* Is the code autogenerated?

* Does it need to be compiled?

* Is there good support for interacting with Numpy arrays?
* Does it support C++?

Before you set out, you should consider your use case. When interfacing with native code, there are usually
two use-cases that come up:

* Existing code in C/C++ that needs to be leveraged, either because it already exists, or because it is faster.

* Python code too slow, push inner loops to native code

Each technology is demonstrated by wrapping the `cos` function from `math.h`. While this is a mostly a trivial example, it should serve us well to demonstrate the basics of the wrapping solution. Since each technique also includes some form of Numpy support, this is also demonstrated using an example where the cosine is computed on some kind of array.

Last but not least, two small warnings:

* All of these techniques may crash (segmentation fault) the Python interpreter, which is (usually) due to bugs in the C code.

* All the examples have been done on Linux, they *should* be possible on other operating systems.

* You will need a C compiler for most of the examples.

## 14.2 Python-C-API

The `Python-C-API` is the backbone of the standard Python interpreter (a.k.a *CPython*). Using this API it is possible to write Python extension module in C and C++. Obviously, these extension modules can, by virtue of language compatibility, call any function written in C or C++.

When using the Python-C-API, one usually writes much boilerplate code, first to parse the arguments that were given to a function, and later to construct the return type.

### Advantages

* Requires no additional libraries

* Lots of low-level control

* Entirely usable from C++

### Disadvantages

* May require a substantial amount of effort

* Much overhead in the code

* Must be compiled

* High maintenance cost

* No forward compatibility across Python versions as C-API changes

* Reference count bugs are easy to create and very hard to track down.

The Python-C-API example here serves mainly for didactic reasons. Many of the other techniques actually depend on this, so it is good to have a high-level understanding of how it works. In 99% of the use-cases you will be better off, using an alternative technique.

Since reference counting bugs are easy to create and hard to track down, anyone really needing to use the Python C-API should read the section about objects, types and reference counts from the official python documentation. Additionally, there is a tool by the name of `cpychecker` which can help discover common errors with reference counting.