# Multiprocessing in Python – all about pickling

Chances are you heard that multiprocessing in Python is hard. That it takes time and, actually, don’t even try because there’s something like global interpreter lock (GIL), so it isn’t even true parallel execution. Well, GIL is true, but the rest is a lie. Multiprocessing in Python is rather easy.

One doesn’t have to look far to find nice introductions into processing in Python [link1, link2]. These are great and I do recommend reading on them. Even first Google result page should return some comprehensible tutorials. However, what I was missing from these tutorials is some information about handling processing within class.

Multiprocessing in Python is flexible. You can either define Processes and orchestrate them as you wishes, or use one of excellent methods herding Pool of processes. By default Pool assumes number of processes to be equal to number of CPU cores, but you can change it by passing processes parameter. Main methods included in Pool are apply and map, which let you run process with arbitrary arguments or execute parallel map, respectively. There are also asynchronous versions of these, i.e. apply_asyncand map_async.

Quick example:

from multiprocessing import Pool
def power(x, n=10):
return x**n

pool = Pool()
pow10 = pool.map(power, range(10,20))
print(pow10)

[10000000000, 25937424601, 61917364224, 137858491849, 289254654976, 576650390625, 1099511627776, 2015993900449, 3570467226624, 6131066257801]


Simple, right? Yes, this is all what’s needed. Now, go and use multiprocessing!

Actually, before you leave to follow your dreams there’s a small caveat to this. When executing processes Python first pickles these methods. This create a bottleneck as only objects that are pickle will be passed to processes. Moreover, Pool doesn’t allow to parallelize objects that refer to the instance of pool which runs them. It sounds convoluted so let me exemplify this:

from multiprocessing import Pool

class BigPow:
def __init__(self, n=10):
self.pool = Pool()

self.n = n

def pow(self, x):
return x**self.n

def run(self, args):
#pows = self.pool.map(ext_pow, args)
pows = self.pool.map(self.pow, args)
return sum(pows)

def ext_pow(x):
return x**10

if __name__ == "__main__":
big_pow = BigPow(n=10)
pow_sum = big_pow.run(range(20))
print(pow_sum)


Code above doesn’t work, unless we replace self.pow with ext_pow. This is because self contains pool instance. We can remove that through removing pool just before pickling through __getstate__ (there’s complimentary function __setstate__ to process after depickling).

from multiprocessing import Pool

class BigPow:
def __init__(self, n=10):
self.pool = Pool()

self.n = n

def pow(self, x):
return x**self.n

def run(self, args):
pows = self.pool.map(self.pow, args)
return sum(pows)

def __getstate__(self):
self_dict = self.__dict__.copy()
del self_dict['pool']
return self_dict

if __name__ == "__main__":
big_pow = BigPow(n=10)
pow_sum = big_pow.run(range(20))
print(pow_sum)


This is good, but sometimes you’ll get an error stating something like “PicklingError: Can’t pickle : attribute lookup __builtin__.instancemethod failed”. In such case you have to update registry for pickle on what to actually goes into pickling.

from multiprocessing import Pool

#######################################
import sys
import types
#Difference between Python3 and 2
if sys.version_info[0] < 3:
import copy_reg as copyreg
else:
import copyreg

def _pickle_method(m):
class_self = m.im_class if m.im_self is None else m.im_self
return getattr, (class_self, m.im_func.func_name)

copyreg.pickle(types.MethodType, _pickle_method)
#######################################

class BigPow:
def __init__(self, n=10):
self.pool = Pool()

self.n = n

def pow(self, x):
return x**self.n

def run(self, args):
pows = self.pool.map(self.pow, args)
return sum(pows)

def __getstate__(self):
self_dict = self.__dict__.copy()
del self_dict['pool']
return self_dict

if __name__ == "__main__":
big_pow = BigPow(n=10)
pow_sum = big_pow.run(range(20))
print(pow_sum)


