Solving problems

Cost of dependencies

We often talk about not re-inventing the wheel and solving a problem yourself can seem useless. However, we must also recognize that pulling in a dependency to solve a problem also carries a huge cost.

I recommend reading this blog post Micro-libraries need to die already about micro-libraries. If you would rather have a live-reading, there is one here by ThePrimeagen.

Another good talk about pitfalls of dependencies is the below from the DevOps conference.

Some common patterns

Below are some common patterns and ideas with respect to solving problems that I think are useful to know. Many of these patterns might look like they are doing the same thing. It looks like this because in many instances they are. For any particular problem the possible ways it can be solved spans a huge space. Hence these patterns will often overlap, not make sense in some situations and fall cleanly into line in others. Learning about these patterns does not mean that you should use them everywhere you can but rather that you now see more possible solutions to your problem and sometimes those fit better with what you are trying to accomplish.

Dependency injection

Dependency injection wikipedia entry

CodeAesthetic - Dependency Injection, The Best Pattern

Composition over inheritance

A very useful concept or pattern is composition over inheritance which solves many of the problems of Object Oriented programming. More on these programming paradigms later.

CodeAesthetic - The Flaws of Inheritance

Adapter

Adapters are useful when you want to plug something external into your own internal logic and this is mostly the case where I have used this pattern.

ArjanCodes - Let’s Take The Adapter Design Pattern To The Next Level

Interface

There are many types of interfaces, the specific one I am thinking about here is using protocols and abstract base classes in Python to specify how something should be implement / behave. I sometimes find the @abstractmethod decorator for the ABC class is quite useful for defining what is needed in case someone wants to extend my code since it will crash directly if not all methods are implemented. Similarly if I am looking for building a very large structure of interaction - I might use Protocol instead.

ArjanCodes - Protocol Or ABC In Python - When to Use Which One?

Another way to think about interfaces is just specifying structure of data rather than code. In this case useful concepts are dataclasses, packing structs, or using Googles protobuf python bindings.

Strategy pattern

I have been using the Strategy pattern for years without even knowing I was doing it: it was just the most modular and readable way to solve many problems I face.

Isaac Harris-Holt - The Strategy Pattern Will Make Your Python Code CLEANER

Bitlytic - Modular Upgrades Made Easy Using the Strategy Pattern

Data oriented design

Although data oriented design is more of a software design philosophy rather than a solution pattern, it is really conducive to thinking about your individual small solutions in a efficient way, and may guide you in a good direction if no clear pattern is visible. The below talk is again at CppCon but the fundamental idea that you should let the data drive your code design rather than how you think about the data, is fundamental to performance. It also ties in to the other patterns listed below.

CppCon 2014 - Mike Acton “Data-Oriented Design and C++”

Vectorisation

Usually when we talk about “vectorisation” we mean one of two different things: SIMD instructions or using vectors in your algorithm.

Translating your problem into tensor/matrix/vector operations

Translating your problem into any vector-like structure and then doing operations on those vectors has special significance in Python as doing loops in Python is very slow due to the nature of interpreted languages. So if we can translate our work into tensor/matrix/vector operations we can offload all looping, processing, casting, etc. to numpy and hence C - no interpretation needed. This also allows for automatic parallelisation across cores, for accessing the speed of C/FORTRAN in Python, and for SIMD (see below) to be utilized fully.

The following is a fantastic introduction into vectorising with numpy and how to use broadcasting to be memory efficient while doing it - a must watch for us scientists!

PyCon UK - Vectorise all the things! How basic linear algebra can speed up your data science code: J Burchell

CPU instruction vectorisation - SIMD

Specific processors can have a special instruction set for doing Single instruction, multiple data (SIMD) operations.

Since we do not have access to direct assembler instructions in Python it does not really have a concept of SIMD instructions, however numpy does! And since a few versions back numpy actually provides a way to look at all the functions it implements see which are optimized: the lib.introspect.opt_func_info function. On my machine I see for example that almost all functions have AVX-512 instruction set available versions, but that all are currently running the regular AVX2 versions.

In C you can pull in immintrin.h and use e.g. the __m256 data type which is a 256-bit vector containing 8 floats instead of using 8 different float data types. In the CPU these 8 floats will be operated on in parallel if compiled towards the correct architecture. More on low-level code inspection in later lectures.

Online-algorithm

An Online algorithm is simply an implementation that can consume the input in sequence without having the entire data available at the start. The concept is incredibly powerful and can be used in everything from something as simple as calculating standard deviations to operating real time systems or handling your computers memory cache.

The following lecture is quite long but introduces the concept of online algorithms and other useful concepts such as cache misses and paging which are important when we look at optimisation, but this is already quite outside the scope of this course so watch at your own peril!

Map-Reduce

The Map-Reduce concept is a way to implement an algorithm such that it is easy to deploy to a super computer cluster trough use of ready made frameworks such as Hadoop and Spark. If a computation is implemented in terms of spark operations (using the Python package), the python code will generate an optimized set of instructions for the super computer cluster to perform on the data. This instruction set can then be submitted to the cluster for execution. Interestingly enough, many algorithms can be implemented in terms of maps and reduces, although this is not guaranteed.

Computerphile - MapReduce

Computerphile - Apache Spark

Visitor

