.. _sympy_lab:

*****************************************
``sympy``: Symbolic Mathematics in Python
*****************************************

.. centered:: Written by Spencer Lyon (5/8/13)

.. contents::
    :backlinks: none

.. math::

    Ax_t + Bx_t-1 + Cy_t + Dz_t = 0

    E\{Fx_{t+1} + Gx_t + Hx_{t-1} + Jy_{t+1} +
       Ky_t + Lz_{t+1} Kz_t \} = 0


============
Introduction
============

Much of the research you will do in economics will be numerical. This is true for a number of reasons:

1. Many economic problems are characterized by systems of equations that cannot be solved analytically
2. Analytical solutions tend to take longer for a computer to find
3. Powerful routines already exist for finding numerical solutions

However, despite these reasons (and possibly others), there are times when it is very useful to understand how to use a symbolics math package. Here is an incomplete list of these situations:

1. You need to do some algebra to get systems of equations in suitable form for numerical computation
2. You need to make substitutions for some variables and don't want to risk a math error by hand
3. You need to perform non-trivial derivatives or integrals
4. You are trying to understand the domain of a new function
5. Most root finding algorithms search for a single root. Often there are more than one
6. Some solutions may be complex. For numerical methods to solve these problems a real and complex version of the objective function would need to be written.

The go-to symbolics math package in python is ``sympy``. The following quote is taken from the main ``sympy`` website:

    SymPy is a Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible. SymPy is written entirely in Python and does not require any external libraries.

The remainder of this lab will be an introduction to ``sympy``. By the end of the lab you should feel comfortable doing algebra, calculus, changing variables, and much more in ``sympy``.

======
Basics
======

For the rest of this document assume that I have run the following code:

.. ipython:: python

    import sympy as sym
    from sympy import pprint

Exact Calculator
----------------

``sympy`` can be used as a calculator that does exact math. This differs from normal calculators or computer packages in that they can only store values up to machine precision. Machine epsilon is :math:`10^-6` and :math:`10^-16` on 32-bit and 64-bit architectures, respectively. Below are some examples taken from the `sympy tutorial <sym_tutorial>`_.

.. ipython:: python

    a = sym.Rational(1, 2)
    a
    a * 2
    small = sym.Rational(2)**50/sym.Rational(10)**50
    small

This can be very important for at least 2 reasons:

1. You are ensured exact accuracy in the solutions given by ``sympy``
2. When asking for a numerical approximation, the user can specify the precision of the result using either the ``n()`` or ``evalf()`` methods with an optional argument for number of significant digits.

.. ipython:: python

    sym.pi.n()
    sym.pi.evalf(100)
    sym.pi.n(100)

An important symbol in math that numerical software packages have a hard time dealing with is :math:`\infty`. ``sympy`` has a built in representation of :math:`\infty` that is accessed through the class ``sympy.oo``, where ``oo`` is the lower-case letter 'o' two times.

.. todo::

    **Problem 1**: find the number of times the digit 7 appears in the first 10,000 digits of pi. This can be very easy or quite difficult. Clever use of python standard library objects makes this a one-liner.

Defining Variables
------------------

Python is a full object-oriented programming language. As such, everything has a type.

.. ipython:: python

    x = 1.2
    y = 'hi'
    z = 1
    q = lambda x: x ** 2
    r = (1, 2)
    s = [1, 2, 'hi']
    vars = {'x':x, 'y':y, 'z':z, 'q':q, 'r':r, 's':s}
    for i in vars.keys():
        print('Variable %s %s' % (i, type(vars[i])))

Because everything has a type, in order for ``sympy`` to work with objects as if they were math symbols, they need to be defined as such. There are three main ways to do this in sympy:

1. One at at time using the ``sympy.Symbol`` function
2. Multiple at a time using ``sympy.symbols``

.. ipython:: python

    x = sym.Symbol('x')
    y, z = sym.symbols('y z')  # Could have done sym.symbols('y, z')
    vars = {'x':x, 'y':y, 'z':z}
    for i in vars.keys():
        print('Variable %s %s' % (i, type(vars[i])))