Yes, this is all because of Pickle. What can you do with it? In a sense, not much, as it’s the default battery-included solution. On the other, pickle is generally slow and now community standard seems to be dill. It would be nice if something was using dill instead of pickle. Right? Well, we are in luck because pathos does exactly that. It has similar interface to multiprocessing is it’s also easy to use. Personally I’m Ok with multiprocessing, as I like not to import too much, but this is a good sign that there are developments towards more flexible solutions.

# Complete Ensemble EMD with Adaptive Noise (CEEMDAN) in Python

CEEMDAN is available in Python through PyEMD.

There are as many empirical mode decomposition (EMD) variations as many teams are working on it. Everyone notices that in general EMD is very helpful method, yet, there’s room for improvement. Few years back I have stopped doing modifications myself in exchange for working on mathematically sound model of coupled oscillator. Nevertheless, since I spent quite a lot of time on EMDs and have enjoy playing with it, from time to time something will catch my eye. This is what happened with Complete Ensemble EMD with Adaptive Noise (CEEMDAN).

As name suggests this is an expansion on the ensemble EMD, which was already covered. What’s the difference? In case of CEEMDAN we’re also decomposing our perturbation to the system, i.e. added noise. Method creates an ensemble of many perturbations, decomposes them using EMD and resulting IMFs are included to evaluate components of the input. I will refer to these components as cIMF. The actual algorithm was first proposed by Torres et. al [1], but shortly after an improvement in efficiency was proposed[2]. These updates refer mainly to noticing that for their purpose one doesn’t need to compute all IMFs and weight parameter can be progressively scaled as well. What exactly is this algorithm?

Let’s define operator IMFi(S) which returns ith IMF (pure EMD) of its input S and M(S) to provide local mean, i.e. M(S) = S – IMF1(S), then the algorithm is as follows:

1. Create Gaussian noise ensemble W={wi}, where $i\in[1..N]$, and decompose them using EMD.
2. For input signal S calculate grand average of local means from signal perturbed by scaled noise first IMF
$R_{1} = \frac{1}{N}\sum_{1}^{N} M(S+ \beta_0 IMF_{1}(w^{i}))$.
3. Assign first cIMF to be: C1 = S – R1.
4. Compute $R_{k}= \frac{1}{N} \sum_{i=1}^{N} M(R_{k-1} + \beta_{k-1} IMF_{k}(w^{i}))$.
5. Calculate kth cIMF as Ck = Rk-1 – Rk.
6. Iterate 4. and 5. until set of {S, Ck} fulfils EMD stopping criteria.

As it can be seen a family of parameters β has been included in the algorithm. These scalars refer to the amount of decomposed noise used to compute cIMFs. This is what the authors refer to as noise adaptive. These parameters are arbitrary, but it’s suggested in improved version [2] to set them as $\beta_{k}= \epsilon_{0} \sigma(R_k)$, where σ is standard divination of argument and ε is another arbitrary parameter. Looking at point 4. one can see that for ith residue we are using ith IMF computed from noise. This is a problem, because EMD decomposes signal into a finite set of components and it can happen that there isn’t ithIMF. In this case authors are suggesting to assume component to be equal 0.

Advantage of this variation comes from the fact that created decomposition {Ci} fully reconstructs input. This is in contrast to EEMD which doesn’t guarantee such completeness. However, with CEEMDAN questions rise regarding the meaning of added scaled IMFs of noise. Augmenting signal with ensemble of pure noise creates perturbations of input without any distinguished direction. As it has been observed by Flandrin et al. [3] when decomposing white noise EMD acts as a dyadic filter bank. This means that extracted IMFs will have preferred structure and adding them to input will be similar to adding vector with random length but particular direction.

Regardless of all, CEEMDAN is definitely an interesting method. Just purely by the number of citations it seems that I’m not the only one thinking that. I’ve included it to my Python PyEMD package, so feel free to play with it and leave some feedback.

References