The Visitor pattern is something that might be much more common in regular software development compared to developing research software, especially physics research. Nevertheless, it is a useful mental model to have for abstracting away the algorithm and operations from the structure of the data that is being processed and shines in things like parsing structured documents such as HTML or XML.

CppCon 2022 - Breaking Dependencies - The Visitor Design Pattern in Cpp - Klaus Iglberger

Honorable mentions

NeetCode - 8 Design Patterns EVERY Developer Should Know

How hash maps work

Last lecture one of the homework assignments were to find duplicate files in a folder. One of the best ways to solve this problem is using hashes and hash maps (or dictionaries as they are called in Python). To build on the solution to that problem, below is a video that goes into detail on how hash maps work and how to optimize them:

strager - Faster than Rust and C++: the PERFECT hash table

Building your own tools

The opening gif was created using this script:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

"""Make a nice gif of a few orbits

pip requirements: numpy matplotlib pyorb
Also requires pillow for saving

"""

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import pyorb

plt.style.use("dark_background")
HERE = Path(__file__).parent

# For reproducibility
np.random.seed(8237487)
num = 100

orbs = pyorb.Orbit(
    M0=pyorb.M_sol,
    num=num,
    a=(1 + np.random.randn(num)*0.2) * pyorb.AU,
    e=np.random.rand(num) * 0.1,
    i=np.random.randn(num) * 40.0,
    omega=np.random.randn(num) * 360.0,
    Omega=np.random.randn(num) * 360.0,
    anom=np.random.randn(num) * 360.0,
    degrees=True,
)

days = 3600 * 24
time_step = 1.0 * days
trail = np.linspace(0, -20 * days, 20, dtype=np.float64)
times = np.arange(60) * time_step

# container for fast calculation
orb_trail = pyorb.Orbit(
    M0=pyorb.M_sol,
    num=len(trail),
    direct_update=False,
    auto_update=False,
    degrees=True,
)


def orb_trail_position(orb):
    orb_trail._kep[:, :] = orb._kep[:, 0][:, None]
    orb_trail.mean_anomaly += orb_trail.mean_motion * trail
    orb_trail.calculate_cartesian()
    r = orb_trail.cartesian[:3, :]
    return r


fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(111, projection="3d")
lines = []
dots = []

# Init lines to draw
for orb in orbs:
    r = orb_trail_position(orb)

    (ln,) = ax.plot(
        r[0, :], r[1, :], r[2, :],
        ls="-", c="#68b4ff", alpha=0.5,
    )
    (dot,) = ax.plot(
        r[0, 0], r[1, 0], r[2, 0],
        ls="none", marker=".", c="#00ff2f", alpha=0.5,
    )
    lines.append(ln)
    dots.append(dot)

# insert some rulers
theta = np.linspace(0, np.pi * 2, 200)
for au_dist in [1, 2, 3, 5]:
    ax.plot(
        au_dist * pyorb.AU * np.cos(theta),
        au_dist * pyorb.AU * np.sin(theta),
        ls="-",
        alpha=0.4,
        c="#ffffff",
    )

ax.grid(False)
ax.axis("off")
ax.view_init(14, -6)

box_size = 1
ax.set_xlim([-box_size * pyorb.AU, box_size * pyorb.AU])
ax.set_ylim([-box_size * pyorb.AU, box_size * pyorb.AU])
ax.set_zlim([-box_size * pyorb.AU, box_size * pyorb.AU])


def animate(ind):
    for line, dot, orb in zip(lines, dots, orbs):
        orb.propagate(times[ind])
        r = orb_trail_position(orb)

        line.set_xdata(r[0, :])
        line.set_ydata(r[1, :])
        line.set_3d_properties(r[2, :])

        dot.set_xdata([r[0, 0]])
        dot.set_ydata([r[1, 0]])
        dot.set_3d_properties([r[2, 0]])
    return lines + dots


ani = animation.FuncAnimation(
    fig,
    animate,
    repeat=True,
    frames=len(times),
    interval=50,
)

writer = animation.PillowWriter(
    fps=15,
    metadata=dict(artist="Daniel Kastinen"),
    bitrate=1800,
)
ani.save(HERE.parents[1] / "static" / "media" / "orbits.gif", writer=writer)

Since I work a lot with orbits - I made a tool that specifically converts orbits from Cartesian coordinates to Kepler elements in an efficient way with lots of helper functions to make working with that type of data easy. Since I divided it up into a package I can depend on that package and the functions in this package are heavily tested so I can thrust the results.

Useful concepts

There are many useful concepts to know about when solving small individual problems, below I list a a few common ones:

• Sidecar files
• More to be added here

Homework assignments

Logistic map

The task is to implement a function that can repeatedly apply the logistic map for a specified number of iterations, with a given starting value and a given value of r. The function should return all the values for each step of the repeated application of the map.

Using this implementation, create a large set of starting values in terms of random x0 between [0, 1] and a regular grid of r between [0, 4]. Iterate these values 50 times trough the logistic map.

The goal is to re-produce a arbitrary frame (or frame 50 in this case) of this animation:

You should see Bifurcation occurring, showing that the dynamical system created by repeated applications of the Logistic Map is in fact chaotic for certain values of the functions parameters. The path to chaos is in this case a infinite production of stable regions, also called a period-doubling bifurcation cascade.

Recommended watching and reading

Ideally, during this course, it would be good if you watched the videos and read the articles linked above.