The syntax for using ``sympy.Symbol`` is ``var = Symbol('var')``. To use ``sympy.symbols`` you list variable names separated by commas before the equal sign and then repeat the same (with or without commas) inside ``symbols``: ``y, z = sym.symbols('y z')`` or ``y, z = sym.symbols('y, z')``. The key part of each method is to make sure the argument to the ``Symbol`` or ``symbols`` function is a string containing the same contents as the variable name on the left of the equal sign.

There are other ways to use the ``sym.symbols`` function, but for the purposes of this introduction we will simply guide the reader to the `sympy documentation <sym_tutorial>`_.

Expressing Expressions
----------------------

Usually when doing math, one deals not with single symbols or variables, but with expressions or compound systems of inter-dependent symbols. With an understanding of how to define the building blocks of these expressions (symbols), we can now move to learn how to combine symbols and express mathematical relationships.

Basic arithmetic operators work on instances of the ``Symbol`` class just as you would

.. ipython:: python

    a, b, c, x, y, z = sym.symbols('a b c x y z')
    ex = a + a + b**c + x - y / c
    pprint(ex)  # pprint 'pretty_print's the object

The object ``ex`` is an expression, or object made up of multiple ``Symbol`` s. There are many methods (functions that act on an instance of a particular class) that can be called on an expression (note that below the ``[::5]`` is a slice object that gets every 5th item in the list of methods. I do this for space reasons because there are 133 methods):

.. ipython:: python

    filter(lambda x: not str.startswith(x, '_'), dir(ex))[::5]

=======
Algebra
=======

We will demonstrate how methods of a ``sympy`` expression can be manipulated in various algebraic tasks.

Useful Expression Methods
-------------------------

The ``subs`` method can be called in a few ways. ``expr.subs(old, new)`` will replace all instances of ``old`` with ``new`` in ``expr``. Or, the user could pass a dictionary of replacement rules. This is better understood by example (see below).

.. ipython:: python

    pprint(ex)
    pprint(ex.subs(c, 42))  # Replace all instances of c with 42
    pprint(ex.subs({x:10, y:z}))  # Replace x with the number 10 and y with z

The ``expand`` method will expand an expression out as far as possible.

.. ipython:: python

    # Notice pprint is not a requirement
    ((x + y + z) ** 2).expand()

``collect`` will gather all instances of a particular ``Symbol`` or expression "up to powers with rational exponents" (quote ``sympy`` documentataion) and group them together.

.. ipython:: python

    # Some expressions to work with
    ex_1 = a*x**2 + b*x**2 + a*x - b*x + c
    ex_2 = ((x + y + z) ** 3).expand()
    ex_3 = a*sym.exp(2*x) + b*sym.exp(2*x)

    # group all powers of x
    ex_1.collect(x)

    # Do multi-variate collection (note it favors expressions earlier in list)
    ex_2.collect([x, z])

    # collect wrt an expression. Note sympy understands exponential properties
    ex_3.collect(sym.exp(x))

There are many other useful methods available to a ``sympy`` expression. Specifically we encourage you to experiment with ``simplify``, ``trigsimp``, ``series``, ``together``, ``apart``, ``match``, ``diff``, and ``integrate``.

Solving Equations
-----------------

As all good CAS systems should, ``sympy`` has the ability to solve equations and systems of equations. The workhorse function used to do this is ``sympy.solve``.
The following is an excerpt from the docstring of the ``solve`` function:

    Algebraically solves equations and systems of equations.

        Currently supported are:
            - univariate polynomial,
            - transcendental
            - piecewise combinations of the above
            - systems of linear and polynomial equations
            - sytems containing relational expressions.

        Input is formed as:
            f
                - a single Expr or Poly that must be zero,
                - an Equality
                - a Relational expression or boolean
                - iterable of one or more of the above

            symbols (Symbol, Function or Derivative) specified as
                - none given (all free symbols will be used)
                - single symbol
                - denested list of symbols
                  e.g. solve(f, x, y)
                - ordered iterable of symbols
                  e.g. solve(f, [x, y])

By way of example, we will attempt to solve :math:`x^3 + 2x^2 + 8`. As a cubic polynomial, we know that 3 unique roots exist. We will plot the function to get an idea of what it looks like.