[1] Torres ME, Colominas MA, Schlotthauer G, Flandrin P. A complete ensemble empirical mode decomposition with adaptive noise. InAcoustics, speech and signal processing (ICASSP), 2011 IEEE international conference on 2011 May 22 (pp. 4144-4147). IEEE.
[2] Colominas MA, Schlotthauer G, Torres ME. Improved complete ensemble EMD: A suitable tool for biomedical signal processing. Biomedical Signal Processing and Control.
[3] Flandrin P, Rilling G, Goncalves P. Empirical mode decomposition as a filter bank. IEEE signal processing letters.

# Python Empirical Mode Decomposition on Image

One of the packages I intend long term maintain and support is Python implementation of Empirical Mode Decomposition (EMD) called PyEMD. I will skip introduction of the method as it has been explained in few other posts [1, 2, 3, …]. This blog entry is more about announcement of new feature which also means new version.

PyEMD version 0.2 is out. This means that PyEMD now supports 2D data (image) decomposition. Other visible improvements include documentation and more thorough testing both of code and data cases. Installation instructions are provided on the project’s webpage.

I am more than happy to include other improvements or suggestions. The next big step will be support for 3D and multi dimensional data. Please get in touch if you feel that there is something missing.

Image decomposition is based on the simple extremum definition: a point that is above (max) or below (min) surrounding. Behind the hood this is done using SciPy’s ndim maximum_filter. These are then connected using SmoothBivariateSpline. Stopping criteria can be chosen to be either based on the number of sifting operations or threshold values for mean and standard deviations.

Below is included exemplary decomposition, with the top image being input and the following two are the outputs. Exact formula with which the image was generated is
$sin(4\pi \cdot x) \cdot \left( cos(8\pi y + 8\pi x) + 3\right) + \left(5x + 2y - 0.4 \right) y + 2$. Python code generating this example is in provided in documentation in Examples/EMD2D.

# Installing Cartopy on Ubuntu 14.04 (or Travis-CI)

I’m becoming increasingly convinced that using GitHub to share projects, even in work-in-progress (WIP) state is quite beneficial. For one, someone can quickly point out that such project exists or suggest a tool which would simplify some operations. What’s more, with tools like Travis-CI or CodeCov it’s easy to follow whether the project is actually buildable or whether updates broke functionally of other features. Real convenience! Although, there is some price to pay for these goods.

One of the projects I’m currently working on is MapViz. Its aim is to provide easy to use visualization and summary for country and other administrative units data, mainly EuroStat dataset. This project heavily depends on Cartopy library – a great tool for handling geographical data, however, not so great when it comes to installing. It comes with many dependencies, some of which have to be installed separately. For example, the newest version of Cartopy requires Proj.4 lib in version 4.9, which is not available for Ubuntu older than 16.04. This is unfortunate, because Travis allows to use Ubuntu in versions 12.04 (precise) or 14.04 (trusty). For these only Proj.4 in version 4.8 is available and that’s not enough. Once these are obtained, installing Cartopy with pip is easy, it’s just:

