Trapezoidal rule in log-log space

Consider the problem described in this StackOverFlow post. You have a function with certain smoothness properties that are apparent on a log-log plot. This is often accompanied by a large domain of integration. It seems worthwhile to "integrate in logspace", whatever that means. 

This Jupyter notebook probes this question and makes some recommendations.

Quicktip: Reindent Python Scripts

Suppose part of a python file uses spaces for indentation, while another part uses tabs. This will throw up exceptions at runtime. So the question is how to fix it.

One answer is to use the python script reindent.py. Stick it in some folder (~/bin/) in the default path and make it executable (chmod +x reindent.py).

The usage is straightforward:

reindent -n file.py

modifies the original file in place.

Matplotlib: Saving TIFF and JPG formats

With pillow installed, on my LinuxMint installation:

import matplotlib
matplotlib.use('TkAgg') # backend

x = np.linspace(0,1)
plt.plot(x, x**2)
plt.savefig('test.tiff', dpi=300, fmt="tiff", pil_kwargs={"compression": "tiff_lzw"})

Parameter Uncertainty in Numpy Polyfit

Say you want to fit a line to (x,y) data. With polyfit, you can say,

coeff = np.polyfit(x, y, 1)

With numpy 1.7 and greater, you can also request the estimated covariance matrix,

coeff, cov = np.polyfit(x, y, 1, cov=True)

The standard error on the parameters is the square-root of the diagonal elements

print(np.sqrt(np.diag(cov)))

This report referenced in the SO page is quite useful!

QuickTip: Math Font in Matplotlib

Matplotlib (v2 and higher) uses "mathtext" to render math by default. It is quite capable, but I don't like the default font, and prefer the classic "Computer Modern" font.

You can fix this globally by modifying the rc file in your custom-style file (use the command matplotlib.get_configdir() to find location) by adding the line:

mathtext.fontset : cm

If you want to render all text using LaTeX (this slows down rendering somewhat), then use:

text.usetex : true

Reshaping NumPy Arrays into Columns

$ x = np.arange(5)
$ x
array([0, 1, 2, 3, 4])

Two different methods:

$ x.reshape(-1,1)

array([[0],

       [1],

       [2],

       [3],

       [4]])

$ np.c_[x]

array([[0],

       [1],

       [2],

       [3],

       [4]])

QuickTip: Default Color Cycle in Matplotlib

For reasonably recent versions of matplotlib [v > 1.5], you can extract the default color scheme into a string array by:

clr = [p['color'] for p in plt.rcParams['axes.prop_cycle']]

Matplotlib Tight Layouts and Plots with Insets

For normal plots or subplots, the tight_layout command does a pretty good job of keeping things from overlapping, and managing the bounding box of the overall figure.

However, if you have plot with insets etc. then tight_layout can throw a tantrum. Something like: "ValueError: max() arg is an empty sequence".

For such plots, if your axis label gets cut, and you don't want to push the label too close to the axis with something like:

ax1.xaxis.labelpad = -10

then, you can use the bbox_inches = "tight" flag when saving your figure to ensure your axis labels are not clipped.

plt.savefig('xyz.pdf',  bbox_inches = "tight")

Image Processing with Python

Basic

• matplotlib can read png and jpg files as numpy objects.
• imageio is a newer library that can read and write to a variety of image formats.
• scipy.ndimage provides some additional functionality for manipulating the image arrays.

Intermediate

• scikit-image is a library that offers a toolbox comparable to Matlab’s image processing toolbox.

Advanced

The following libraries provide more advanced functions for image manipulation.

• pillow
• openCV,
• scikit-image (to some extent)

Some resources on using Python’s basic image processing capabilities.

• Image Processing with Numpy
• SciPy Lectures on Image Manipulation
• scikit-image standard examples and introductory paper.
• SciPy2017 scikit-image video

Jupyter Notebooks: Interacting with Python Files

You can make functions defined inside a python file [myFunctions.py] visible in a Jupyter Notebook by simply importing it as a module.

For example suppose myFunctions.py contains:

$ cat myFunctions.py

def func1():

def func2():

etc.

You can use the functions by importing the python file as a module.