.. ipython:: python

    import matplotlib.pyplot as plt
    import numpy as np

    x = np.linspace(-4, 2, 200)  # 200 evenly spaced points on [-4, 2]

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.axhline(linewidth=2, color='k')
    ax.axvline(linewidth=2, color='k')
    plt.title('Cubic polynomial')
    @savefig CubicPolynomial.png width=6in
    ax.plot(x, x**3 + 2*x**2 + 8)

The plot reveals that this particular function only has one real root near -3. Let's see if ``sympy`` can find all three roots for us.

.. ipython:: python

    x = sym.Symbol('x')

    # Syntax is expression that should be 0 and variable to solve for
    roots = sym.solve(x**3 + 2*x**2 + 8, x)
    print(roots)
    pprint(roots)

If you parse through the printed items above you will see that ``sympy`` did find all 3 roots. Let's ask ``sympy`` to give us the final answer with 5 significant digits.

.. ipython:: python

    for i in range(len(roots)):
        print("Root %i = %s" % (i + 1, str(roots[i].n(5))))

Just for fun let's see if we can use ``scipy.optimize.fsolve`` to find the same roots:

.. ipython:: python

    from scipy.optimize import fsolve

    def opt_func(x):
        return x**3 + 2*x**2 + 8

    fsolve(opt_func, 0., full_output=1)
    fsolve(opt_func, -1., full_output=1)
    fsolve(opt_func, -2., full_output=1)

As you can see ``fsolve`` was able to find the correct real root, but it took a bit of work on our end to come up with a sufficiently good starting guess. It is also very important to note that we cannot solve for the imaginary roots using ``fsolve`` (or any other root finding routine in ``scipy``) because the solutions form a hyperplane.

Application
^^^^^^^^^^^

We will use ``sympy`` to help us come up with an algorithm for quadratic extrapolation. One possible case where this could be useful is taking numerical derivatives on an evenly spaced grid. By the definition of the numerical derivative you will be unable to compute the derivative at the end point. One common technique used to do this is to use some sort of extrapolation to reach beyond the grid one point and use the extrapolated data to compute the derivative at the end points.

We will start with a generic quadratic polynomial.

.. math::

    y = a + bx + cx^2

In order to solve for the 3 coefficients :math:`a`, :math:`b`, and :math:`c` we will need at least three data points: call them :math:`(x_1, y_1)`, :math:`(x_2, y_2)`, and :math:`(x_3, y_3)`. To determine the constants, you set up three equations that force the parabloa to match the data points, like this:

.. math::

    y_j = a + bx_j + cx_j^2

with :math:`j = 1, 2, 3`, and then solve for :math:`a`, :math:`b`, and :math:`c`.

.. todo::

    Assuming points :math:`x_1`, :math:`x_2`, and :math:`x_3` are equally spaced with a separation of :math:`h` (so that :math:`x_{j + 1} = x_j + h`), solve the system of three quadratic equations for the parameters :math:`a`, :math:`b`, and :math:`c` in terms of :math:`h`, :math:`x_1`, :math:`y_1`, :math:`y_2`, and :math:`y_3`. Hint: look at the sympy documentation for how to solve systems of equations.

.. todo::

    Using the coefficients :math:`a`, :math:`b`, and :math:`c` found in the previous problem, estimate the value of :math:`y(x)` at three places in terms of :math:`y_1`, :math:`y_2`, and :math:`y_3` (no :math:`x` or :math:`h` should appear in your solutions):

    1. At :math:`x _3 + h` (one point beyond the grid)
    2. Half way between :math:`x_1` and :math:`x_2`
    3. Half way between :math:`x_2` and :math:`x_3`.


========
Calculus
========

Derivatives
-----------

The derivative function in ``sympy`` can be accessed as an instance method for a ``sympy`` expression as in ``expr.diff(*vars)``. It can also be accessed as a function ``sympy.diff`` where the first argument is the expression as in ``sympy.diff(expr, *vars)``.

Notice the seemingly odd notation ``*vars``. This was done intentionally to emphasize the flexibility provided by the ``diff`` function. ``*vars`` be best understood by example. In the following list we explore some of the things ``*vars`` can be. The syntax of the list is ``expr: what derivative operation is to take place``.