$pip install cartopy  but the trick is to install all dependencies. The workaround is dirty, but it works. One can directly download and install Proj.4 and its dependencies (libproj9) from 16.04 (xenial) repo. For exact links go to libproj9 and libproj-dev, select your architecture and then click on any mirror link. (Just a note that by default Travis-CI uses amd64.) So, download it with wget and then install with dpkg. Here’s an example: $ wget http://es.archive.ubuntu.com/ubuntu/pool/universe/p/proj/libproj9_4.9.2-2_amd64.deb
$sudo dpkg -i libproj9_4.9.2-2_amd64.deb  For comparison, here is my whole .travis.yml file: dist: trusty language: python python: - 2.7 - 3.4 - 3.5 before_install: - sudo apt-get -qq update - sudo apt-get install -y libproj-dev proj-bin proj-data - sudo apt-get install -y python-pyproj - sudo apt-get install -y libc6 - wget http://es.archive.ubuntu.com/ubuntu/pool/universe/p/proj/libproj9_4.9.2-2_amd64.deb - sudo dpkg -i libproj9_4.9.2-2_amd64.deb - wget http://es.archive.ubuntu.com/ubuntu/pool/universe/p/proj/libproj-dev_4.9.2-2_amd64.deb - sudo dpkg -i libproj-dev_4.9.2-2_amd64.deb install: - pip install -r requirements.txt script: - python setup.py install  # Kuramoto in Stan (PyStan) tl;dr: Project on github: https://github.com/laszukdawid/pystan-kuramoto Stan is a programming language focused on probabilistic computations. Although it’s a rather recent language it’s been nicely received in data science/Bayesian community for its focus on designing model, rather than programming and getting stuck with computational details. Depending what is your background you might have heard about it either from Facebook and its Prophet project or as a native implementation for Hamiltonian Monte Carlo (HMC) and its optimised variation – No-U-Turn Sampler for HMC (NUTS). For ease of including models in other programmes there are some interfaces/wrappers available, including RStan and PyStan. Stan is not the easiest language to go through. Currently there are about 600 pages of documentation and then separate “documentations” for wrappers, which for PyStan isn’t very helpful. Obviously there’s no need for reading all of it, but it took me a while to actually understand what goes where an why. The reward, however, is very satisfying. Since I’ve written a bit about Kuramoto models on this blog, it’s consistent if I share its implementation in Stan as well. Pystan-kuramoto project uses PyStan, but the actual Stan code is platform independent. Currently there are two implementations (couldn’t come up with better names): • All-to-all, where Kuramoto model fit is performed to phase vector $\vec{\Phi}$ with distinct oscillators, i.e. $\vec{\Phi}_{N}(t) = \{\phi_1(t), \phi_2(t), \dots, \phi_N(t)\}$. • All-to-one, where the model fits superposition (sum) of oscillators to phase time series $\Phi_{N}(t) = \sum_{n=1}^{N} \phi_n(t)$. In all honesty, this seems to be rather efficient. Optimisation is performed using penalized maximum likelihood estimation with optimization (L-BFGS). Before using it I wasn’t very knowledgeable in the algorithm, but now I’m simply amazed with its speed and accuracy. Thumbs up. # Kuramoto in Python Code for Kuramoto in Python is available here or from code subpage. Explanation on how to use it is on the bottom of this post. ## Tiny introduction Kuramoto[1, 2] is probably one of the most popular and successful models for coupled oscillators. There is plenty of information about it, but in brief summary it models oscillators’ phases to be dependent on scaled phase differences for all pairs of oscillator. A while ago I have actually wrote on a specific method that I used to find optimal parameters for a Kuramoto model given some observation. (See Bayesian Dynamic Inference here.) Mathematically speaking this model is defined by a set of coupled ordinary differential equations (ODEs). Given N oscillators, dynamics for each oscillator’s phase $\phi_i$ is defined as $i\in N$ is defined by $\dot\phi_i = \omega_i + \sum_{j=0}^N k_{ij} \sin(\phi_i-\phi_j)$, where the summation is over all others oscillators. Note that one can leave coupling with itself, because in such case $\sin(\phi_i - \phi_i) = 0$, thus it doesn’t input anything into the final solution. Quick note that Kuramoto coupling model can be used to reference model with coupling function that can be represented as a sum of harmonic functions. For example, second order means that we are including both $k \sin(\Delta\phi)$ and $k^{(2)}\sin(2\Delta\phi)$. In general Kuramoto model of order $M$ would describe $\dot\phi_i = \omega_i + \sum_{m=0}^{M} \sum_{j=0}^N k_{ij}^{(m)} \sin(m (\phi_i-\phi_j))$. ## Examples Below are two examples. Both use the same core values with the exception that for the second system the coupling function defined by two harmonic terms. This naturally changes dynamics, but it shouldn’t change much average value, which is close to intrinsic frequency. Exact values for the first experiment are presented in table below. Each row is a different oscillator with initial phase Φ0 and intrinsic frequency Ω. The following columns denote coupling strength between respective pairs of oscillators. Φ0 ω k.1 k.2 k.3 1 0 28 0.2 1.1 2 π 19 0.5 -0.7 3 0 11 0.3 0.9 Phase dynamics, i.e. time derivative of obtained phases, are presented in Fig. 1. One can see that all plots are centred on intrinsic frequency with some modulations. Fig. 1. Phase dynamics, i.e. time derivative, in simple Kuramoto model. Table below shows values used in the second simulation. Extra 3 columns denote scaling values used in second harmonic. Φ0 ω k(1).1 k(1).2 k(1).3 k(2).1 k(2).2 k(2).3 1 0 28 0.2 1.1 -0.5 0.2 2 π 19 0.5 -0.7 -0.4 1.0 3 0 11 0.3 0.9 0.8 0.8 Again phase dynamics have been presented in a form of plot (Fig. 2). I think the difference is clear. Despite having similar mean values (approximately equal to intrinsic frequency), their modulations have change. Not only the frequency content of these modulations has changed, but also their amplitude. Fig. 2. Phase dynamics, i.e. time derivative, in Kuramoto model with included second harmonic. ## Using code Except for downloading code from either github or code subpage, one is expected to have SciPy module. Kuramoto uses it to solve differential equations. Script below shows an example of how one can use the Kuramoto module. When you run this, make sure that kuramoto.py is either in your path. Note also that most of the code is to just defining initial parameters and plotting the results. Actual execution of the module are two lines. Fig. 1 is the expected output. import numpy as np import pylab as plt from kuramoto import Kuramoto # Defining time array t0, t1, dt = 0, 40, 0.01 T = np.arange(t0, t1, dt) # Y0, W, K are initial phase, intrinsic freq and # coupling K matrix respectively Y0 = np.array([0, np.pi, 0]) W = np.array([28, 19, 11]) K1 = np.array([[ 0, 0.2, 1.1], [0.5, 0, -0.7], [0.3, 0.9, 0]]) K2 = np.array([[ 0, -0.5, 0.2], [-0.4, 0, 1.0], [ 0.8, 0.8, 0]]) K = np.dstack((K1, K2)).T # Passing parameters as a dictionary init_params = {'W':W, 'K':K, 'Y0':Y0} # Running Kuramoto model kuramoto = Kuramoto(init_params) odePhi = kuramoto.solve(T) # Computing phase dynamics phaseDynamics = np.diff(odePhi)/dt # Plotting response nOsc = len(W) for osc in range(nOsc): plt.subplot(nOsc, 1, 1+osc) plt.plot(T[:-1], phaseDynamics[osc]) plt.ylabel("$\dot\phi_{%i}\$" %(osc+1))
plt.show()


