Arrays

In the previous lesson the environment was pinned, and the notebook now gives the same number on two machines. The calculation itself did not change, though: 3,780 cells were walked one by one with nested loops, each cell was held as a separate object, and for every question the data was read again from start to finish. This lesson does the same job with an array interface.

The lesson does not merely show a convenience; it sets up the unit of measure of the course. From here on, two numbers will stand side by side in every lesson: what the abstraction gains and what the abstraction hides. The gain is measured here with three units — how many lines of code, how many full passes over the data, how many bytes. The hidden default is the decision the array makes for you without asking. No library is called; the array behavior is modeled with the standard library, because the aim is not to imitate a tool but to make its behavior measurable.

• NC3. An array interface fixes three things together: contiguous memory, a single dtype and a shape. The gain is a result of this triple and shows up in bytes and in the number of lines at the call site; it does not show up in the number of passes over the data.
• NC4. An array has no empty value. A missing cell has to turn into a number, and which number it turns into is a decision. This decision moves some summaries and does not move others at all.

The Same Job, Two Structures

# arrays.py — MODEL: the three properties of an array interface are written with the standard library:
# contiguous memory, a single dtype, shape. No library is called; the aim is not to imitate
# a tool but to make its behavior measurable. The dataset is the fictional one prepared by the Data Collection and Preparation course.
import inspect
import sys
from array import array

def rng(seed):
    s = seed & 0xFFFFFFFF
    def next_value():
        nonlocal s
        s = (1664525 * s + 1013904223) & 0xFFFFFFFF
        return s / 4294967296
    return next_value

SUBSCRIBERS, PERIODS = 1260, 3
r = rng(20260218)
RAW = [round(4 + r() * 56, 2) for _ in range(SUBSCRIBERS * PERIODS)]
r2 = rng(9110723)
p = [r2() for _ in range(SUBSCRIBERS * PERIODS)]
EMPTY = set(sorted(range(len(p)), key=lambda j: p[j])[:581])
READINGS = [None if j in EMPTY else RAW[j] for j in range(SUBSCRIBERS * PERIODS)]
TABLE = [READINGS[i*PERIODS:(i+1)*PERIODS] for i in range(SUBSCRIBERS)]

PASSES = 0  # one full scan over the data from start to finish

def by_hand(t):
    global PASSES
    subscriber_totals = []
    for s in t:
        a = 0.0
        for x in s:
            if x is not None:
                a += x
        subscriber_totals.append(a)
    overall, n = 0.0, 0
    for s in t:
        for x in s:
            if x is not None:
                overall += x
                n += 1
    above = 0
    for s in t:
        for x in s:
            if x is not None and x > 40.0:
                above += 1
    PASSES += 3
    return subscriber_totals, overall / n, above

class Array:
    def __init__(self, data, shape, dtype="d"):
        self.d, self.shape, self.dtype = array(dtype, data), shape, dtype
    @property
    def size(self):
        return self.shape[0] * self.shape[1]
    def axis_sum(self, axis):
        global PASSES
        PASSES += 1
        a, b = self.shape
        if axis == 1:
            return [sum(self.d[i*b:(i+1)*b]) for i in range(a)]
        return [sum(self.d[i*b+j] for i in range(a)) for j in range(b)]
    def total(self):
        global PASSES
        PASSES += 1
        return sum(self.d)
    def above(self, threshold):
        global PASSES
        PASSES += 1
        return sum(1 for x in self.d if x > threshold)

def with_array(a):
    return a.axis_sum(1), a.total() / a.size, a.above(40.0)

def list_bytes(t):
    b = sys.getsizeof(t)
    for s in t:
        b += sys.getsizeof(s) + sum(sys.getsizeof(x) for x in s if x is not None)
    return b

FLAT = Array([0.0 if x is None else x for x in READINGS], (SUBSCRIBERS, PERIODS))
PASSES = 0
E = by_hand(TABLE)
pass_hand = PASSES
PASSES = 0
D = with_array(FLAT)
pass_arr = PASSES
lines = lambda f: len(inspect.getsource(f).strip().splitlines())
bytes_list, bytes_arr = list_bytes(TABLE), sys.getsizeof(FLAT.d)

row_fmt = lambda w: lambda h: "".join(str(v).ljust(w[0]) if j == 0 else str(v).rjust(w[j])
                                      for j, v in enumerate(h))
y = row_fmt([20, 20, 16, 10])
print(f"shape {FLAT.shape}, element count {FLAT.size}, dtype 64-bit float, "
      f"empty cells {len(EMPTY)}")
