DEV Community: Justin Swaney

Computing pretrained image embeddings in Python

Justin Swaney — Mon, 29 Aug 2022 13:17:00 +0000

Photo by Matthew Ansley on Unsplash

In this tutorial, we'll create a Python package for computing image embeddings using a pretrained convolutional neural network (CNN) from the command line.

An image embedding is a numerical representation of an image that "embeds" it within a high-dimensional vector space. Image embeddings are useful for a variety of applications, including content-based information retreival and image classifiation. Furthermore, there are many off-the-shelf pretrained CNN architectures that you can use whenever you need. Machine learning engineers often use these pretrained models for benchmarks or as useful starting points.

After completing this tutorial, you will know how to:

Load image data using numpy and scikit-image
Compute linear image embeddings using PCA from scikit-learn
Compute deep image embeddings using a pretrained ResNet model from Keras
Remove any number of layers from a pretrained model in Keras
Visualize deep image embeddings using UMAP

We are going to make a command line tool called embedding to do all of this. As a final demonstration, we will compare PCA with ResNet image embeddings of MNIST handwritten digits.

Setting up the `embedding` package

We're going to skip scaffolding out the embedding Python package and only walkthrough the development. First, you'll need to clone this repository:

git clone https://github.com/jmswaney/blog.git
cd blog/07_pretrained_image_embeddings

Within this folder, you'll see a pyproject.toml file as well as an environment.yml file for poetry and conda, respectively. You'll need to have conda installed on your system, and you can either follow my previous tutorial to install poetry globally or just use the one included in the embedding virtual environment. To install the virtual environment, you can run:

conda env create -f environment.yml
conda activate embedding

Note: I've included cudatoolkit and cudnn in the environment.yml dependencies for GPU support. You can remove these dependencies if you don't want GPU acceleration.

Once installed and activated, we just need to install the Python dependencies, including the embedding CLI itself:

# from 07_pretrained_image_embeddings/
poetry install

This will install our machine learning dependencies like scikit-learn and tensorflow. It will also create an executable console script with its entrypoint connected to embedding:cli.

Developing the `embedding` CLI

The embedding folder has the completed source code for this tutorial, so if you want to follow along with the development steps yourself, you may want to rename it and create an empty folder for your work.

Creating a simple entrypoint

The first thing we need to implement is a simple entrypoint function called cli and expose it from the package __init__.py. We can create cli.py and add a simple click command like in my CLI tutorial.

# in cli.py
import click

@click.command()
def cli():
    print("Hello, embedding!")

To expose this as embedding:cli, we just need to import it within our package's __init__.py.

# in __init__.py
__version__ = '0.1.0'

from .cli import *

Now we can check that our entrypoint is working from the terminal,

embedding
Hello, embedding!

By default, poetry installs the current package in "editable" mode, which means our source changes will take effect immediately in our CLI.

Loading image data

Next, we need to be able to load an image from a file. The best way I've found to do this is to use the imread function from scikit-image, which offers great speed and support for many lossless and lossy image formats. We can update the cli function to open an input image file.

import click
from skimage.io import imread

@click.command()
@click.argument("input")
def cli(input):
    print(main(input))

def main(input):
    return imread(input).shape

Testing this with a picture of a dog, it prints the shape of the image.

embedding data/dog.jpg
(333, 516, 3)

Now, we may also want to load data from Numpy arrays, so let's factor this out into a module specifically for I/O operations called io.py.

# in io.py
import os
import numpy as np
from skimage.io import imread

NUMPY_EXT = ".npy"

def load(filename: str) -> np.ndarray:
    ext = os.path.splitext(filename)[-1]
    if ext == NUMPY_EXT: return np.load(filename)
    return imread(filename)

Then loading .npy and .npz files should work as well,

embedding data/dog.npz
(333, 516, 3)

embedding data/dog.npy
(333, 516, 3)

Linear embeddings with PCA

Before we get to using pretrained image models, let's set the stage with principal component analysis (PCA). PCA is an unsupervised learning algorithm that will calculate a linear projection to a lower-dimensional space that preserves as much of the variance in the original dataset as possible.

We'll need multiple images for this, so let's randomly sample square patches from our example image. Then, in a new embed.py module, we can use sklearn to calculate our PCA embedding.

# in cli.py
from embedding.io import load
from embedding.embed import pca
import numpy as np
from random import random
from math import floor

def main(input):
    data = load(input)

    size = 64
    start = lambda: (floor(random()*(data.shape[0] - size)), floor(random()*(data.shape[1] - size)))
    x = np.asarray([data[y:y+size, x:x+size] for (y, x) in [start() for _ in range(32)]])
    print("Patches shape", x.shape)

    x_pca = pca(x.reshape(len(x), -1), 4)
    print("PCA embedded shape", x_pca.shape)

Since PCA applies to vector inputs, we needed to reshape our images into vectors. When we run this, we can see the data dimensionality is reduced, as expected.

embedding data/dog.jpg

Patches shape (32, 64, 64, 3)
PCA embedded shape (32, 4)

Here, we've sampled 32 image patches, each (64, 64, 3), and we have chosen to use 4 principal components.

Deep embeddings using ResNet50V2

Now that we have our PCA image embeddings working, we can move to a pretrained neural network. Keras has several CNNs available with weights obtained from ImageNet pretraining here. To keep things simple, we'll just use the ResNet50V2 achitecture for demonstration.

Within our embed.py module, we can import the model and instantiate it with the ImageNet weights. The first time we do this, Keras will download the weights and keep them locally cached. We need to make sure that we instantiate the model with include_top=False so that we get the last hidden representation in the network (as opposed to class label probabilities from the final softmax layer).

import numpy as np
from sklearn.decomposition import PCA
import tensorflow as tf
from tensorflow.keras.applications.resnet_v2 import preprocess_input

def pca(x: np.ndarray, n: int) -> np.ndarray:
    return PCA(n_components=n).fit_transform(x)

def resnet(x: np.ndarray) -> np.ndarray:
    model = tf.keras.applications.ResNet50V2(include_top=False)
    maps = model.predict(preprocess_input(x))
    if np.prod(maps.shape) == maps.shape[-1] * len(x):
        return np.squeeze(maps)
    else:
        return maps.mean(axis=1).mean(axis=1)

Notice that we also used the Keras utility function preprocess_input for the resnet_v2 family of models. This is important to make sure that we are scaling our images to intensity values between -1 and 1, as the model was trained on images that were scaled similarly. We also perform global average pooling to return a vector representation, effectively removing any remaining spatial information left in the activation maps.

When we call resnet with our image patches in cli.py, we get the following,

embedding data/dog.jpg

Patches shape (32, 64, 64, 3)
PCA embedded shape (32, 4)
ResNet embedded shape (32, 2048)

This shows that the ResNet50V2 architecture provides a 2048-dimensional embedding from the last layer of the model.

Note: if you encounter GPU memory allocation issues, then you may need to allow memory growth using this suggestion.

Saving image embeddings

Now that we can compute our image embeddings, we need a way for our CLI to save them to a file. Back in our io.py module, we can use numpy again to save our image embeddings with optional lossless compression.

# in io.py
def save(filename: str, x: np.ndarray, compress: bool) -> None:
    np.savez_compressed(filename, x) if compress else np.save(filename, x)

At this point, it would be helpful to refactor our main function into a slightly more scoped embed function that uses a given model to embed images as numpy arrays.

# in embed.py
def embed(x: np.ndarray, model: str, n: int):
    match model:
        case "pca":
            return pca(x.reshape(len(x), -1), n)
        case "resnet":
            if x.ndim == 3:
                return resnet(np.repeat(x[..., np.newaxis], 3, axis=-1))
            return resnet(x)
        case _:
            raise Exception("Invalid model provided")