[1] Y. Kuramoto, “Chemical Oscillations, Waves, and Turbulence,” 1984, doi: 10.1007/978-3-642-69689-3.
[2] Steven H. Strogatz, “From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators,” 2000, doi: 10.1016/S0167-2789(00)00094-4.

# Update: Particle Swarm Optimisation in Python

It came to my attention that my PSO for Python is actually quite popular. Only after few people contacted me I’ve noticed that the public version was not the same that I’ve been using. It wasn’t bad, but definitely not as good. Thus, obviously, I’ve updated the version and cleaned it a bit.

Update programme is available from my github or from Code subpage.

What’s the difference? There are few.
– Initial values, unless provided, are psuedo-random generated using Halton sequence. This prevents from artificial clustering.
– Perturbation in form of a Gaussian noise should mitigate false local minima by forcing particles to search surrounding area.
– Added max_repetition threshold, which states the maximum number of obtaining the same optimum value. Upon reaching threshold program finishes.
– General improvement in performance.
– Improved usage documentation within the file.
– Program is now compatible with Python3.

Feel free to request any features.

There is an idea of adding progressive save, which would quit, resume and modify parameters at any point of computation.

# Common dtype for NumPy arrays

Recent challenge I had was to convert two given numpy arrays such it’ll be possible to have them in common type without losing information. One can follow this link for stackoverflow entry.