from myFunctions import *

This makes func1() and func2() visible inside the Jupyter notebook.

This feature can be helpful in reducing clutter by moving large walls of code out of the notebook, which can then retain a simpler look and feel.

There are two potential issues:

• when you make changes to python file, they are not immediately reflected in the notebook
• if your python file has any "script"-like commands outside the function definitions, they are executed when the file is imported as a module

The first issue can be taken care of by the magic command %autoreload.

In [1]: %load_ext autoreload

In [2]: %autoreload 2

In [3]: from myFunctions import func1

In [4]: func1()
Out[4]: 42

In [5]: # open myFunctions.py in an editor and change func1() to return 43

In [6]: func1()
Out[6]: 43

The second issue can be resolved by decorating the "script" commands with an appropriate if statement, which ensures that those commands are not executed unless the file is executed directly.

$ cat myFunctions.py

def func1():

def func(2):

#  
# Main Driver
# This part is not run when imported as a module
#
if __name__ == '__main__':
    some script commands
    print('something')

Plotting CDF: Note to Self

Consider the histogram of samples from a normal distribution:

x = np.random.normal(0., 1., size=10000)
pdf, bins = np.histogram(x, normed=True)

The size of the array "bins" is not equal to the size of "pdf". Consecutive elements of "bins" specify the left and right edges of a particular bin. Thus, by default in python, "bins" array has 11 elements, while "pdf" has 10 elements.

Note that the matplotlib command "hist" is identical in this regard.

Now suppose you want to compare the histogram with the theoretical PDF (Gaussian). Using the histogram, one could construct an equivalent line chart by taking the mid point of each bin.

# the histogram of the data
pdf, bins, patches = plt.hist(x, 30, normed=1, facecolor='green', alpha=0.4)
xpdf = (bins[1:]+bins[:-1])/2 # midpoints
plt.plot(xpdf, pdf, 'o-')

# theoretical curve
xi = np.linspace(-4, 4)
gx = 1/np.sqrt(2.*np.pi)*np.exp(-xi**2/2)
plt.plot(xi, gx, 'k--')

Everything looks fine.

Now let's consider the CDF, and plot it against the theoretical CDF. If I use bin midpoints to plot the empirical CDF I get something funky.

from scipy.special import erf

cdf  = np.cumsum(pdf)

cdf  = cdf/cdf[-1]

plt.plot(xpdf, cdf, 'o')

gcdf = 0.5*(1 + erf(xi/np.sqrt(2.)))

plt.plot(xi, gcdf)

There is a visible offset.

Instead of using bin midpoints, I should use the right limits when plotting the CDF (this makes sense upon a moments reflection!).

xcdf = bins[1:]
plt.plot(xcdf, cdf, 'o')

gcdf = 0.5*(1 + erf(xi/np.sqrt(2.)))
plt.plot(xi, gcdf)

Python Default Search Path

Suppose you have a bunch of handy python utility functions. As as example consider the matrixTeX function, to write numpy matrices in LaTeX format.

I use this function all the time, would like it to be visible, no matter where I fire my python session from. To do this we need to modify python's default search path.

Here is a three-step setup process:

1. Store the Utility Functions in a Subdirectory

There are several ways to do this.

• create a sub-directory called 'myPython' in my home directory ('/home/sachins/')
• in '/home/sachins/myPython/' create a file myPyUtils.py
• put all your utility functions in this file

> cat /home/sachins/myPython/myPyUtils.py

def matrixTeX():
    function definition

def util2():
    function definition

def util3():
    function definition

2. Point Python to the Subdirectory

You can do this temporarily or permanently.

To set it up temporarily, fire up a python session, and append the sys.path variable to expand the default search path. In our example, I would do the following:

> python3
>>> import sys
>>> sys.path.append('/home/sachins/myPython')

This setting is forgotten when you terminate your python session. Thus, you have to invoke this pair of commands, every time you start a session.

To avoid this, you can set things up so that they are more permanent.

Open your .bashrc file (if you are using bash) and modify the PYTHONPATH variable.

In this case, I  add the following lines:

# set python path to look for scripts

