I wanted to investigate the warmup behavior of the PyPy interpreter, so I wrote a somewhat arbitrary microbenchmark:

all_results = set()
num = int(sys.argv[1])

class A(object):
    pass

def main():
    res = 0
    for i in range(num):
        a = A()
        a.x = i
        d = {"a": a.x}
        l = [0, 1, d["a"]]
        res += l[-1]
    all_results.add(res)

This function is the ideal case for PyPy's JIT compiler, as it has a tight loop with many object allocations that all have predictable lifetimes. There is no control flow inside the loop. The JIT compiler can optimize this loop body extremely well, to just an integer addition.

However, I wanted to investigate the warmup behavior of the JIT compiler. PyPy's tracing JIT will trace this loop at its 1041th iteration. Tracing a loop is rather expensive in PyPy, due its meta-tracing architecture. When the loop is traced, there is actually a stack of two interpreters: the meta-tracing interpreter will execute the Python interpreter, while the latter executes one iteration of the Python loop. The meta-tracing interpreter executes the loop, but also keeps a log of all the operations that were executed by PyPy's Python interpreter. This log – the trace – is then optimized and compiled to machine code. From then on, that chunk of machine code will be executed instead of interpreting the loop body.

What really interests me here is just how costly the double-interpretation overhead of the meta-tracing interpreter is. I had no intuition about the slowdown factor of it, is it 10x, 100x, 1000x, 10,000x, 100,000? My goal is to understand the order of magnitude of that slowdown factor.

All the numbers in this blog post are run on PyPy 2.7.18 (because I had a nicely instrumented version of PyPy2 around), I expect the results to be similar (ie the same order of magnitude at least) on PyPy 3.11. Timing results are from my AMD Ryzen 7 PRO 7840U, running Ubuntu Linux 24.04.2.

Estimating the slowdown factor of the meta-tracing interpreter

To estimate the slowdown factor we can run this benchmark with num set to exactly 1041. Since the JIT will trace the 1041th iteration, this is the maximally bad case for the JIT. The JIT will trace the loop, and the meta-tracing interpreter will execute the Python interpreter for one iteration of the loop, which is exactly what we want to measure. However, since we don't execute another iteration, we never execute the compiled machine code, so we can measure the overhead of the meta-tracing interpreter. If we set the environment variable PYPYLOG=jit-summary:outputfilename, we can see how much time the meta-tracing interpreter spent tracing the loops.

To make the numbers slightly larger, I actually add and run 100 copies of the above main function, to reduce the noise of the measurement.

from __future__ import print_function
import sys

all_funcs = []
all_results = set()
s = """
class A%s(object):
    pass

def main%s():
    res = 0
    for i in range(num):
        a = A%s()
        a.x = i
        d = {"a": a.x}
        l = [0, 1, d["a"]]
        res += l[-1]
    all_results.add(res)
all_funcs.append(main%s)
"""
fullstr = "\n".join(s % (i, i, i, i) for i in range(100))
exec(fullstr)

import time

num = int(sys.argv[1])

t1 = time.time()
for fn in all_funcs:
    fn()
t2 = time.time()
print(num, t2 - t1)

Now we can run this benchmark 100 times with bash to get a good average across many process execution:

for i in {0..100..1};
    do PYPYLOG=jit-summary:out%d pypy -S x.py 1041 >> times;
done

Afterwards, we can look at the out%d files to see how much time the meta-tracing interpreter spent tracing the loops (%d is expanded to the process id by the PyPy logging system). The content of these files will look roughly like this:

[c8ca76f3db881] {jit-summary
Tracing:        100     0.033022
Backend:        100     0.009439
TOTAL:                  0.072409
ops:                    66200
heapcached ops:         20200
recorded ops:           11200
  calls:                2000
guards:                 2200
opt ops:                9700
opt guards:             1800
opt guards shared:      1300
forcings:               0
abort: trace too long:  0
abort: compiling:       0
abort: vable escape:    0
abort: bad loop:        0
abort: force quasi-immut:       0
abort: segmenting trace:        0
virtualizables forced:  0
nvirtuals:              2500
nvholes:                0
nvreused:               1100
vecopt tried:           0
vecopt success:         0
Total # of loops:       100
Total # of bridges:     0
Freed # of loops:       0
Freed # of bridges:     0
[c8ca76f3e9618] jit-summary}