- ``x``: derivative with respect to x
- ``x, x, x``: third derivative with respect to x
- ``x, 3``: third derivative with respect to x
- ``x, y, x``: derivative with respect to x, then y, then x again
- ``cos(x)``: derivative with respect to cos(x)
- ``x, 2, y, 2``: second derivative with respect to x, then second derivative with respect to y

The general rule of thumb is that any combination of expressions or symbols can be used to specify the derivative operations. There are various subtleties associated with how ``sympy`` actually carries the operations out and in order to get a more complete understanding of them; see the documentation.

To demonstrate some of these derivative operations, we will use the expressions defined above.

.. ipython:: python

    sym.pprint(ex_1)
    sym.pprint(ex_2)
    sym.pprint(ex_3)

    sym.pprint(ex_1.diff(x))
    sym.pprint(ex_1.diff(x, 2))
    sym.pprint(ex_1.diff(x, x))

    sym.pprint(ex_2.diff(x, y, z))
    sym.pprint(ex_2.diff(z, y, x))

    sym.pprint(sym.diff(ex_3, a, x, 4))


Integrals
---------

The integration functionality in ``sympy`` is also very powerful. ``sympy`` can provide definite and indefinite integrals for polynomial, rational, and exponential-polynomial. In addition, ``sympy`` can easily deal with integrals that, if they were to be solved by hand, would require multiple applications of integration by parts or u-substitutions.

Like the ``diff`` function, the ``integrate`` function can be accessed either as a method for an expression or as ``sympy.integrate``. The call signature is a bit more simple for ``sympy.integrate``, however, in that you need to give an expression and then one of two things:

- a variable of integration by itself
- a 3-element tuple (or list) ``(var, lower, upper)``, where ``lower`` and ``upper`` specify the lower and upper bounds of integration, respectively.

.. ipython:: python

    sym.integrate(x**2, x)

    # Notice the precision
    sym.integrate(x**2, (x, -2, 2))

    # Using method now
    ex_2.integrate(z)

Differential Equations
----------------------

Dynamical systems are often represented as differential equations. ``sympy`` has robust support for the following types of ordinary differential equations (ODEs):

- 1st order separable differential equations
- 1st order differential equations whose coefficients or dx and dy
  are functions homogeneous of the same order.
- 1st order exact differential equations.
- 1st order linear differential equations
- 1st order Bernoulli differential equations.
- 2nd order Liouville differential equations.
- nth order linear homogeneous differential equation with constant
  coefficients.
- nth order linear inhomogeneous differential equation with constant
  coefficients using the method of undetermined coefficients.
- nth order linear inhomogeneous differential equation with constant
  coefficients using the method of variation of parameters.

In addition to ODE support, ``sympy`` can even solve separable (either multiplicative or additive) partial differential equations (PDEs).

The main functionality for ODE solving in sympy is the ``sympy.dsolve`` function. The call signature for this function is ``sympy.dsolve(eq, goal, **kwargs)``, where ``eq`` is the differential equation, ``goal`` is the function you want to end up with, and ``**kwargs`` is a placeholder for other arguments that could be passed to the function to help it out a bit. For more detail on what ``**kwargs`` could be, see the documentation. We will demonstrate its use by solving a simple differential equation:

.. math::

    \frac{d^2f(x)}{dx^2} + 9 f(x) = 0

.. ipython:: python

    x, n = sym.symbols('x n')
    f = sym.Function('f')
    de = sym.diff(f(x), x, x) + 9*f(x)
    sym.pprint(de)

    soln = sym.dsolve(de, f(x))
    sym.pprint(soln)


Application
-----------

Suppose you are given a set of demand functions.

.. math::

    Q_d = c + b P

    Q_s = g + hP

where :math:`b < 0`  and :math:`c, g, h > 0` are constants, :math:`Q_d, Q_s` and :math:`P` are the quantity demanded, quantity supplied and price, respectively. Also suppose that price is time varying and is a positive function of excess demand:

.. math::

    \frac{dP}{dt} = m (Q_d - Q_s) \quad m > 0