Initially I thought about comparing dtypes, because in Python 2 it is allowed to do something like:

>>> np.float32 < np.float64
True
>>> np.int32 < np.float32
True
>>> np.int64 > np.float16
False


… which kind of makes sense(?). Well, except that int64 vs float16 which looks suspicious. It turns out that these are type objects and Python is simply comparing their locations in memory[Citation needed] and that is obviously not reliable. Actually, in Python 3 such comparison is forbidden and it fails.

As the answer to my stackoverflow question suggests one could try to use dtype.kind and dtype.itemsize to create own ordering. This is not difficult, but it should contain all types, such as (unsigned) ints, floats, bools, complex…

Fortunately, for my purposes, there is NumPy’s method which does what I want, i.e. numpy.find_common_type. It determines common type following standard coercion rules. With a help of this function my common conversion looks like:

import numpy as np

def common_dtype(x, y):
dtype = np.find_common_type([x.dtype, y.dtype], [])
if x.dtype != dtype: x = x.astype(dtype)
if y.dtype != dtype: y = y.astype(dtype)

return x, y


What to expect from such function?
float32 for float16 and float32
float32 for int16 and float32
int64 for int32 and uint16
float64 for int32 and float16

Behaves as it should.

# Halton sequence in Python

Sometimes when we ask for random we don’t actually mean random by just random. Yes, pseudo-random.

Consider unitary distribution with ranges 0 and 1. Say you want to draw 5 samples. Selecting them at random would mean that we might end up with set of {0, 0.1, 0.02, 0.09, 0.01} or {0.11, 0.99, 0.09, 0.91, 0.01}. Yes, these values don’t seem very random, but that’s the thing about randomness, that it randomly can seem to not be random.

Depending on the purpose of our selection, these values might be just OK. After all, they came from that distribution. However, if our goal is to reconstruct the distribution, or extract information about with limited number of samples, it is often better to draw those samples in pseudo-random way. For example, in accordance to van der Corput sequences for 1D distributions or its generalized version Halton sequence.

The best practice for sampling N dimensional distribution is to use different prime numbers for each dimension. For example, when I need to sample a 5 dimensional unitary distribution, or search space, I will use bases of (5, 7, 11, 13, 17). This is to prevent periodic visits of the same position.

In case you are wondering what’s the difference between actual random and pseudo-random, here is a gist:

Both are good, but the actual random can produce many empty holes. What we like to have is a fair representation of all areas of our search space.

Thus, without further ado, here are some code snippets.

This is a definition of my prime generating generator:

def next_prime():
def is_prime(num):
"Checks if num is a prime value"
for i in range(2,int(num**0.5)+1):
if(num % i)==0: return False
return True

prime = 3
while(1):
if is_prime(prime):
yield prime
prime += 2


As for Halton sequence, as mentioned before it uses van der Corput sequence. Again, here is the definition:

def vdc(n, base=2):
vdc, denom = 0, 1
while n:
denom *= base
n, remainder = divmod(n, base)
vdc += remainder/float(denom)
return vdc


And finally, definition for the Halton sequence:

def halton_sequence(size, dim):
seq = []
primeGen = next_prime()
next(primeGen)
for d in range(dim):
base = next(primeGen)
seq.append([vdc(i, base) for i in range(size)])
return seq


To use all of this simply call halton_sequence(size, dim). These variables refer to the number of size of sample poll and the dimension of your problem. So if one wants to sample 3 dimensional space with 10 samples each it would be called as below. (Notice: first dimension has prime value 5, then it’s 7, 11, and following prime values.)

>>> halton_sequence(10, 3)
[
[0, 0.2, 0.4, 0.6, 0.8, 0.04, 0.24000000000000002, 0.44, 0.64, 0.8400000000000001],
[0, 0.14285714285714285, 0.2857142857142857, 0.42857142857142855, 0.5714285714285714, 0.7142857142857143, 0.8571428571428571, 0.02040816326530612, 0.16326530612244897, 0.30612244897959184],
[0, 0.09090909090909091, 0.18181818181818182, 0.2727272727272727, 0.36363636363636365, 0.45454545454545453, 0.5454545454545454, 0.6363636363636364, 0.7272727272727273, 0.8181818181818182]
]