export PYTHONPATH=$PYTHONPATH:/home/sachins/myPython

Now you are all set. We can import the utility functions and use them from anywhere.

3. Import the Utilities Functions

Fire up a python session. Once you have expanded the search path (temporarily or permanently), you can import the utility file (myPyUtils.py) as a module, and access its individual functions. As an example,

>>> import numpy as np

>>> from myPyUtils import *

>>> print(matrixTeX(np.array([1,2])))

\begin{bmatrix}

  1 & 2\\

\end{bmatrix}

Notes:

• the first time you import the file, it will create a ".pyc" file in the subdirectory
• you can be more careful with the namespace during importing

Links to matplotlib Resources

I wanted to pull together a list of matplotlib resources that I need to consult frequently.

1. SciPy Lectures: The entire series is great, including the introduction to matplotlib.

2. Tutorials from J.R. Johansson and Nicholas P. Rougier

3. A couple of my own Jupyter notebooks on customizing styles, and multiplots.

Google Colaboratory

If you need to use and interact with a jupyter notebook on a computer that does not have it installed, Google Colaboratory seems like a great in-browser solution. I learned about it from a student earlier this semester from a student.

The best part is that you don't need to install any software locally on the machine. The standard scientific/data science python stack (numpy, scipy, sympy, pandas) is available, and you can even "install" some additional on the fly using pip install.

It works more or less like Google Docs, in that you documents are saved on Google Drive, and you can collaborate with others in much the same way.

Check it out!

Exporting Numpy Arrays and Matrices to LaTeX

Over the past couple of years, a lot of my "numerical experimentation" work has moved from Octave to python/numpy.

I incorporate a lot of this work into my classes and presentations (made using beamer), and having a script to translate vectors and matrices to LaTeX format is handy.

In the past, I shared a Matlab/Octave script which does this.

Here is a python/numpy script which does something similar. The script

• autodetects integers and floats
• allows you to control the number of decimals for floats
• allows you optionally render floats in scientific format
• right-justify using the bmatrix* environment (good for -ve numbers)
• suppress small values near zero (~ 1e-16)

Python: Orthogonal Polynomials and Generalized Gauss Quadrature

A new (to me) python library for easily computing families of orthogonal polynomials.

Getting standard (generalized) Gauss quadrature schemes is extremely simple. For example to get 13 nodes and weights for Gauss-Laguerre integration, correct up to 50 decimal places:

pts,wts = orthopy.schemes.laguerre(13, decimal_places=50)

The numpy Polynomial package provides similar functionality:

pts, wts = numpy.polynomial.laguerre.laggauss(13)

A nice feature (besides arbitrary precision) is that you can derive custom orthogonal polynomials and quadrature rules. All you need to provide is a weight function and domain of the polynomials. From the project webpage:

import orthopy
moments = orthopy.compute_moments(lambda x: x**2, -1, +1, 20)
alpha, beta = orthopy.chebyshev(moments)
points, weights = orthopy.schemes.custom(alpha, beta, decimal_places=30)

This generates a 10-point scheme for integrating functions over the interval [-1, 1], with weight function \(w(x) = x^2\).

Matplotlib: Subplots, Inset Plots, and Twin Y-axes

This jupyter notebook highlights ways in which matplotlib gives you control over the layout of your charts. This is intended as a personal cheatsheet.

Matplotlib Styles

I created a jupyter notebook demonstrating the use of built-in or customized styles in matplotlib, mostly as a bookmark for myself.

PyCon 2017 Talks

Some interesting Python talks (links to YouTube videos) from this year's PyCon.

1. Jake Vanderplas: The Python Visualization Landscape

2. Chistopher Fonnesbeck: PyMC3

3. Eric Ma: Bayesian analysis

4. Alex Orlov: Cython

5. Bret Cannon: What new is python 3.6?

QuickTip: Sorting Pairs of Numpy Arrays

Consider the two "connected" numpy arrays:

import numpy as np

x = np.array([1992,1991,1993])

y = np.array([15, 20, 30])

order = x.argsort()

x     = x[order]

y     = y[order]

x = array([1991, 1992, 1993])

y = array([20, 15, 30])