.. todo::

    Use ``sympy.dsolve`` to solve the differential equation above.

.. todo::

    Analyze the solution you just found by using values ``b, c, g, h, m = (-1.2, 2, 0.1, 1.1, 1.1)``. Solve for the constant of integration in two different cases by applying two different initial conditions :math:`P(0) = .5` and :math:`P(0) = 1.3`, where t runs from 0 to 1. Plot the time series of prices. Hint: You should see that this economy is not very stable. Depending in the initial condition it will experience exponential growth or decay in prices.

Other
-----

It is worth noting that ``sympy`` has the ability to do a lot more calculus based tasks such as:

- Limits
- Summations
- Series Expansions
- Special functions (Bessel, harmonic, spherical harmonic, gamma, zeta, factorials, polynomial families)

We will not cover these topics in detail here, but the excellent `sympy documentation <sym_tutorial>`_ is a great place to start looking into these topics.

==============
Linear Algebra
==============

``sympy`` also includes a ``Matrix`` class that can be used to perform symbolic linear algebra. There are a number of ways to instantiate an object of the ``sympy.Matrix`` class and we will show a few of them below.

.. ipython:: python

    # numpy-style. Specify lists of lists, one row at a time
    sym.Matrix([[1, 2], [2, 3]])

    # Using sym.eye, sym.ones, or sym.zeros
    sym.ones(2)  # Note one argument creates square matrix
    sym.zeros((2, 4))
    sym.eye(3)

    # Using a function to specify i, j entries
    def x_mat(i, j):
        return sym.Symbol('x%i%i' % (i, j))

    sym.Matrix(4, 4, x_mat)  # give rows, cols, then function

Accessing elements of a ``Matrix`` object is more similar to MATLAB syntax than numpy in that you can specify i, j pairs for a specific element, or give a single index and access the data as if it were a flattened list. Note that like in MATLAB or numpy, slices can be passed to get chunks of matrices at at time. See the examples for clarity.

.. ipython:: python

    my_mat = sym.Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 100]])
    my_mat
    my_mat[3]

    # This should be the same as before
    my_mat[1, 0]
    my_mat[1, 1:4]

All standard matrix arithmetic works on ``Matrix`` objects. Note that unlike in numpy, the ``*`` operation does matrix multiplication, not element-wise multiplication.

There are many useful methods available to ``Matrix`` objects. Among them are:

.. hlist::
    :columns: 3

    - Some standard expression methods like ``subs``
    - ``norm``: Return a specified norm. Wraps ``numpy.linalg.norm``
    - ``det``: matrix determinant
    - ``inv``: find matrix inverse via specific algorithm`
    - ``QRdecomposition``
    - ``LUsolve``: solve linear system
    - ``GramSchmidt``: Orthogonalize columns of a matrix
    - ``LUdecomposition``
    - ``applyfunc``: apply a function element-wise

===============
Tips and Tricks
===============

Printing Formats
----------------

All expressions or symbols in ``sympy`` can be represented in various formats. Among the more useful are ``sympy.latex(ex)`` and ``sympy.ccode(ex)``, which export the expression ``ex`` directly into latex format. This can be a big time saver, especially when working with matrices.

.. ipython:: python

    print(sym.latex(my_mat))
    print(sym.ccode(ex_2))

For more information on the printing capabilities of ``sympy`` see `this page <http://docs.sympy.org/0.7.2/modules/printing.html>`_ in the documentation.

.. Make problems todo items


.. _sympy: http://docs.sympy.org/dev/install.html
.. _sym_tutorial: http://docs.sympy.org/0.7.2/tutorial.html#introduction

..    Use iPython for everything!

.. .. ipython:: python
..    :suppress:

..    import matplotlib.pyplot as plt
..    import pandas as pd
..    import numpy as np
..    np.random.seed(123456)
..    plt.close('all')
..    pd.options.display.mpl_style='default'

.. .. ipython:: python

..    from numpy.random import randn
..    ts = pd.Series(randn(1000), index=pd.date_range('1/1/2000', periods=1000))
..    ts = ts.cumsum()

..    @savefig series_plot_basic.png width=6in
..    ts.plot()