We've repeated black and white images in the last dimension to match the expected input dimensionality for ResNet50V2.

We can use this in cli.py which uses io.py to perform side-effecting operations.

# in cli.py
MODEL_TYPES = ['pca', 'resnet']

@click.command()
@click.argument("input")
@click.argument("output")
@click.option("--model", "-m", required=True, type=click.Choice(MODEL_TYPES))
@click.option("--dimensions", "-d", "n_dim", default=32, type=int)
@click.option("--compress", "-c", is_flag=True)
def cli(input, output, model, n_dim, compress):
    x = sample_patches(load(input), 32, 64)
    x_embed = embed(x, model, n_dim)
    save(output, x_embed, compress)

We've added arguments to write the output file with optional compression as well as choose which model to use and the embedding dimension (for PCA). We've also refactored the image patch sampling code into the sample_patches function.

Loading MNIST data

Since it's uncommon to be embedding a single image, we can remove our sample_patches function. If the input is a numpy array, we assume it has the proper dimensionality for the downstream models. If the input is a folder, we should assume it is a flat folder of only images which should all be embedded in alphanumerical order. We'd also like the ability to load the MNIST training images. To perform all these loading operations, we can edit the load function in io.py,

from typing import Optional
import os
import numpy as np
from skimage.io import imread
from tensorflow.keras.datasets.mnist import load_data

MNIST_NAME = "mnist"
MNIST_CACHE = "mnist.npz"
NUMPY_EXT = ".npy"

def load(name: str) -> tuple[np.ndarray, Optional[np.ndarray]]:
    if name == MNIST_NAME: 
        return load_data(MNIST_CACHE)[0] # Only training set
    ext = os.path.splitext(name)[-1]
    if ext == NUMPY_EXT: return np.load(name)
    return np.asarray([imread(os.path.join(name, x)) for x in os.listdir(name)])

So if we pass "mnist" as the first argument to embedding, we'll get the image data and labels from the MNIST dataset in Keras. Otherwise, we'll load the given matrix or folder of images.

Generating UMAP visualizations

We often want to visualize high-dimensional embeddings on a 2D plot. The most common tools for this are TSNE and UMAP, so for our example, we'll add an option to generate a UMAP plot of the image embeddings.

Back in our embed.py module, we'll add the following:

from umap import UMAP

def umap(x, seed = None):
    if seed: np.seed(seed)
    return UMAP().fit_transform(x)

We can call this in cli.py to save our results to a file.

from embedding.io import load, save
from embedding.embed import embed, umap
from embedding.plot import plot_umap, plot_pca
import os
import click

MODEL_TYPES = ['pca', 'resnet']

@click.command()
@click.argument("input")
@click.argument("output")
@click.option("--model", "-m", required=True, type=click.Choice(MODEL_TYPES))
@click.option("--dimensions", "-d", "n_dim", default=32, type=int)
@click.option("--compress", "-c", is_flag=True)
@click.option("--plot", "-p", is_flag=True)
@click.option("--seed", "-s", type=int, default=None)
def cli(input, output, model, n_dim, compress, plot, seed):
    x, y = load(input)
    print("Input shape", x.shape)
    x_embed = embed(x, model, n_dim)
    print("Embedded shape", x_embed.shape)
    save(output, x_embed, compress)
    if plot:
        root, ext = os.path.splitext(output)
        name = os.path.basename(root)
        if x_embed.shape[-1] == 2:
            plot_pca(x_embed, y, f"{name}_pca.png")
        else:
            x_umap = umap(x_embed, seed)
            print("UMAP shape", x_umap.shape)
            save(f"{name}_umap{ext}", x_umap, compress)
            plot_umap(x_umap, y, f"{name}_umap.png")