print("job: total per subscriber, overall mean, count of cells above 40 m3\n")
print(y(["measure", "hand-written loop", "array interface", "ratio"]))
print(y(["line count", lines(by_hand), lines(with_array), f"{lines(by_hand)/lines(with_array):.1f}x"]))
print(y(["pass count", pass_hand, pass_arr, f"{pass_hand/pass_arr:.1f}x"]))
print(y(["bytes allocated", bytes_list, bytes_arr, f"{bytes_list/bytes_arr:.1f}x"]))
print(y(["bytes per element", f"{bytes_list/FLAT.size:.2f}", f"{bytes_arr/FLAT.size:.2f}", f"{bytes_list/bytes_arr:.1f}x"]))
print(f"\nsame job, two results: mean by hand {E[1]:.6f} m3, with the array interface {D[1]:.6f} m3; "
      f"cells above 40 m3 on both paths {E[2]}")

shape (1260, 3), element count 3780, dtype 64-bit float, empty cells 581
job: total per subscriber, overall mean, count of cells above 40 m3

measure                hand-written loop array interface     ratio
line count                            22               2     11.0x
pass count                             3               3      1.0x
bytes allocated                   188864           30320      6.2x
bytes per element                  49.96            8.02      6.2x

same job, two results: mean by hand 31.895661 m3, with the array interface 26.993180 m3; cells above 40 m3 on both paths 1116

The definition of the three measures stays fixed throughout the course. The line count is the number of lines of code written at the call site, and it is read from the code, not counted by hand. The pass count is the number of full scans made over the data from start to finish; in the model each scan raises a counter by one, so it is a number that is measured, not claimed. The allocated bytes is the space taken up by the structure that holds the data and by the objects inside it. None of the three changes from run to run; time does change, which is why time is not measured. The gain of an abstraction is defined by which of these three numbers it lowers, and its cost by which it does not.

The first row of the table is the expected gain. Written by hand, the same three questions take 22 lines; with the array interface they take 2. The difference is not just brevity of writing: in the hand-written version the missing-cell check is repeated three times, the loop bounds are written three times, and at each repetition there is a chance of leaving one of them out. The abstraction moves these three repetitions to a single place.

The third row is the second gain, and its cause is that the decision about shape and dtype was made once. In nested lists each of the 3,780 numbers is an object of its own; the object itself takes about 24 bytes, the pointer that points to it 8 bytes, and each of the 1,260 small lists that hold the rows takes space as well. The total is 188,864 bytes, 49.96 per element. In the array the 3,780 numbers sit in one contiguous block of memory, exactly 8 bytes per element, and the shape is nothing more than two integers: 30,320 bytes, 8.02 per element. The ratio is 6.2 times.

The second row, however, is not a gain. The pass count is 3 on both paths, because there are three questions and each question reads the data from start to finish once. The array interface lifts the loop out of the call site, not out of the data. This is the most often misread number of the course: because two lines of code are written, a single pass is assumed. Every chained call adds its own pass, and that is exactly what the next lesson measures. In the hand-written version, reducing these three passes to one means merging three loops into one, and writing that takes two lines; in the array interface, making the same merge requires leaving the call chain and going back to the loop. Brevity and pass count stand against each other here.

The last row points to the hidden default. The same job, the same data, two different means: 31.895661 and 26.993180 m³. The 4.90 m³ between them is not a calculation error but a decision that comes from the structure of the array.

Shape and Axis

The Array class in the model carries two fields, and keeping the two apart is the most defining property of an array interface. The d field holds the 3,780 numbers in a single flat block; the shape field is nothing more than two integers and says how that flat block is to be read. The data is not two-dimensional — its reading is.

The direct consequence of this separation is that changing the shape does not copy the data at all. Writing (3, 1260) in place of (1260, 3) moves none of the 30,240 bytes; it only changes the slice bounds in the axis_sum function. This cheapness is the most useful side of the abstraction and at the same time its quietest side: a wrong shape does not produce an error, it only produces another number.

Axis is a direction within that shape, and the numbering in the model runs from left to right: axis 0 is the subscriber direction, axis 1 is the period direction. axis_sum(1) adds the three periods of each subscriber and returns 1,260 values; axis_sum(0) adds all the subscribers of each period and returns 3 values. Both carry the name “sum”, both distribute the same 102,034.22 m³, but one describes subscriber behavior and the other period behavior. In the hand-written loop the direction of the sum is read from the loop lines; in the array interface it is read from a single integer, and that integer is often left as the default. In this lesson the axis is written out explicitly; the rule to be kept until the next lesson is this: an aggregation with no written axis is an aggregation whose direction was not chosen.

An Array Has No Empty Value

# missing.py — MODEL: an array has no "empty" value. The 581 empty cells have to turn into a number
# and which number they turn into is a decision. Three decisions, three means.
DECISIONS = [
    ("fill with zero", lambda v: [0.0 if x is None else x for x in v]),
    ("sentinel -1", lambda v: [-1.0 if x is None else x for x in v]),
    ("filled cells only", lambda v: [x for x in v if x is not None]),
]
z = row_fmt([24, 12, 16, 12])
print(z(["581 empty cells", "elements", "mean (m3)", "over 40 m3"]))
for name, f in DECISIONS:
    a = Array(f(READINGS), (len(f(READINGS)), 1))
    print(z([name, a.size, f"{a.total()/a.size:.6f}", a.above(40.0)]))