# Matrix Multiplication with Python 3.5

Only recently I have started to use Python 3. It’s been out for good 8+ years and all these excuses about incompatibility with some packages were just lazy. Most packages I use are already ported and if I ever find that something is incompatible… well, I’ll think then. But for now let me pat myself on the back for this great leap, because:

In Python 3.5.3 (released today) there is an operator for matrix multiplication! Check out: PEP 465 — A dedicated infix operator for matrix multiplication. The choice of operator, @, is a bit unfortunate, because of the decorators and general association with reference/internet, but seeing how few possibilities are left it’s probably the best choice.

Yes, this is big news for me. The number of times I confused myself with my own matrix operations is just too damn high! I cannot agree more with the author of the PEP 465, so let my shamelessly copy&paste (paraphrased) his reasoning. Behold!

(…) encounter many mathematical formulas that look like:

S = ( H β r ) T ( H V H T ) − 1 ( H β r )

Here the various variables are all vectors or matrices (details for the curious: [5] ).

Now we need to write code to perform this calculation. In current numpy, matrix multiplication can be performed using either the function or method call syntax. Neither provides a particularly readable translation of the formula:

import numpy as np
from numpy.linalg import inv, solve

# Using dot function:
S = np.dot((np.dot(H, beta) - r).T,
np.dot(inv(np.dot(np.dot(H, V), H.T)), np.dot(H, beta) - r))

# Using dot method:
S = (H.dot(beta) - r).T.dot(inv(H.dot(V).dot(H.T))).dot(H.dot(beta) - r)


With the @ operator, the direct translation of the above formula becomes:

S = (H @ beta - r).T @ inv(H @ V @ H.T) @ (H @ beta - r)


Notice that there is now a transparent, 1-to-1 mapping between the symbols in the original formula and the code that implements it.

Of course, an experienced programmer will probably notice that this is not the best way to compute this expression. The repeated computation of H β r should perhaps be factored out; and, expressions of the form dot(inv(A), B) should almost always be replaced by the more numerically stable solve(A, B) . When using @ , performing these two refactorings gives us:

# Version 1 (as above)
S = (H @ beta - r).T @ inv(H @ V @ H.T) @ (H @ beta - r)

# Version 2
trans_coef = H @ beta - r
S = trans_coef.T @ inv(H @ V @ H.T) @ trans_coef

# Version 3
S = trans_coef.T @ solve(H @ V @ H.T, trans_coef)


Notice that when comparing between each pair of steps, it’s very easy to see exactly what was changed. If we apply the equivalent transformations to the code using the .dot method, then the changes are much harder to read out or verify for correctness:

# Version 1 (as above)
S = (H.dot(beta) - r).T.dot(inv(H.dot(V).dot(H.T))).dot(H.dot(beta) - r)

# Version 2
trans_coef = H.dot(beta) - r
S = trans_coef.T.dot(inv(H.dot(V).dot(H.T))).dot(trans_coef)

# Version 3
S = trans_coef.T.dot(solve(H.dot(V).dot(H.T)), trans_coef)


Readability counts! The statements using @ are shorter, contain more whitespace, can be directly and easily compared both to each other and to the textbook formula, and contain only meaningful parentheses. This last point is particularly important for readability: when using function-call syntax, the required parentheses on every operation create visual clutter that makes it very difficult to parse out the overall structure of the formula by eye, even for a relatively simple formula like this one. Eyes are terrible at parsing non-regular languages. I made and caught many errors while trying to write out the ‘dot’ formulas above. I know they still contain at least one error, maybe more. (Exercise: find it. Or them.) The @ examples, by contrast, are not only correct, they’re obviously correct at a glance.

Again: yes!