We added a _umap suffix to the output plot and matrix filenames for UMAP results (if the embedding dimension isn't already 2). If the embedding dimension is 2, then we just assume it's from PCA, and we plot it without performing UMAP.

We can write the plot_pca and plot_umap function in plot.py.

from typing import Optional
from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np

sns.set(style="white", context="poster", rc={ "figure.figsize": (14, 10), "lines.markersize": 1 })

def make_scatter(x: np.ndarray, y: Optional[np.ndarray]):
    if y is None:
        plt.scatter(x[idx, 0], x[idx, 1])
    else:
        for label in np.unique(y):
            idx = np.where(y == label)[0]
            plt.scatter(x[idx, 0], x[idx, 1], label=label)
            plt.legend(bbox_to_anchor=(0, 1), loc="upper left", markerscale=6)

def plot_umap(x: np.ndarray, y: Optional[np.ndarray], filename: str):
    make_scatter(x, y)
    plt.xlabel("UMAP 1")
    plt.ylabel("UMAP 2")
    plt.title(filename)
    plt.savefig(filename)

def plot_pca(x: np.ndarray, y: Optional[np.ndarray], filename: str):
    make_scatter(x, y)
    plt.xlabel("PC 1")
    plt.ylabel("PC 2")
    plt.title(filename)
    plt.savefig(filename)

In make_scatter, we use seaborn and matplotlib to scatter plot all 2D image representations along with class labels if we have them.

Removing any number of layers

It's often helpful to remove more layers than just the top softmax layer, and this is an area that you can experiment with. Since convolutional architectures are typically structured in blocks, it's often worthwhile to remove entire blocks when trying to use pretrained models in a new image domain.

We can add a command-line argument to specify how many layers we'd like to keep (at which point in the original model do we stop and use as our image embedding). If this argument is not specified, then we take the enitre model without the top.

# Add option to cli.py
...
@click.option("--layers", "-l", type=int, default=None)
def cli(input, output, model, n_dim, layers, compress, plot, seed, batch_size):
    ...
    x_embed = embed(x, model, n_dim, batch_size, layers)
    ...

# Remove layers in embed.py
def resnet(x: np.ndarray, batch_size = None, layers = None) -> np.ndarray:
    full = tf.keras.applications.ResNet50V2(include_top=False)
    model = full if layers is None else tf.keras.models.Model(full.input, full.layers[layers].output)
    ...

Now if we pass --layers 100, we'll get the first 100 layers of the ResNet50V2 model.

Results

Now that we've completed the initial development of our embedding tool, it's time to use it. First, let's look at a simple 2D PCA image embedding.

# Embed MNIST using 2D PCA
embedding mnist mnist_pca_2d.npy --model pca -d 2 -p

The first principal component looks like it corresponds roughly to having a hole in the center of the digit (zero-like).

Now let's compare this to our ResNet50V2 image embedding visualized with UMAP.

# Embed MNIST using ResNet50V2
embedding mnist mnist_resnet.npy --model resnet -p -s 582

It seems that in both image embeddings, 0 and 1 are the most obvious digits to identify. While the ResNet image embeddings have 0 and 1 separated slightly better than the PCA plot, it's surprisingly not that much better. Perhaps we can use fewer layers on this simpler dataset.

After some experimentation, I've found that removing blocks 4 and 5 from the ResNet50V2 model give image embeddings that are closer to linearly separable. Removing blocks 4 and 5 corresponds to taking the first 82 layers of the model.

# Embed MNIST using 82 layers of ResNet50V2
embedding mnist mnist_resnet_82layers.npy --model resnet -p -s 582 -l 82

Final thoughts

Making use of pretrained models is an essential tool in the toolbox for machine learning engineers. Pretrained image embeddings are often used as a benchmark or starting point when approaching new computer vision problems. Several architectures that have performed extremely well at ImageNet classification are avilable off-the-shelf with pretrained weights from Keras, Pytorch, Hugging Face, etc.

Here, we've created a simple Python package to help in the computation of image embeddings from the command-line. We used our embedding CLI to compute and visualize MNIST image embeddings using PCA and a ResNet50V2 model. In this relatively simple dataset, we found that using a subset of layers from the ResNet50V2 model performed qualitatively better than PCA at distinguishing between handwitten digits (without supervised training).

We might have expected the earlier layers of the ResNet50V2 model to be more applicable to the MNIST dataset because earlier layers in CNN's are known to extract features related to edges and textures. The higher-level visual features needed to classify the thousands of RGB natural images in ImageNet may not be as helpful when dealing with black and white images of handwritten digits.

Source availability

All source materials for this article are available here on my blog GitHub repo.

Tips for using Jupyter Notebooks with GitHub

Justin Swaney — Mon, 22 Aug 2022 15:07:00 +0000

Photo by Snowscat on Unsplash

To machine learning engineers and data scientists, Jupyter notebooks have become a workhorse tool in the toolbox. Notebooks allow researchers to quickly prototype new solutions and make interactive visualizations in a remote computing environment. Unfortunately, the ease-of-use and interactivity of notebooks comes with some trade-offs when using version control.

Version control tools like git are powerful collaboration tools that track changes to source code and synchronize local and remote copies of a shared codebase. They allow developers to work together on the same codebase and seamlessly merge their improvements together. These tools also work out of the box with services like GitHub, which provide hosting for the shared codebase. Unfortunately, most of the user experience of using git and GitHub has been designed around text-based source code, but Jupyter notebooks are saved with embedded output media in a JSON format within .ipynb files.

As an example, let's imagine that we are quants at a new hedge fund called "Patagonia Capital". We decide it would be a good idea to create a repository to hold all of our exploratory data analysis (EDA) which will include shared Jupyter notebooks.

After completing this guide, you will learn how to:

Setup auto-formatting in Jupyter notebooks
Use VSCode's built-in tools for viewing notebook diffs
Use nbstripout with git to only commit changes to notebook source
Use papermill to create a shared library of executable notebooks
And fetch some stock market data in Python!

After creating the "patagonia" Github repository for all the quants at Patagonia Capital, our next task is to setup the patagonia environment.

Tip #1: Auto-formatting in Jupyter notebooks

In our newly created repository, we can create the patagonia conda virtual environment. It's helpful to setup code formatting to provide some reasonable formatting standards across our company.

conda create -y -n patagonia python=3.8 
conda activate patagonia
conda install -y -c conda-forge pylint black jupyterlab jupyterlab_code_formatter

We can start jupyter lab to check that the code formatting is working. When you create a new notebook, you should see a Format notebook option in the notebook's toolbar.

Since it's somewhat inconvenient to have to click the toolbar to perform code formatting, we can configure our formatter to format on save. Open up the Advanced Settings Editor via the Settings dropdown (or by pressing Cmd-, on Mac). Here, you will see the "Jupyterlab Code Formatter" settings for customization. At the bottom of these options, you can turn on the Auto format config checkbox.

Note: if you are using an older version of Jupyter, you may not see a checkbox for this setting. In this case, you can add the following under "User Preferences":
{
  "formatOnSave": true
}

When you return to your notebook, you should see that simply saving the notebook automatically formats your code. This will help make our changes easier to read within version control because our code will be committed with consistent formatting that is also designed to make diffs as small as they can be.

Tip #2: VSCode notebook diffs

Now that we have our auto-formatting setup, let's plot the closing price of NVDA stock over the last 5 years. We'll first need to install some new dependencies.

conda install -y -c conda-forge pandas yfinance seaborn

Next, we can use Yahoo Finance along with Seaborn to make our plot of the NVDA ticker history in a Jupyter notebook:

# in plot_history.ipynb
import seaborn as sns
import yfinance as yf
from matplotlib import pyplot as plt

sns.set_theme("paper")

ticker = yf.Ticker("NVDA")
hist = ticker.history(period="5y")

plt.figure()
sns.lineplot(x="Date", y="Close", data=hist)
plt.savefig("NVDA_close_5y.png")
plt.show()

This will save the plot as a PNG image embedded in the notebook as well as in a new file called NVDA_close_5y.png. Since we're just getting started, let's imagine we directly commit this notebook and push it into the patagonia repository.

When we try to change something simple in our notebook like the ticker symbol from NVDA to AAPL, this is the diff we get in GitHub.

On the top, you can see the plot image changed (presumably both base-64 encoded). On the bottom, you can see the only change to the code, but it is inconveniently embedded in a string with escape characters.

If you were to look at the diff in VSCode, however, you would be pleasantly surprised to see this:

In VSCode, comparing different versions of Jupyter notebooks will show more informative diffs with your source code as well as changes in the rendered output.

Tip #3: nbstripout as git filter

Using auto-formatting and VSCode notebook diffs are really helping the devs at Patagonia Capital work togther, but sometimes we'd like to be a little more strict about what gets committed to our repository.

At Patagonia Capital, we've had some issues with new developers pushing massive notebooks to GitHub. To solve this type of problem, we can use nbstripout to automatically remove all notebook outputs before committing to our repository. To setup nbstripout, we need to do the following within our patagonia repository:

conda install -y -c conda-forge nbstripout
nbstripout --install --attributes .gitattributes

This will add a .gitattributes file (similar to a .gitignore file) which will tell git to apply the nbstripout filter to all .ipynb files. When this filter is applied to our notebooks, GitHub will no longer have the output plots.

If you'd like to automatically remove empty / tagged cells or retroactively apply this filter to your git history, you can read the nbstripout documentation on GitHub.

Tip #4: Executable notebooks with Papermill

While nbstripout has really helped Patagonia Capital avoid a bloated repository, we still sometimes would like to generate reports with standard visualizations from our growing library of notebooks.

Papermill allows running a Jupyter notebook as if it were an executable script. By designating certain cells as "parameters", papermill provides a command-line interface (CLI) to execute your notebooks and to specify parameter values via command-line arguments.

Let's install papermill and create an output/ folder for our "rendered" notebooks.

conda install -y -c conda-forge papermill
mkdir output

To avoid committing redundant rendered notebooks, we can add output/ to our .gitignore.

In Jupyterlab, we can bring our symbol, stat, and period constants to a separate cell. With that new cell active, we can open the Property Inspector and add a tag called "parameters":

The values in this tagged cell will be default values to the notebook. We can execute our notebook to create a rendered notebook for the last 10 years of AAPL using the papermill CLI:

papermill plot_history.ipynb \
    output/plot_history_rendered.ipynb \
    -p symbol AAPL \
    -p period 10y

Papermill can also target cloud storage outputs for hosting rendered notebooks, execute notebooks from custom Python code, and even be used within distributed data pipelines like Dagster (see
Dagstermill). For more information, see the papermill documentation.

Final thoughts

Jupyter notebooks and git are powerful tools for machine learning engineers and data scientists to prototype solutions and collaborate on a shared codebase. We discussed some tips and tricks for getting the best of both worlds from notebooks and version control in our simple financial example from Patagonia Capital.

Source availability

All source materials for this article are available here on my blog GitHub repo.

Setting up Python Linting and Auto-formatting in VSCode

Justin Swaney — Mon, 07 Feb 2022 02:31:23 +0000

In this tutorial, we will setup Python linting and automatic code formatting in VSCode. We will use pylint to highlight linting errors and black for auto-formatting our code on save.

Linting and auto-formatting are commonly used in Javascript projects and help to enforce consistent code standards across a team of developers working together in a single codebase. For anyone using Python on a team, pylint and black could help speed up development, reduce time spent upskilling new hires, and make your diffs easier to read for code reviews.

Setup

First, let's start by making a new virtual environment for this tutorial called python-linting.

conda create -y -n python-linting python=3.8
conda activate python-linting

Next, we'll make sure that we have our linter and formatter installed.

conda install -y pylint black

Now that we have our tools installed, we can move on to setting them up in VSCode.

Linting with `pylint`

To enable linting for Python in VSCode, we need to make sure some things are enabled in our settings. Be sure to make these changes to the "User" level (as opposed to the Workspace or folder only) so that it takes effect in any project.

First, open your VSCode Settings UI and search for "Python linting". Make sure that the following setting is checked:

Python > Linting: Enabled
[x] Whether to lint Python files.

Also specifically allow pylint linting by making sure the following setting is checked:

Python > Linting: Pylint Enabled
[x] Whether to lint Python files using pylint.

You can optionally suppress "convention"-level linting messages if you want to only highlight more egregious issues. To disable convention messages, add --disable=C to the pylint arguments:

Python > Linting: Pylint Args
--disable=C

Only do this part if you know that you don't want convention linting errors shown.

Auto-formatting with `black`

To setup auto-formatting, we need to tell VSCode to use black as our formatter for Python code. We can do this in our settings by selecting black from the following dropdown selector:

Python > Formatting: Provider
black

Next, we need to tell VSCode that we'd like to auto-format our code on save of the file. We can do this with the Format on Save option:

Editor: Format on Save
[x] Format a file on save. A formatter must be available, the file must not be saved after delay, and the editor must not be shutting down.

Result

If everything went well, you should get helpful linting messages while you code and automatic code formatting each time you save your file.

Final thoughts

In this tutorial, we setup pylint and black for Python linting and auto-formatting in VSCode. Linting and auto-formatting can help speed up development by removing time spent debugging and thinking about formatting. They can also help a team of developers work together by enforcing a reasonable set of standards that strive to make their codebase and pull requests easier to read and review.

Source availability

All source materials for this article are available here on my blog GitHub repo.

Poetry: a yarn-like package manager for Python

Justin Swaney — Fri, 04 Feb 2022 06:18:28 +0000

In this tutorial, we will create a Python package called math-demo using the poetry Python package manager. We will use poetry to scaffold out the math-demo package, install our dependencies in a conda virtual environment, and run scripts for unit testing and making plots.

After completing this tutorial, you will learn how to:

Install the poetry package manager
Setup a new Python package with poetry
Use poetry with conda virtual environments
Run unit tests and other scripts with poetry
And some math!

Installing the `poetry` package manager

The installation instructions for poetry are available here. To install poetry on macOS, we run:

curl -sSL https://raw.githubusercontent.com/python-poetry/poetry/master/get-poetry.py | python -

We can also install poetry tab completions to help us remember all subcommands and options from the command line. The instructions for most popular shells are here. If you are using Oh-My-Zsh, then you'd run:

# Only for Oh-My-Zsh
mkdir $ZSH_CUSTOM/plugins/poetry
poetry completions zsh > $ZSH_CUSTOM/plugins/poetry/_poetry

Once this is done, we need to add poetry to our list of plugins in our ~/.zshrc file.

After restarting the terminal, running poetry --version prints our version of poetry. We now have tab completion listing out possible poetry commands as well. It's useful to note that poetry is installed globally similar to tools like conda and yarn.

Creating our `math-demo` package

Now that we have poetry installed globally, we can use it to scaffold out a new Python package for us:

poetry new math-demo
cd math-demo

This will build out the basic structure of a Python package in poetry's standard format in the math-demo folder. We should find a subfolder named math_demo which, similar to other Python packages, is where our source code will go. There should also be a tests folder for unit tests and a curious new pyproject.toml file at the top-level.

Historically, Python packages would have a setup.py file and optionally a requirements.txt file which can be used to install the Python package locally. Instead of writing the code to install the Python package, our pyproject.toml file is just a file with data in it that declares information about the Python package in a standard way.

If you are familiar with the package.json in web-development, the pyproject.toml file is very similar. A package.json specifies package dependencies, development dependencies, and scripts for running tests and creating minified builds. Some rough analogies between Python and Javascript would be:

poetry = yarn / npm
pyproject.toml = package.json
poetry.lock = yarn.lock / package-lock.json

If we open up the pyproject.toml file, we can see that package dependencies and development-only dependencies go under [tool.poetry.dependencies] and [tool.poetry.dev-dependencies], respectively.

But what about `conda`?

If you are familiar with conda, then you might think the pyproject.toml file is equivalent to an environment.yaml file. The environment.yaml file was created specifically for recreating consistent runtime environments across machines, whereas the pyproject.toml file was first introduced in PEP 518 as a way to move Python builds to a community standard. Both of these tools have merits, and they are not mutually exclusive.

Here, we will use conda to handle our virtual environment as we normally would, and then poetry will be our main package manager for dependences within that environment. This has the added benefit of having your IDE automatically recognize your new virutal environment if you already have conda setup.

Let's start by adding a bare-bones environment.yaml next to our pyproject.toml.

# environment.yaml
name: math-demo

channels:
  - default
  - conda-forge

dependencies:
  - python=3.8

Running conda env create -f environment.yaml will then create the math-demo virtual environment. We can then activate it in your IDE or by running conda activate math-demo.

Now that we have our virtual environment setup, we can use poetry to add our dependencies to it. As mentioned in the documentation, poetry can manage environments itself, but it will also take your lead as to which environment it should use. That's why this works with the active conda environment.

Installing our dependencies

First of all, let's double-check that we are in the correct environment before installing anything:

poetry env info
# => Virtualenv
# => Path:  /Users/{username}/anaconda3/envs/math-demo

From the path, we can see that poetry is using the correct virtual environment. Let's add plotly, numpy, and pandas so we can do some math!

poetry add plotly numpy pandas

Note: If you are using VSCode on macOS and do not want to use conda for virtual environments, you can go into your settings and add Library/Caches/pypoetry/virtualenvs to your Python: Venv Folders. This will allow you to select the correct Python interpreter and, importantly, let VSCode find the packages you have just installed for linting.

Catenary curves and their parabolic approximations

Think of a power line hanging between two telephone poles. What shape does the hanging wire make? It kinda looks like a really wide parabola. It turns out that shape is called a catenary. It can be derived using calculus of variations after the simple observation that the curve minimizes the total gravitational potential energy along its arc length. That is, if we were to make a small variation in the way the wire was hanging, it would settle back down to the catenary eventually. You can try this out with string; don't touch power lines.

The formula for a catenary curve centered at x = 0 is given by:

where a is a constant.

Since we are curious how closely a parabola approximates this curve, we can write the second-order Maclaurin series for it.

Aren't you excited to plot these equations?! Or are the power lines looking pretty good right now?

Implementing the exact and approximate catenary calculations

Next, we are going to add code to calculate f(x) and g(x) to our math-demo package.

Writing a failing test

First, we can write a test knowing that the curve should be symmetric about the y-axis (since it is centered at x = 0). To do this, start by making a new Python module called hyperbolic.py inside math_demo next to __init__.py.

Let's write the function we want to test and start with something that will definitely fail our symmetry test:

# in hyperbolic.py
def catenary(x):
    return x

Great, now let's change our test suite to use this:

# in test_math_demo.py
from math_demo import hyperbolic
import numpy as np

N_POINTS = 10

def test_catenary_symmetric():
    x = np.asarray([2 ** i for i in range(N_POINTS)])
    assert np.all(hyperbolic.catenary(x) == hyperbolic.catenary(-x))

Here, we are asserting that f(x) = f(-x) for 10 choices of x. Running poetry run pytest should show this test fails as expected.

Making the tests pass

Now we can work on making our test pass. Adding the exact catenary formula, we get:

# in hyperbolic.py
import numpy as np

def catenary(x, a = 1):
    return a * np.cosh(x / a)

Now running poetry run pytest shows that our symmetry test passes!

Adding a parabolic approximation

Below the catenary function, let's add our parabolic approximation in a new function called catenary_approx.

# in hyperbolic.py
def catenary_approx(x, a = 1):
    return a + x ** 2 / (2 * a)

Now we are ready to see how well our approximation does for various values of a.

Visualizing exact and approximate catenary curves

We can use plotly to visualize the exact and approximate catenary curves. We can add a main function to hyperbolic.py and only run it if the module is executed directly (as opposed to imported).

# in hyperbolic.py
def main():
    import plotly.express as px
    import pandas as pd

    x = np.linspace(-2, 2)
    a = [1, 1.4, 1.8]

    df = pd.DataFrame()

if __name__ == "__main__":
    main()

Next, we need to build a long-form DataFrame with all (x, y) coordinates, line group, color, and line style labels. We will color curves by the value of the a parameter and make our approximations dashed lines.

# in hyperbolic.py
def main():
    from itertools import cycle
    import plotly.express as px
    import pandas as pd

    x = np.linspace(-2, 2)
    a = [1, 1.4, 1.8]

    df = pd.DataFrame()
    cols = [f"${func}_{{{ai:0.1f}}}$" for func in ["f", "g"] for ai in a]
    for col, ai in zip(cols, cycle(a)):
        df2 = pd.DataFrame({"x": x})
        df2["y"] = caternary(x, ai) if "f" in col else caternary_approx(x, ai)
        df2["line"] = col
        df2["a"] = ai
        df2["form"] = "Exact" if "f" in col else "Approx"
        df = pd.concat([df, df2])

    fig = px.line(
        df,
        x="x",
        y="y",
        line_group="line",
        color="a",
        line_dash="form",
        title="Catenary Demo",
    )
    fig.update_yaxes(range=[0.5, 4])
    fig.show()

if __name__ == "__main__":
    main()

To see the results, we can run:

poetry run python math_demo/hyperbolic.py

This will open the plot in your browser to show our catenary curves and decent parabolic approximations to them.

Final thoughts

In this tutorial, we learned how to install and use the poetry package manager to make a new Python package. We separated our pytest development dependencies from our package dependencies using the pyproject.toml file. We were also still able to still use a conda virtual environment with poetry.

We also did a little bit of math about catenary curves! If you want to see the details of how the exact catenary formula can be derived, here is a post from the University of Virginia.

Source Availability

All source materials for this article are available here on my blog GitHub repo.

Creating beautiful command-line interfaces in Python

Justin Swaney — Thu, 04 Mar 2021 13:45:26 +0000

Photo by Ryan Quintal on Unsplash

In this tutorial, we will create a Python package called starfox with a command-line interface (CLI) that simulates conversations between players from Star Fox 64. We will create the CLI using the click package and add a wizard using questionary for conveniently choosing which characters to include in our simulated conversation.

After completing this tutorial, you will learn how to:

Setup a Python package with custom entrypoints
Use editable mode to see the effects of our changes as we develop
Define a simple CLI with built-in help messages
Build a command-line wizard to easily run CLI commands

Setting up the `starfox` Python package

First, create a folder for this tutorial along with a new file called setup.py inside. We'll use the setuptools package to help describe our new starfox package.

"""starfox setup.py
"""
from setuptools import setup, find_packages

setup(
    name="starfox",
    version="0.1.0",
    description="Do a barrel roll!",
    packages=find_packages(),
    install_requires=[
        "click",
        "questionary"
    ],
    entry_points={
        'console_scripts': ['starfox=starfox.main:main']
    },
    author="Justin Swaney",
    license="MIT"
)

Some things to note about this setup function call:

The find_packages function will auto-discover our package
The install_requires argument specifies minimum dependencies to run
The entry_points argument allows us to name our CLI starfox and bind it to a function called main within the starfox/main.py module (which we still need to write)

Make a folder called "starfox/" next to our setup.py file. This folder will hold all of the source code for our new package. Inside the starfox folder, create an empty file called __init__.py. The existence of this file will cause Python to consider our starfox folder as a package.

This may seem dubious with our simple example, but an __init__.py file allows us to namespace our package modules and execute common initialization code when imported.

Next to the __init__.py file in the starfox folder, create a file called main.py and define a function called main inside.

def main():
    print("Do a barrel roll!")

The starfox.main:main part of the console_scripts option in setup.py is the syntax for referring to this main function.

Installing our new package

We still need to edit our starfox code, but it would be nice if we could test if we've set things up correctly first. Let's just install it in editable mode. This means the package metadata and compiled byte code will be stored within our current folder instead of a distant site-packages folder with the rest of the installed packages on the system.

It's helpful to add .egg-info and __pycache__ to a .gitignore file when working in editable mode so that they do not pollute the git history.

To install starfox in the base environment, we can use pip.

# In the same folder as setup.py
$ pip install -e .

If you want to test this out in a clean environment, you can create a new conda environment and install starfox there. Just remember you'll have to activate your new environment if you want to use the starfox package.

# In the same folder as setup.py
$ conda create -y -n starfox python=3.7
$ conda activate starfox
(starfox) $ pip install -e .

Finally, we can test if our starfox console script works.

$ starfox
Do a barrel roll!

Defining the `starfox` CLI

Now that we have our entry point set up, let's build a CLI that will print random quotes from a given character in the game Star Fox 64. To keep things simple, we will only consider the following characters:

Create a file called quotes.py next to our main.py module. Fill it with these character quotes:

"""Quotes from characters in Star Fox 64"""

FOX = [
    "All aircraft report!",
    "I'll go it alone from here",
    "Sorry to jet, but I'm in a hurry"
]

FALCO = [
    "Enemy group behind us!",
    "AaawwwwWW man, I'm gonna have ta BACK OFF",
    "Hey Einstein, I'm on yourrr siiide!!"
]

SLIPPY = [
    "Don't worry, Slippy's here!",
    "Hold A to charge your laser",
    "This baby can take temperatures up to 9000 degrees!"
]

PEPPY = [
    "It's quiet, TOO quiet...",
    "Do a barrel roll!!!",
    "You've got an enemy on your tail!"
]

GENERAL = [
    "It's about time you showed up, Fox. You're the only hope for our world!",
    "Recover our base from the enemy army",
    "Star Fox, we are in your debt"
]

ANDROSS = [
    "Ahhh, the son of James McCloud",
    "I've been waiting for you, Star Fox",
    "Only I have the brains to rule Lylat!"
]

QUOTES = dict(
    fox=FOX,
    falco=FALCO,
    slippy=SLIPPY,
    peppy=PEPPY,
    general=GENERAL,
    andross=ANDROSS
)

Inside main.py, let's import QUOTES and print QUOTES.keys() instead. We can test that the characters show up in the terminal.

# Inside main.py
from starfox.quotes import QUOTES

def main():
    print(QUOTES.keys())

# In the terminal
$ starfox
dict_keys(['fox', 'falco', 'slippy', 'peppy', 'general', 'andross'])

Now that we have QUOTES available, it would be nice if we could select which character we would like to quote directly from the terminal. We can use the click package to add a character argument to our main function and bind it to a command-line argument.

# Inside main.py
import click
from starfox.quotes import QUOTES

@click.command()
@click.argument('character')
def main(character):
    for k, q in enumerate(QUOTES[character.lower()]): 
        print(f'Quote #{k + 1}: {q}')

If we try running starfox in the terminal, we will get an error because the character argument is a required positional argument (named arguments have dashes in front of them). If we also indicate a character in the terminal, we'll see all quotes for that character:

$ starfox falco
Quote #1: Enemy group behind us!
Quote #2: AaawwwwWW man, I'm gonna have ta BACK OFF
Quote #3: Hey Einstein, I'm on yourrr siiide!!

The last thing we need to do is randomly sample from these quotes. We can use the random.sample built-in to do this.

# Inside main.py
from random import sample
import click
from starfox.quotes import QUOTES

@click.command()
@click.argument('character')
def main(character):
    qs = QUOTES[character.lower()]
    q = sample(qs, k=1)[0]
    print(f'{character.upper()} - "{q}"')

Now we'll get random quotes from the given character.

$ starfox slippy
SLIPPY - "Don't worry, Slippy's here!"
$ starfox slippy
SLIPPY - "This baby can take temperatures up to 9000 degrees!"

Creating a command-line wizard

Now that we can print random quotes, let's simulate a conversation between Star Fox 64 characters that the user selects. Ideally, we would re-use the starfox CLI code for this conversation feature and still support the previous random quote feature. Fortunately we can do this with @click.option, which will make our character input an optional named argument. This is convenient because we can alter the behavior of the starfox CLI based on whether or not the character option is provided (not None).

For our simulated conversation, we can prompt the user for which characters to include and how many quotes from each to sample. With that in mind, it would be convenient to re-use the random sampling logic we have already written for this as well as the above mentioned case when character is provided. Putting this all together, our main.py file would look like this.

from random import sample
import click
import questionary
from starfox.quotes import QUOTES

def quote_character(character):
    """Prints a random quote from a given Star Fox character"""
    qs = QUOTES[character.lower()]
    q = sample(qs, k=1)[0]
    print(f'{character.upper()} - "{q}"')

@click.command()
@click.option('-c', '--character')
def main(character):
    """The `starfox` CLI"""
    if character:
        quote_character(character)
        return
    answers= questionary.form(
        characters = questionary.checkbox("Select characters", choices=QUOTES.keys()),
        n_iter = questionary.text("How long do you want the conversation to be? (int)")
    ).ask()
    for _ in range(int(answers['n_iter'])):
        for c in answers['characters']: 
            quote_character(c)

Notice that we have factored out the random quote logic to a function called quote_character. We then use this function when the user provides a character as well as in our simulated conversation. We have changed our @click.argument to a @click.option to make it an optional named argument. We also return before getting to the simulated conversation when character is specified. If the user were to run starfox in the terminal, character would be None and we would skip right to the questionary form.

questionary is a Python package for making command-line wizards. Here, we use it to prompt the user for two pieces of information: which characters to inlcude in our conversation and how many quotes to sample from each character. Once we get this information, we simply use our quote_character function in a loop that goes around the horn.

Let's check that we didn't break anything from before.

$ starfox -c general
GENERAL - "It's about time you showed up, Fox. You're the only hope for our world!"

Looks good. Now let's try our new wizard.

$ starfox

? Select characters (Use arrow keys to move, <space> to select, <a> to toggle, <i> to invert)
   ○ fox
   ● falco
   ● slippy
   ● peppy
   ○ general
 » ● andross

? How long do you want the conversation to be? (int) 2

FALCO - "Hey Einstein, I'm on yourrr siiide!!"
SLIPPY - "This baby can take temperatures up to 9000 degrees!"
PEPPY - "Do a barrel roll!!!"
ANDROSS - "Ahhh, the son of James McCloud"
FALCO - "Hey Einstein, I'm on yourrr siiide!!"
SLIPPY - "This baby can take temperatures up to 9000 degrees!"
PEPPY - "You've got an enemy on your tail!"
ANDROSS - "Only I have the brains to rule Lylat!"

Nice! I suggest reading your results out loud to get the full effect.

Final thoughts

In this tutorial, we learned how to setup a Python package with a CLI called starfox. We used click to bind command-line arguments to function arguments. We also used questionary to create a command-line wizard that simulates a conversation between characters in Star Fox 64.

Do a barrel roll! You deserve it!

Source Availability

All source materials for this article are available here on my blog GitHub repo.

Speeding up scientific computing with multiprocessing in Python

Justin Swaney — Mon, 01 Mar 2021 17:09:51 +0000

Photo by Christian Wiediger on Unsplash

In this tutorial, we will look at how we can speed up scientific computations using multiprocessing in a real-world example. Specifically, we will detect the location of all nuclei within fluorescence microscopy images from the public MCF7 Cell Painting dataset released by the Broad Institute.
After completing this tutorial, you will know how to do the following:

Detect nuclei in fluorescence microscopy images
Speed up image processing using multiprocessing.Pool
Reduce serialization overhead using copy-on-write (Unix only)
Add a progress bar to parallel processing tasks

Download the dataset

Note: The imaging data is nearly 1 GB for just one plate. You will need a few GB of space to complete this tutorial.

First, download the imaging data for the first plate of the Cell Painting dataset here to a folder called images/ and extract it. You should get a folder called images/Week1_22123/ filled with TIFF images. The images are named according to the following convention:

# Template
{week}_{plate}_{well}_{field}_{channel}{id}.tif

# Example
Week1_150607_B02_s4_w1EB868A72-19FC-46BB-BDFA-66B8966A817C.tif

Find all DAPI images

DAPI is a fluorescent dye that binds to DNA and will stain the nuclei in each image. Typically, DAPI signal will be captured in the first channel of the microscope. We'll first need to parse the image names to find images from the first channel (w1) so that we can detect the nuclei from only the DAPI images.

from glob import glob

paths = glob('images/Week1_22123/Week1_*_w1*.tif')
paths.sort()
print(len(paths))  # => 240

Here, we used glob and a pattern with wildcards to find the paths to all DAPI images. Apparently there are 240 DAPI images that we need to process! We can load the first image to see what we are working with.

from skimage import io
import matplotlib.pyplot as plt

img = io.imread(paths[0])
print(img.shape, img.dtype, img.min(), img.max())
# => (1024, 1280) uint16 176 10016

plt.figure(figsize=(6, 6))
plt.imshow(img, cmap='gray', vmax=4000)
plt.axis('off')
plt.show()

We can see that the images are (1024, 1280) arrays of 16-bit unsigned integers. This particular image has pixel intensities ranging from 176 to 10,016.

Using `skimage` to detect nuclei

Scikit-image is a Python package for image processing that we can use to detect nuclei in these DAPI images. A simple and effective approach is to smooth the image to reduce noise and then find all local maxima that are above a given intensity threshold. Let's write a function to do this on our example image.

from skimage.filters import gaussian
from skimage.feature import peak_local_max

def detect_nuclei(img, sigma=4, min_distance=6, threshold_abs=1000):
    g = gaussian(img, sigma, preserve_range=True)
    return peak_local_max(g, min_distance, threshold_abs)

centers = detect_nuclei(img)
print(centers.shape)  # => (214, 2)

plt.figure(figsize=(6, 6))
plt.imshow(img, cmap='gray', vmax=4000)
plt.plot(centers[:, 1], centers[:, 0], 'r.')
plt.axis('off')
plt.show()

Now we're getting somewhere! The detect_nuclei function takes in an img array and returns an array of nuclei (x, y)-coordinates called centers. It first smooths the img array by applying a gaussian blur with sigma = 4. The preserve_range = True argument prevents skimage from normalizing the image so that we get back a smoothed image g in the same intensity range as the input. Then, we use the peak_local_max function to detect all local maxima in our smoothed image. We set min_distance = 4 to prevent any two nuclei centers from being unrealistically close together and threshold_abs = 1000 to avoid detecting false positives in dark regions of the image.

Note: We will use this detect_nuclei function in all subsequent experiments.

Method 1 - List comprehension with IO

Now that we have a function to detect nuclei coordinates, let's apply this function to all of our DAPI images using a simple list comprehension.

from tqdm.notebook import tqdm

def process_images1(paths):
    return [detect_nuclei(io.imread(p)) for p in paths]

meth1_times = %timeit -n 4 -r 1 -o centers = process_images1(tqdm(paths))
# => 18 s ± 0 ns per loop (mean ± std. dev. of 1 run, 4 loops each)

In process_images1, we take in a list of image paths and use a list comprehension to load each image and detect nuclei. Using the %timeit magic command in Jupyter, we can see that this approach has an average execution time of ~18 seconds.

Method 2 - multiprocessing.Pool with IO

Let's see if we can speed things up with a processing Pool. To do this, we need to define a function to wrap detect_nuclei and io.imread together so that we can map that function over the list of image paths.

import multiprocessing as mp

def _process_image(path):
    return detect_nuclei(io.imread(path))

def process_images2(paths):
    with mp.Pool() as pool:
        return pool.map(_process_image, paths)

meth2_times = %timeit -n 4 -r 1 -o centers = process_images2(paths)
# => 5.54 s ± 0 ns per loop (mean ± std. dev. of 1 run, 4 loops each)

Much better. This was run on a machine with cpu_count() == 8. Now the same task has an average execution time of ~5.5 seconds (>3x speed up).

Method 3 - List comprehension from memory

What would happen if we read all the images into memory before detecting nuclei? Presumably, this approach would be faster because we wouldn't have to read the image data from disk during the %timeit command. Let's try it.

import numpy as np

images = np.asarray([io.imread(p) for p in tqdm(paths)])

def process_images3(images):
    return [detect_nuclei(img) for img in images]

meth3_times = %timeit -n 4 -r 1 -o centers = process_images3(tqdm(images))
# => 17.7 s ± 0 ns per loop (mean ± std. dev. of 1 run, 4 loops each)

Now the list comprehension is slightly faster, but only by about 0.3 seconds. This suggests that reading image data is not the rate-limiting step in Method 1, but rather the detect_nuclei computations.

Method 4 - multiprocessing.Pool from memory

Just for completeness, let's try using multiprocessing.Pool to map our detect_nuclei function over all the images in memory.

def process_images4(images):
    with mp.Pool() as pool:
        return pool.map(detect_nuclei, images)

meth4_times = %timeit -n 4 -r 1 -o centers = process_images4(images)
# => 6.23 s ± 0 ns per loop (mean ± std. dev. of 1 run, 4 loops each)

Wait, what? Why is this slower than also reading the images from disk with multiprocessing? If you read the documentation for the multiprocessing module closely, you will learn that data passed to workers of a Pool must be serialized via pickle. This serialization step creates some computational overhead. That means for Method 2, we were pickling path strings, but in this case, we were pickling entire images.

Method 5 - multiprocessing.Pool with serialization fix

Fortunately, there are several ways to avoid this serialization overhead when using multiprocessing. On Mac and Linux systems, we can take advantage of how the operating system handles process forks to efficiently process large arrays in memory (sorry Windows 🤷‍♂️). Unix-based systems use copy-on-write behavior for forked processes. Loosely, copy-on-write means the forked process will only copy data when attempting to modify the shared virtual memory. This all happens behind the scenes, and copy-on-write is sometimes referred to as implicit sharing.

def _process_image_memory_fix(i):
    global images
    return detect_nuclei(images[i])

def process_images5(n):
    with mp.Pool() as pool:
        return pool.map(_process_image_memory_fix, range(n))

meth5_times = %timeit -n 4 -r 1 -o centers = process_images5(len(paths))
# => 5.31 s ± 0 ns per loop (mean ± std. dev. of 1 run, 4 loops each)

Nice! This is now slightly faster than Method 2 (as we may have originally expected). Notice that we use a global statement here to explicitly declare that this function uses the images array defined globally. While this is not strictly required, I have found it helpful to indicate that the _process_image_memory_fix function depends on some images array being available. We then map over all indices and allow each process to access the image it needs by indexing into the images array. This approach will only pickle integers instead of the images themselves.

Results

Let's compare the average execution times of all 5 methods.

Overall, multiprocessing drastically reduced the average execution times for this nuclei detection task. Interestingly, a naive approach to multiprocessing with images in memory actually increased the average execution time. By taking advantage of copy-on-write behavior, we were able to remove the significant serialization overhead of incurred when pickling whole images.

To show off more of our results, let's show a random sample of 16 DAPI images sorted by the number of detected nuclei in each image.

It appears that our simple nuclei detection strategy was effective and has allowed us to organize the DAPI images according to the observed cell density in each field of view.

Conclusion

From Spaceballs the movie

I really hope this article was helpful, so please let me know if you enjoyed it. The multiprocessing techniques presented here have helped me achieve ridiculous speed when processing images for our SCOUT paper. Perhaps in the future we can revisit this example but with GPU acceleration to finally achieve ludicrous speed.

Bonus: Using progress bars with multiprocessing.Pool

In the above multiprocessing examples, we do not have a nice tqdm progress bar. We can get one if we use pool.imap instead, which is useful for long-running computations.

def _process_image(path):
    return detect_nuclei(io.imread(path))

def process_images_progress(paths):
    with mp.Pool() as pool:
        return list(tqdm(pool.imap(_process_image, paths), total=len(paths)))

centers = process_images_progress(paths)

References

We used image set BBBC021v1 [Caie et al., Molecular Cancer Therapeutics, 2010], available from the Broad Bioimage Benchmark Collection [Ljosa et al., Nature Methods, 2012].

Source Availability

All source materials for this article are available here on my blog GitHub repo.

10 powerful built-in functions from the Python standard library

Justin Swaney — Sun, 28 Feb 2021 22:26:43 +0000

Photo by Divide By Zero on Unsplash

There are a ton of awesome packages available in the Python ecosystem, but sometimes all you need is a built-in function. It might not be as glamorous without the branding of an exciting new package, but you can go a long way with the standard library. If nothing else, learning about the following Python built-ins could help you speed up your development and improve code readability.

Note: this tutorial uses some Python 3.7+ syntax

Globally available functions

#1 - enumerate

The enumerate function is handy when you want to keep track of an index as you loop through an iterable. This can remove counter variables you might have been using if coming to Python from another language.

abcs = list('ABCDEF')

# Without enumerate
k = 0
for letter in abcs:
    print(k, letter.lower())  # => 0 a, 1 b ...
    k += 1

# With enumerate
for k, letter in enumerate(abcs): 
    print(k, letter.lower())  # => 0 a, 1 b ...

The way enumerate works is by returning an enumerate object with a __next__ method that returns a tuple of (count, value). You can even use destructuring to enumerate more complex values.

abcs = [('a', 'A'), ('b', 'B'), ('c', 'C')]

for k, (lower, upper) in enumerate(abcs):
    print(k, lower, upper)  # => 0 a A, 1 b B ...

#2 - zip

The zip function is useful when you want to combine multiple iterables together. A common use-case is to loop over two lists simultaneously without direct indexing.

ABCs = list('ABCDEF')
abcs = list('abcdef')

# Without zip
for upper in ABCs:
    idx = ABCs.index(upper)
    print(upper, abcs[idx])  # => A a, B b ...

# With zip
for upper, lower in zip(ABCs, abcs):
    print(upper, lower)  # => A a, B b ...

Removing the manual indexing makes the code easier to understand because it removes a step of thinking about what the index corresponds to. If you're feeling fancy, you can combine enumerate with zip because enumerate doesn't care what kind of iterable it recieves.

for k, (upper, lower) in enumerate(zip(ABCs, abcs)):
    print(k, upper, lower)  # => 0 A a, 1 B b ...

#3 - map

The map function is powerful in its generality, and understanding it can help you speed up computations with multiprocessing. Conceptually, map applies a given function to each element in an iterable and returns the results in order. We can use lambda to quickly create the function f(x) = abs(x ** 2 - 4 * x) and map it over some values of x.

func = lambda x: abs(x ** 2 - 4 * x)
xs = range(10)

# Without map
results = [func(x) for x in xs]
print(results)  # => [0, 3, 4, 3, ...]

# With map
results = list(map(func, xs))
print(results)  # => [0, 3, 4, 3, ...]

With map, it's clear that we are trying to apply func to all xs but without specifing exactly how to do it. For example, if xs were a large vector, we could speed up this computation using the analogous map method on a multiprocessing.Pool object.

from multiprocessing import Pool, cpu_count

def func(x): 
    return abs(x ** 2 - 4 * x)

xs = range(1_000_000)
with Pool(cpu_count()) as pool:
    results = pool.map(func, xs)
print(len(results))  # => 1000000

Notice we are using the with keyword to create a context for managing resouces for the processing pool. The cpu_count function will return the number of cores available on the machine, and pool.map will apply func to items in xs in parallel. This demonstrates that map also provides a useful syntax for expressing a single program, multiple data (SPMD)-style computation.

#4 - dir

The dir function is helpful to perform introspection on Python objects (to see what attributes they have). For example, we can use dir to see what's defined in the os module.

import os
print(dir(os))  # => ['CLD_CONTINUED', 'CLD_DUMPED', ...]

We can also use dir to observe the special methods that underlie the Python data model.

print(dir(list('ABCDEF')))  # => ['__add__', '__class__', ...]

Not only can we see methods like append and sort on our list, but we can also see the __iter__ method, which returns an iterator object when involked.

The string module

#5 - string.ascii_uppercase

The string module contains useful constants that can save you some typing. For example, if you ever need all the letters of the English alphabet, you can use string.ascii_lowercase or string.ascii_uppercase.

import string
print(string.ascii_lowercase)
print(string.ascii_uppercase)

Personally, I have used string.ascii_uppercase when working with 384-well plates that have well names from "A1" to "P24"

The itertools module

As the name suggests, the itertools module contains utilities for working with iterables. These iterator building-blocks not only provide memory efficiency through lazy execution, but they also provide a clean syntax for several common transformations.

#6 - itertools.cycle

The itertools.cycle function takes an iterable and returns an iterator that will endlessly loop through the original input values.

import numpy as np
import matplotlib.pyplot as plt
from itertools import cycle

n_lines = 10
n_pts = 10
styles = [ 'k-', 'b--', 'r-', 'g--' ]

plt.figure()
for k, style in zip(range(n_lines), cycle(styles)):
    x = np.arange(n_pts)
    y = np.sqrt(x) + np.random.rand(n_pts)
    plt.plot(x, y, style, label=f'Line {k}: {style}')
plt.legend()
plt.show()

As you can see, this can be useful for cycling through plotting styles where it's cumbersome to specify a ton of styles and okay to repeat some. Notice that cycle was used inside zip, so the cycled styles iterator didn't exhaust even though len(styles) < n_lines.

#7 - itertools.product

The itertools.product function returns the Cartesian product of the input iterables, which is just a fancy way of saying it makes a grid.

import string
from itertools import product

row_names = string.ascii_uppercase[:16]
col_names = range(1, 25)

wells = list(product(row_names, col_names))
print(wells)  # => [('A', 1), ('A', 2), ...]

Here, we get a list of 384 tuples (16 rows * 24 columns) that cycle like an odometer through the input iterables.

The functools module

The functools module contains functions that act on other functions. You can find lots of great functional programming utilities here.

#8 - functools.partial

The functools.partial function allows us to Curry a function, which means to specify certain arguments as fixed values ahead of time. Currying a function returns a new function that takes fewer arguments.

from functools import partial

# f takes 2 arguments
def f(x, y):
    return abs(x ** 2 - 4 * y)

# g takes 1 argument
g = partial(f, y=0)

print(g(4))  # => 16

This will print 16 because g(x) = f(x, y=0) = abs(x ** 2). This is particularly useful when working with higher-order functions that assume a certain signature for an input function (ie. using pool.map for multiprocessing)

Useful IO modules

Including these built-in modules might be cheating since they are not functions themselves. However, they provide simple interfaces for reading and writing different types of data, so it seems more logical to consider these functions together.

#9 - json

The json module provides a built-in interface for reading and writing JSON data.

import json

data = {
    'name': 'Bia',
    'color': 'Black',
    'breed': 'Labrador mix'
}

with open('data.json', 'w') as fd:
    json.dump(data, fd)  # => data.json JSON file created

with open('data.json', 'r') as fd:
    data_json = json.load(fd)

print(data_json)  # => {'name': 'Bia', ...}

This is particularly useful when dealing with data from the web, doing your own scraping, or working with other languages that work well with JSON.

#10 - pickle

The pickle module provides a built-in interface for reading and writing pickled objects. Pickling is the process of serializing Python objects into a flat, binary structure that can be written to disk. The pickle module is useful, therefore, to save an object in its current state for later use.

import pickle

model = {
    'method': 'logistic',
    'weights': [-0.34, 0.45, 0.72, -1.21]
}

with open('model.pkl', 'wb') as fd:
    pickle.dump(model, fd, pickle.HIGHEST_PROTOCOL)  # => model.pkl created

with open('model.pkl', 'rb') as fd:
    model_pkl = pickle.load(fd)

print(model_pkl)  # => {'method': 'logistic', ...}

This use-case is common in machine learning applications where models are Python objects and weights are internal state that must be saved.

Final thoughts

The functions discussed here are just a few examples of the awesome built-in functions available in Python. Next time you find yourself trying to figure out how to write something in a clear and concise way, give the Standard Library some thought. It has helped me many times, and I hope it can help you too.

DEV Community: Justin Swaney

Computing pretrained image embeddings in Python

Setting up the embedding package

Developing the embedding CLI

Creating a simple entrypoint

Loading image data

Linear embeddings with PCA

Deep embeddings using ResNet50V2

Saving image embeddings

Loading MNIST data

Generating UMAP visualizations

Removing any number of layers

Results

Final thoughts

Source availability

Tips for using Jupyter Notebooks with GitHub

Tip #1: Auto-formatting in Jupyter notebooks

Tip #2: VSCode notebook diffs

Tip #3: nbstripout as git filter

Tip #4: Executable notebooks with Papermill

Final thoughts

Source availability

Setting up Python Linting and Auto-formatting in VSCode

Setup

Linting with pylint

Auto-formatting with black

Result

Final thoughts

Source availability

Poetry: a yarn-like package manager for Python

Installing the poetry package manager

Creating our math-demo package

But what about conda?

Installing our dependencies

Catenary curves and their parabolic approximations

Implementing the exact and approximate catenary calculations

Writing a failing test

Making the tests pass

Adding a parabolic approximation

Visualizing exact and approximate catenary curves

Final thoughts

Source Availability

Creating beautiful command-line interfaces in Python

Setting up the starfox Python package

Installing our new package

Defining the starfox CLI

Creating a command-line wizard

Final thoughts

Source Availability

Speeding up scientific computing with multiprocessing in Python

Download the dataset

Find all DAPI images

Using skimage to detect nuclei

Method 1 - List comprehension with IO

Method 2 - multiprocessing.Pool with IO

Method 3 - List comprehension from memory

Method 4 - multiprocessing.Pool from memory

Method 5 - multiprocessing.Pool with serialization fix

Results

Conclusion

Bonus: Using progress bars with multiprocessing.Pool

References

Source Availability

10 powerful built-in functions from the Python standard library

Globally available functions

#1 - enumerate

#2 - zip

#3 - map

#4 - dir

The string module

#5 - string.ascii_uppercase

The itertools module

#6 - itertools.cycle

#7 - itertools.product

The functools module

#8 - functools.partial

Useful IO modules

#9 - json

#10 - pickle

Final thoughts

Setting up the `embedding` package

Developing the `embedding` CLI

Linting with `pylint`

Auto-formatting with `black`

Installing the `poetry` package manager

Creating our `math-demo` package

But what about `conda`?

Setting up the `starfox` Python package

Defining the `starfox` CLI

Using `skimage` to detect nuclei