DTYPES = [("16-bit integer", "h"), ("32-bit integer", "i"),
          ("32-bit float", "f"), ("64-bit float", "d")]
LITERS = [int(round((0.0 if x is None else x) * 1000)) for x in READINGS]
z2 = row_fmt([18, 15, 15, 30])
print("\n" + z2(["dtype", "bytes/element", "3780 elements", "representation range"]))
for name, t in DTYPES:
    b = array(t).itemsize
    span = (f"-{2**(8*b-1)} .. {2**(8*b-1)-1}" if t in ("h", "i")
            else ("7 significant digits" if t == "f" else "15 significant digits"))
    print(z2([name, b, b * FLAT.size, span]))

TOTAL = FLAT.total()
small = array("f", [TOTAL])[0]
print(f"\nlargest reading {max(LITERS)} liters; a 16-bit integer cannot hold this value")
print(f"dataset total {TOTAL:.2f} m3; stored as {small} in the 32-bit float, "
      f"difference {abs(small-TOTAL):.5f} m3 and above the 0.0005 m3 threshold of the Working Environment lesson")

581 empty cells             elements       mean (m3)  over 40 m3
fill with zero                  3780       26.993180        1116
sentinel -1                     3780       26.839476        1116
filled cells only               3199       31.895661        1116

dtype               bytes/element  3780 elements          representation range
16-bit integer                  2           7560               -32768 .. 32767
32-bit integer                  4          15120     -2147483648 .. 2147483647
32-bit float                    4          15120          7 significant digits
64-bit float                    8          30240         15 significant digits

largest reading 60000 liters; a 16-bit integer cannot hold this value
dataset total 102034.22 m3; stored as 102034.21875 in the 32-bit float, difference 0.00125 m3 and above the 0.0005 m3 threshold of the Working Environment lesson

The three rows are three defensible decisions, and they move the mean between 26.839476 and 31.895661 m³. Filling with zero counts a missing reading as “no water was used”, and because it raises the denominator without raising the numerator it lowers the mean by 4.90 m³. The sentinel value is chosen to mark the missing cell, but if excluding the marker is forgotten, 581 m³ is subtracted from the total and the mean falls even further. The filled cells only decision uses 3,199 as the denominator and gives the number the previous course measured. All three have a reason; the only thing without a reason is leaving the decision unwritten.

The fourth column shows where this decision does not show. The count of cells above 40 m³ is 1,116 under all three decisions, because all three fill values lie below the threshold. An audit that looks only at this number cannot tell the three decisions apart. The number the decision moves and the number that makes the decision invisible sit in the same table, and which of the two is reported matters as much as the decision itself.

Notice that the hand-written version never had to make this decision. In nested lists the empty cell stays as None; this value is not a number, it does not enter the sum, and the is not None check leaves it out explicitly inside the loop. The decision is there but it is visible, because it was written one by one in every loop. In the array the decision is never written; a fill value is chosen on the line that builds the array, and that value, never mentioned again, mixes into all the summaries. What the abstraction hides is not the decision itself but where the decision is made.

The lower table measures the dtype choice. The same 3,780 elements take between 7,560 and 30,240 bytes depending on the chosen dtype; the fourfold difference is paid not in memory but in representation range. The largest reading is 60,000 liters and exceeds the 32,767 limit of the 16-bit integer — choosing this dtype is not a fourfold memory gain but a silent overflow. The 32-bit float does not overflow but holds roughly seven significant digits; because the dataset total of 102,034.22 m³ needs eight significant digits, it is stored as 102,034.21875. The difference is 0.00125 m³ and above the 0.0005 m³ threshold of the previous lesson, which means this dtype choice alone is enough to make two runs differ.

Summary

• An array interface fixes contiguous memory, a single dtype and shape together; here the shape is (1260, 3) and the element count is 3,780.
• The same three questions take 22 lines by hand and 2 lines with the array interface; the allocated bytes fall from 188,864 to 30,320, from 49.96 to 8.02 bytes per element.
• The pass count is 3 on both paths: the abstraction lifts the loop out of the call site, not out of the data.
• Because an array has no empty value, 581 empty cells turn into a number; the mean comes out 26.993180 with zero fill, 26.839476 with a sentinel value, and 31.895661 m³ with filled cells only.
• The same three decisions do not change the count of cells above 40 m³ at all (1,116); the dtype choice keeps the 3,780 elements between 7,560 and 30,240 bytes, and a 16-bit integer cannot hold the 60,000-liter reading.

Next Step

The unit of measure is set up, and one of the three numbers came out contrary to expectation: the pass count did not fall. The next lesson takes up exactly this number. The same derivation job can be done in a single loop with one pass, or it can be written as a chain of vectorized calls; the second way adds one pass and one intermediate array with every call. The same lesson also counts the two decisions vectorization gives silently: the promotion of the dtype and overflow in a narrow dtype. In the last table of this lesson it was seen that a 16-bit integer cannot hold the 60,000-liter reading; the next lesson measures what happens when this limit is exceeded.