Let's extract the tracing times with grep:

grep Tracing out* | grep -o "0[.].*" > tracingtimes

Now we can compute the slowdown factor like this:

>>> with open('tracingtimes') as f:
...     f = [float(l) for l in f]
...     tracingtime = sum(f) / len(f)
...     
>>> with open('times') as f:
...     f = [float(l.split()[1]) for l in f]
...     exectime = sum(f) / len(f)
...     
>>> regulartime = (exectime - tracingtime) / 1040
>>> print(regulartime)
3.550879183375572e-05
>>> print(tracingtime)
0.03153696039603962
>>> slowdown = tracingtime / regulartime
>>> print(slowdown)
888.1451259645404

(random side-note: I want to complain about the fact that the pygments lexer does not support the >>>> syntax of PyPy's interactive interpreter.)

This means that the meta-tracing interpreter is about ~900x slower than the regular interpreter for this benchmark. (Note that the tracing time alone is not the full warmup time of the JIT, since it doesn't include the time spent optimizing the trace and producing machine code.)

How many instructions does the meta-tracing interpreter execute?

Let's understand the work the Python interpreter does in one iteration of the loop. We can disassemble one of the main function to see what it does:

>>> import dis
>>> dis.dis(main0)
  6           0 LOAD_CONST               1 (0)
              3 STORE_FAST               0 (res)

  7           6 SETUP_LOOP              87 (to 96)
              9 LOAD_GLOBAL              0 (range)
             12 LOAD_GLOBAL              1 (num)
             15 CALL_FUNCTION            1
             18 GET_ITER            
        >>   19 FOR_ITER                73 (to 95)
             22 STORE_FAST               1 (i)

  8          25 LOAD_GLOBAL              2 (A0)
             28 CALL_FUNCTION            0
             31 STORE_FAST               2 (a)

  9          34 LOAD_FAST                1 (i)
             37 LOAD_FAST                2 (a)
             40 STORE_ATTR               3 (x)

 10          43 BUILD_MAP                0
             46 LOAD_FAST                2 (a)
             49 LOAD_ATTR                3 (x)
             52 LOAD_CONST               2 ('a')
             55 STORE_MAP           
             56 STORE_FAST               3 (d)

 11          59 LOAD_CONST               1 (0)
             62 LOAD_CONST               3 (1)
             65 LOAD_FAST                3 (d)
             68 LOAD_CONST               2 ('a')
             71 BINARY_SUBSCR       
             72 BUILD_LIST               3
             75 STORE_FAST               4 (l)

 12          78 LOAD_FAST                0 (res)
             81 LOAD_FAST                4 (l)
             84 LOAD_CONST               4 (-1)
             87 BINARY_SUBSCR       
             88 INPLACE_ADD         
             89 STORE_FAST               0 (res)
             92 JUMP_ABSOLUTE           19
        >>   95 POP_BLOCK           

 13     >>   96 LOAD_GLOBAL              4 (all_results)
             99 LOOKUP_METHOD            5 (add)
            102 LOAD_FAST                0 (res)
            105 CALL_METHOD              1
            108 POP_TOP             
            109 LOAD_CONST               0 (None)
            112 RETURN_VALUE  

The loop body starts at bytecode offset 19 and ends at bytecode offset 95. The loop body consists of 25 Python bytecode instructions.

When the JIT compiler traces this loop, it will execute the Python interpreter for one iteration of the loop on top of the meta-tracing interpreter. I built a special version of PyPy that counts the number of opcodes executed by the meta-tracing interpreter while tracing. Here are the results:

The names of these functions aren't too important, what matters is that the meta-tracing interpreter executes 3675 operations to trace one iteration of the loop in main, which is a lot. It also traces 103 operations, most of which are then removed again by the trace optimizer, leaving a trace with just 22 operations:

label(p0, p1, p6, p7, i45, i32, p20, p31, p35, i50, i36)
debug_merge_point(0, 0, '<code object main0. file '<string>'. line 5> #19 FOR_ITER')
i52 = int_lt(i50, 0)
guard_false(i52)
i53 = int_ge(i50, i36)
guard_false(i53)
i55 = int_add(i50, 1)
debug_merge_point(0, 0, '<code object main0. file '<string>'. line 5> #25 LOAD_GLOBAL')
# update iterater object
setfield_gc(p20, i55, descr=<FieldS pypy.objspace.std.iterobject.W_AbstractSeqIterObject.inst_index 8>)
guard_not_invalidated()
debug_merge_point(0, 0, '<code object main0. file '<string>'. line 5> #88 INPLACE_ADD')
# perform the addition, checking for machine integer overflow
i56 = int_add_ovf(i45, i50)
guard_no_overflow(descr=<Guard0x729830a52170>) [p0, p6, p7, p20, p1, i50, i32, i45]
debug_merge_point(0, 0, '<code object main0. file '<string>'. line 5> #92 JUMP_ABSOLUTE')
# check whether ctrl-c was pressed
i58 = getfield_raw_i(125998063394752, descr=<FieldS pypysig_long_struct.c_value 0>)
i60 = int_lt(i58, 0)
guard_false(i60)
debug_merge_point(0, 0, '<code object main0. file '<string>'. line 5> #19 FOR_ITER')
jump(p0, p1, p6, p7, i56, i50, p20, p31, p35, i55, i36)

Note how in particular, no allocations are left in this trace. The optimize recognized that neither the A instance, the dictionary, nor the list need to be allocated.

How long does it take for the JIT compilation time to be amortized?

To find out after how many iterations the JIT compilation time is amortized, we can run the benchmark with varying num values and measure the execution time, both with and without the JIT. We do this by starting many processes using more bash for loops. We can then plot the execution time against the number of iterations to see when the time with the JIT becomes lower than the time without the JIT.

for i in {0..20000..1}; do pypy -S x.py "$i" >> datajit ; done; for i in {0..20000..1}; do pypy --jit off -S x.py "$i" >> datanojit ; done

And we can also run the same script on CPython, while we are at it:

for i in {1..10000..1}; do python3.13 -S x.py "$i" >> datacpy ; done;

(this is a bit apples against oranges because I am comparing PyPy2 against CPython 3.13 but oh well.)

Both will run for a while and then we can plot the data.

# Co-created by Copilot
import matplotlib.pyplot as plt
import numpy as np

THRESHOLD = 1041  # JIT threshold

def plot_trace_costs(maximum, save_as, baseline_filename='datanojit', baseline_legend='PyPy No JIT'):
    with open('datajit', 'r') as f:
        datajit = [tuple(map(float, line.strip().split())) for line in f.readlines()]
    with open(baseline_filename, 'r') as f:
        datanojit = [tuple(map(float, line.strip().split())) for line in f.readlines()]
    # scale the y values by 1000
    datajit = [(x, y * 1000) for x, y in datajit]
    datanojit = [(x, y * 1000) for x, y in datanojit]
    origdatajit = datajit.copy()
    # sort the data by x value
    datajit.sort(key=lambda x: x[0])
    datanojit.sort(key=lambda x: x[0])
    # filter out lines with x value greater than maximum
    datajit = [x for x in datajit if x[0] <= maximum]
    datanojit = [x for x in datanojit if x[0] <= maximum]
    # create new figure
    plt.figure(figsize=(10, 6))
    plt.clf()

    # unpack data
    jit_x, jit_y = zip(*datajit) if datajit else ([], [])
    nojit_x, nojit_y = zip(*datanojit) if datanojit else ([], [])

    # plot the data points, using a scatter plot. use circle size of 2
    plt.scatter(jit_x, jit_y, label='PyPy w/ JIT', alpha=0.3, s=1)
    plt.scatter(nojit_x, nojit_y, label=baseline_legend, alpha=0.3, s=1)
    # add thin vertical line at x=THRESHOLD
    plt.axvline(x=THRESHOLD, color='gray', linestyle='--', linewidth=1)
    # add text to the left of the line, near the top saying "JIT threshold"
    all_y = list(jit_y) + list(nojit_y)
    if all_y:
        plt.text(THRESHOLD, max(all_y) * 0.9, 'JIT threshold', color='gray', fontsize=10, ha='left')

    # fit a line to the JIT data, starting from x=THRESHOLD
    jit_x_fit = [x for x, _ in origdatajit if x >= THRESHOLD]
    jit_y_fit = [y for x, y in origdatajit if x >= THRESHOLD]
    assert len(jit_x_fit) == len(jit_y_fit), "JIT x and y data lengths do not match"
    assert jit_x_fit and jit_y_fit
    jit_fit = np.polyfit(jit_x_fit, jit_y_fit, 1)
    jit_fit_line = np.polyval(jit_fit, jit_x)
    plt.plot(jit_x, jit_fit_line, color='black', linestyle='--', linewidth=1)
    assert nojit_x and nojit_y
    nojit_x_fit = [x for x in nojit_x if x >= THRESHOLD]
    nojit_y_fit = [y for x, y in zip(nojit_x, nojit_y) if x >= THRESHOLD]
    assert nojit_x_fit and nojit_y_fit
    nojit_fit = np.polyfit(nojit_x_fit, nojit_y_fit, 1)
    nojit_fit_line = np.polyval(nojit_fit, nojit_x)
    plt.plot(nojit_x, nojit_fit_line, color='black', linestyle='--', linewidth=1)

    # compute intersection point of the two lines
    intersection_x = (nojit_fit[1] - jit_fit[1]) / (jit_fit[0] - nojit_fit[0])
    intersection_y = jit_fit[0] * intersection_x + jit_fit[1]
    # also compute the speedup of the JIT and No JIT fits by dividing the slopes
    jit_speedup = 1 / (jit_fit[0] / nojit_fit[0])

    # add text box to the upper right corner with the fit parameters
    fit_text = f'PyPy w/ JIT fit: y = {jit_fit[0]:.6f}x + {jit_fit[1]:.6f}\n{baseline_legend} fit: y = {nojit_fit[0]:.6f}x + {nojit_fit[1]:.6f}'
    speedup_text = f'Speedup: {jit_speedup:.2f}x\nIntersection: ({intersection_x:.2f}, {intersection_y:.2f})'
    fit_text += "\n" + speedup_text
    plt.text(0.95, 0.95, fit_text, transform=plt.gca().transAxes, fontsize=10,
             verticalalignment='top', horizontalalignment='right',
             bbox=dict(facecolor='white', alpha=0.5, edgecolor='none'))
    # compute the y-values of both fits at the JIT threshold
    jit_y_at_threshold = jit_fit[0] * THRESHOLD + jit_fit[1]
    nojit_y_at_threshold = nojit_fit[0] * THRESHOLD + nojit_fit[1]
    diff = jit_y_at_threshold - nojit_y_at_threshold
    diff_text = f'Difference at threshold: {diff:.2f} ms'
    print(diff_text)
    plt.xlabel('Iterations')
    plt.ylabel('time taken (microseconds)')
    plt.grid()
    plt.title(f'Warmup behaviour: PyPy w/ JIT vs {baseline_legend}')
    # legend on the top left corner
    plt.legend(loc='upper left')
    plt.savefig(save_as)

if __name__ == '__main__':
    plot_trace_costs(maximum=20000, save_as='trace_costs_start.svg')
    plot_trace_costs(maximum=2000, save_as='trace_costs_very_start.svg')
    plot_trace_costs(maximum=10000, save_as='trace_costs_cpy.svg', baseline_filename='datacpy', baseline_legend='CPython 3.13')

This is what the plot looks like:

The plot shows that the JIT compilation time is amortized after about 2600 iterations. The JIT compilation time is about 39ms (which consists of the meta-tracing time plus code generation time). If the program runs for a very long time, the speedup of the JIT will approach >200x (but the example was picked to make this easy for the JIT, so it's not indicative of the performance of realistic Python code).

We can zoom in the at the first 2000 data points:

And here's the comparison with CPython 3.13:

Compared to CPython the speedup is ~160x and PyPy catches up with CPython only at around 4800 iterations. I don't know why CPython's performance is so much more noisy across the 10,000 iterations.

Conclusion

After these quick measurements we conclude that the meta-tracing interpreter is three orders of magnitude slower than the regular interpreter. We (Yusuke Izawa and I) plan to improve the performance of the meta-tracing interpreter in the future, and I look forward to seeing how the numbers change then.

Acknowledgements

Thanks to Yusuke Izawa for his work on instrumenting the meta-tracing interpreter, and for our previous collaboration.

I mainly post on the PyPy blog, but there are sometimes things I write that don't really fit there. In 2018 and 2020 I wrote two essays about (computational) art which lived on glitch.com. Now that that is going to shut down, I decided to migrate and backdate them here, together with some shorter backdated draft posts that never really made it anywhere (or only to social media).

We don't do this regularly or anything, but every couple of years I look at how the PyPy binary sizes have developed in the meantime. Looks like an ok balance between occasionally cleaning something up and thereby shrinking the binary; and then slow growth or explicit time/binary-size-tradeoffs.

This is the size of PyPy 2.7, btw. I can't use the 3.x variant, because the size of that changes due to newer Python versions implemented. Interestingly enough PyPy 3.x is mostly smaller than 2.7. PyPy 3.11 is 58.8 MiB, despite having a bunch more features. I never quite figure out why. Just for comparison, CPython 3.11 is 25.2 MiB.

Here's a plot of CPython/PyPy 3.x binary sizes:

I really enjoyed the two posts by Christophe Grand about Datalog, Writing the Worst Datalog Ever in 26loc and Half Dumb Datalog in 30 loc (even more than the first one), so I ended up porting them to Python.

There's a bunch of differences to the Clojure versions:

• they are roughly twice as long
• use generators and thus have to environment threading somewhat differently
• a custom class for variables (no symbols in Python)
• some dictionary copying is going on because of mutability (and I didn't want to think too much about when mutation is happening).

The semi-naive implementation ended up being roughly 50 lines, which is pretty reasonable I thought.

The code for them both is in this Gist:

I was inspired by Phil Zucker's blog post Naive Automata Minimization. On Twitter he wrote:

Maybe a way to actualize this is using knot tied rational trees instead of integer state ids. But this has its own confusions.

I really wanted to try to implement this in SWI-Prolog.

?- X = f(X).
X = f(X).

?- X = f(X), X = f(f(f(f(f(f(f(A))))))).
X = A, A = f(A).

>>> l = []
>>> l.append(l)
>>> l
[[...]]

>>> import copy
>>> copy.deepcopy(l)
[[...]]
>>> copy.deepcopy(l) == l
Traceback (most recent call last):
  File "<python-input-9>", line 1, in <module>
    copy.deepcopy(l) == l
RecursionError: maximum recursion depth exceeded

dfa_a_plus(S1) :-
    S1 = state(false, S2, none),
    S2 = state(true, S2, none).

dfa_a_plus_inefficient(S1) :- % less efficient example automaton for a+ with extra states
    S1 = state(false, S2, none),
    S2 = state(true, S3, none),
    S3 = state(true, S4, none),
    S4 = state(true, S2, none).

?- dfa_a_plus(S1), dfa_a_plus_inefficient(S2), S1=S2.
 S1 = S2, S2 = state(false, _S1, none), % where
     _S1 = state(true, state(true, state(true, _S1, none), none), none).

Matching

So far we only defined automata, but we can also write a predicate that matches one against a list of a or b.

match([], state(true, _, _)).
match([a | Rest], state(_, StateA, _)) :- match(Rest, StateA).
match([b | Rest], state(_, _, StateB)) :- match(Rest, StateB).

Some example calls:

?- dfa_a_plus(S), match([a, a, a, a], S). % succeeds

?- dfa_a_plus(S), match([], S).
false.

?- dfa_a_plus(S), match([a, a, a, b], S).
false.

?- dfa_a_plus(S), match([b], S).
false.

Let's look at a more complicated DFA example, the one from the blog post:

dfa_blog_post(S1) :-
    S1 = state(false, S2, S3),
    S2 = state(false, S4, S3), % states S2 and S3 are equivalent
    S3 = state(false, S5, S3),
    S4 = state(true, S5, S4), % states S4 and S5 are equivalent
    S5 = state(true, S4, S4).

(there's a nice diagram in Phil's post.)

If we want to minimize it, we can assign a number to each unique state reachable from the starting state. We can use an association list, i.e. a list where the entries are terms State/Number. States that are observationally equivalent unify, so we can just use Prolog unification (via the "member" predicate which checks whether a term is an element of a list).

number_states(D, L) :-
    % number the states, the first number is 0
    number_states_helper(D, 0, _, [], LRev),
    % reverse the result, because it is built in reverse order
    reverse(LRev, L).

number_states_helper(none, Num, Num, L, L). % none doesn't need a number
number_states_helper(State, Num, Num, L, L) :-
    % if we already gave a number to the state, we are done.
    member(State/_, L).
number_states_helper(State, NumIn, NumOut, LIn, LOut) :-
    % if we've never seen the state, we assign it NumIn
    not(member(State/_, LIn)),
    L1 = [State/NumIn | LIn],
    State = state(_, NextStateA, NextStateB),
    succ(NumIn, Num1),
    number_states_helper(NextStateA, Num1, Num2, L1, L2),
    number_states_helper(NextStateB, Num2, NumOut, L2, LOut).

We can apply this predicate from the non-minimal DFA:

:- dfa_blog_post(D), number_states(D, L), length(L, Length), write(Length), nl, write(L), nl.
3
@([state(false,state(false,S_1,S_2),S_2)/0,state(false,S_1,S_2)/1,S_1/2],[S_1=state(true,state(true,S_1,S_1),S_1),S_2=state(false,state(true,S_1,S_1),S_2)])

The list L of numbered states is kind of hard to grasp due to the messy way that rational tree terms gets printed, but the important part is that we found out that three states is the minimal size.

After we numbered the states we can use that numbering to produce a graphviz graph for the DFA:

to_dot(D) :-
    number_states(D, L),
    write("digraph G {"), nl,
    write('start [label="", shape=none]'), nl,
    write("start -> 0"), nl,
    to_dot_list(L, L),
    write("}"), nl.

to_dot_list([], _).
to_dot_list([State/NumState | Rest], FullList) :-
    State = state(Accept, StateA, StateB),
    to_dot_node(Accept, NumState),
    to_dot_edge(StateA, NumState, FullList, a),
    to_dot_edge(StateB, NumState, FullList, b),
    to_dot_list(Rest, FullList).

to_dot_node(true, Num) :-
    write(Num), write(" [shape=doublecircle]"), nl.
to_dot_node(false, Num) :-
    write(Num), write(" [shape=circle]"), nl.

to_dot_edge(none, _, _, _).
to_dot_edge(State, NumPrev, FullList, Char) :-
    member(State/NumState, FullList),
    write(NumPrev), write(" -> "), write(NumState), write(" [label="), write(Char), write("]"), nl.

If we apply this to the example we get:

digraph G {
start [label="", shape=none]
start -> 0
0 [shape=circle]
0 -> 1 [label=a]
0 -> 1 [label=b]
1 [shape=circle]
1 -> 2 [label=a]
1 -> 1 [label=b]
2 [shape=doublecircle]
2 -> 2 [label=a]
2 -> 2 [label=b]
}

Which looks like this:

Somebody was Wrong on the Internet about the performance of Python, using computing Collatz sequence as a benchmark. Therefore I had to try to measure how long this takes on PyPy, with and without adding a cache for already seen numbers. Takes about 4min on my laptop for 1 billion numbers, or 30s with a cache. Code is here:

For a lesson about JITs in a colleague's compiler course I wrote a small just-in-time compiler for deterministic finite automata in Python, using PeachPy. I planned to write a complete blog post about them, but it doesn't seem to quite be happening, so below is the code. It has a few comments at least.

The idea is to show a few of the typical things a JIT can do that a static compiler can't: - switching between interpreter and compiled code - profiling to find out which part of the program are executed most commonly - deoptimization back to the interpreter - patching of existing code if new code is added

Here's the code:

(if you don't want to read anything just click this button to play around: TLDR)

Last year I had the chance to experience the choreographer William Forsythe's video installation City of Abstracts in the Folkwang museum in Essen. Here's a video of a few people in front of it:

The installation consists of an open space, a camera and a big screen. The camera films the people that are moving and standing in the space, the screen shows a distorted version of the filmed images. The distortion seems to have the following properties: people that are standing still appear only minimally distorted on the screen (for example the person on the left in the back, in the first part of the video). People that are moving across the field of view of the camera turn into diagonally stretched versions of themselves (for example the person moving from right to left at around 0:20). The head is ahead (hrm) in their direction of movement. If they then stand still, their heads stand still first, while the rest of their bodies slowly catch up, from top to bottom.

While I was standing in the museum playing around, I started to wonder how the installation worked. At first I thought it must be a fairly complicated effect, that somehow analyzed motion and did a complicated transformation. However, after some more playing, running back and forth etc. I became convinced that the program is actually quite simple. In the following I want to explain how it works and recreate the effect in Javascript.

A Miniature City

Because it is quite hard to reason about what happens in a big video, let's start by looking at a miniature version of the problem. The following is a 5x5 pixel input video where a stack of pixels first moves left, then right, and then back and forth for a bit:

Start

We want to find out how this input video has to be transformed into an output video in a way that recreates the effect of Forsythe's installation.

Above we saw that the bottom of a moving object is somehow further back in time than the top. This could mean that the output video uses older and older lines of pixels, from top to bottom. The top row of pixels would be from the current camera input, the second row from one frame back, the third row from two frames back, and so on. Let's try this for the miniature input video:

Start Pause

This seems to recreate the effect! When the stack of pixels moves left, it is transformed into a bottom left to top right diagonal, when it moves back it slowly becomes a top left to bottom right diagonal. When it stands still, the top pixels stand still first in the output, the lower pixels slowly catch up. Below you can look at the output pixel video:

Start

The Real Thing

Now that we know how the effect works, we can apply it to a real video, using a webcam. If you click Start and give permission, you can try it yourself (The webcam video is not transmitted anywhere, just shown below). The left video is the webcam output, the right video is the transformed version. If you are interested in the Javascript code, you can use the Glitch button in the top right corner to inspect (and edit) it.

The first video is directly the output of the webcam, the second video shows the "timeshifted" output as explained above. To see anything interesting happening, you need to move around in front of the webcam. Sideways motions are easiest to understand. Remember that the second video takes a while to catch up with your movements. Some fun things to try: walking across the room, slowly shaking your head, jumping, rotating your arms...

Start

Fullscreen

(Aside: I am quite impressed that it works just fine to implement this piece of video processing in Javascript and have it run relatively fluidly in a browser on a phone.)

One parameter we can play with in the above transformation is we can make the lower pixels catch up faster with the upper pixels by doing the transformation block-wise, instead of a single row of pixels at a time. If you increase the block size below, the output video will become less laggy, but on the other hand more blocky.

Relationship to Rolling Shutter

An amusing sidenote is that we can use the transformation to understand the rolling shutter effect. Rolling shutter happens when taking a picture of a fast moving object with a smartphone camera. Those cameras don't record the whole image at once, but do it line by line. And if the object moves very fast relative to the time the phone takes to record the whole image, different parts of the final image show the object at different points in time (and thus space). What the effect above does is artificially slow down the creation of a picture. Therefore we can recreate the effects that rolling shutter gives with objects that move much more slowly, for example by rotating an arm.

LIT Fischer developed his light-graphics step by step, an oscillating rhythmically-flowing, undulating movements on the screen – "heterophonics" – as he calls them. The variety of the undulating lightwaves, often brought about by minute deviations from the original program, fascinate the artist. With his light-movement studies he actually succeeds in visualizing series of musical tones, sound sequences and meditative song. Dr. Eva Karcher

I've been fascinated by the light art installations of artist Jürgen LIT Fischer since a while, e.g. "Fraktal" as part of the Tetraeder in Bottrop:

Photo CC-BY Mirek Claßen

Later I found out that LIT Fischer used to do computer art in the 1980s using Fortran and a plotter. Here are two examples:

"Mit Zwischenraum, hindurchzuschaun" ("with space between to look through"), Jürgen LIT Fischer, 1986. Aside: the title is a part of a poem by Christian Morgenstern, "Der Lattenzaun" (click through for an English translation)

"Intervalle" ("intervals"), Jürgen LIT Fischer, 1987

Dreiklang

"Dreiklang" ("triad"), Jürgen LIT Fischer, 1987

"Dreiklang", top three panels, 1987

Looking through some of the books about his exhibitions I found that one of them contained fragments of the Fortran code that he used as the basis of one of the installations:

Despite not knowing Fortran I decided to still try to port the code to Javascript (which I of course also don't really know, but at least it's easier to Google and more importantly, to run). I guessed that the Z* functions were responsible for communicating with the plotter. In particular that the ZDRAW function had to move the plotter pen somehow.

After some amount of fiddling I managed to produce the output below. I used different parameters than in the Fortran source, but reading some more about his working style, Fischer anyway played around with the parameters a lot (in fact, his final output is usually a layer of many different semi-transparent versions of the outputs of the same program, using different colors).

What we see is basically a series of sine waves, plotted at a fixed distance below each other. They all have the same frequency but the amplitude goes up from 0 to some value, and then down to 0 again. What I find cool is how simple the ingredients are: just lots of sine waves, but the result is still quite an interesting effect.

Compared to the Fortran code I changed things somewhat: I fixed the code duplication between the code that computes the upper and the lower half of the waves and I removed the extra level of loops responsible for increasing the stroke width by drawing the same line several times with a tiny offset (that's much easier in svg!). I also tried to translate the variable names from abbreviated German to English. If you want to read or modify the source code, use the Glitch button in the upper right corner!

The Discontinuity in the Middle

In the picture of the real artwork above, those waves are split up into three acrylic glass panels. In the top middle panel there is what looks like some strange discontinuity in the middle.

"Dreiklang", top middle panel, rotated

I suppose that was done intentionally, but I still wanted to try to figure out where it came from. And indeed, playing with the parameters some more, we can reproduce it:

To understand this effect, let's look at an even more exaggerated view:

I have added in red to the left the size of the amplitude of each sine wave. On the peak of every wave, the lines above bunch together, and the lines below are spread apart. The eye notices the sharp jump between the distances of the lines in the middle. The reason for that jump is the fact that the amplitudes of the stacked waves vary linearly up, and then down. But since that interpolation function is piecewise linear, it has a non-smooth part in the middle.

This is going beyond understanding the original code, but we can fix that! Instead of interpolating the amplitudes linearly, we can do it smoothly, e.g. again with a sine wave:

Interactive Version

If you want to play around with the parameter space too, here's an interactive version:

right-click and save as to download svg