A dataflow language for deep weather forecasting: Riverflow | Training and preprocessing

RiverFlow data loader, part 2 of 2. Previous: A small language for loading messy data

At the end of part 1,every variable was a clean NumPy array: a column of rainfall values, a column of dates, a 1839 × 24 matrix of hourly river levels. A forecasting model needs something else. It needs samples, and each sample combines several variables over a stretch of time: the river over the last month, the rain over the last week and the next few days, and as target, the river over the coming week.

Building those samples by hand is where off-by-one errors live. The sources also start on different days. The Lozzolo weather station begins on 1 October 2001, the river series on 23 March 2018. Every window has to be measured from the same day in every file, and every sample has to have all its windows complete. The .make section describes this in a few lines, and its interpreter logs enough numbers that you can check the result by hand. This post does that check for three examples from the repository [1].

A plan describes one sample around a moving day

Here is the plan from PlotData.makefile, with its commented-out lines removed:

.make

plan_file NewFile = EXAMPLESROOT + '/River Height/savefileN'
plan_file TestFile = EXAMPLESROOT + '/River Height/testN'
log_file Logs = EXAMPLESROOT + '/River Height/savefileN.txt'

begin plan NewPlan expecting attempt_recovery
{
    align NormalizedRes, NormalizedResT against RData, RData as date with format %Y-%M-%D
    align PNove against Data as date with format %Y-%M-%D
    consider x

    take 30 before x from NormalizedRes
    take 7 before x and 6 after x from PNove

    make NormalizedResT the target and take 6 after y from NormalizedResT
    pair x and target
}
end plan

split NewPlan into NewFile, TestFile as 75, 25
compile NewPlan
log NewPlan into Logs

Read it as a description of one sample. align says which column gives each variable its position in time: the river variables are keyed by RData, the rainfall PNove (Italian "precipitazione dalle 9 alle 9", rain from 9:00 to 9:00) by Data. consider x introduces the moving day. Each take line then asks for a window relative to that day. make ... the target marks the variables that go into y instead of x. NormalizedResT is a copy of the normalized river series made in .act, because a variable can only have one window per plan.

The last three lines act on the finished plan. split assigns output files and proportions, compile builds the arrays and saves them, and log writes a description of the result. pair x and target reads well but does nothing: the parser only checks that consider x came before it. The pairing happens in any case.

The window statements are parsed by splitting on keywords:

if 'and' in statement:
    request[0] = int(statement.split('take')[1].split('before')[0].strip())
    request[1] = int(statement.split('and')[1].split('after')[0].strip())
elif 'before' in statement:
    request[0] = int(statement.split('take')[1].split('before')[0].strip())
elif 'after' in statement:
    request[1] = int(statement.split('take')[1].split('after')[0].strip())
else:
    self.aligner.singleton(var_name)
    return

self.aligner.wants(var_name, request)

So a window is two integers, [before, after]. take x from V with no count is [0, 0], one value per sample.

Alignment finds a shared start, then counts rows

When the parser reaches consider x, it builds an Aligner from all the align statements. The first job is to find where the sources overlap:

align_data = [self.var_vec.get_variable(var) for var in self.alignment]
intersection = align_data[0]

independent_vars = [self.var_vec.get_variable(var) for var in set(self.alignment)]
for data in independent_vars:
    intersection = intersection[np.in1d(intersection, data)]
...
bottom_aligns = [np.nonzero(data1 == intersection[0])[0][0] for data1 in align_data]
self.init_align = max(bottom_aligns)

The intersection holds the keys that appear in every alignment column. The aligner takes the first shared key, 23/03/2018 in the river example, and records at which row it appears in each source. From then on it works with row numbers only. Each variable gets a window, the range of rows inside the overlap, and a budget, the number of rows it has outside the overlap on each side.

The budget is what makes the counting interesting. When a variable asks for 7 days before x, the aligner first pays from that variable's budget. Only the part the budget cannot cover shrinks the usable range for every variable. In the plan above:

The river series starts on the first shared day, so its backward budget is zero. The 30-day history costs 30 samples. The Lozzolo series starts in 2001, so its budget before the overlap is thousands of rows, and the 7-day rain history costs nothing. At the other end, Lozzolo stops on 31 December 2022 and the river continues to 4 April 2023. The 6 days of future rain cost 6 samples, and the 6-day river target costs nothing.

This is the reason the design uses budgets instead of a plain intersection. A long history on one station does not have to be cut to the shortest file. Only the windows that actually reach past a file's edge reduce the dataset.

Windows are sliding views over the aligned rows

After alignment, each variable has a start row, an end row and a window. get_convolution trims the array to the usable range and hands it to NumPy:

slider = self.top_slide(var_name) + self.bot_slide(var_name) + 1
# Trim data up to usable window in order to slide window over it and collect n samples
trimmed_data = self.var_vec.get_variable(var_name)[l_bot:l_top]
var = self.var_vec.get_variable(var_name)

if not np.iterable(var[0]):
    # Var is a 1-d numpy vector (var.ndim==1)
    conv_data = np.lib.stride_tricks.sliding_window_view(trimmed_data, slider)
else:
    e_size = np.size(var[0])
    views = np.lib.stride_tricks.sliding_window_view(trimmed_data, (slider, e_size))
    conv_data = np.array([np.concatenate(*w) for w in views])

sliding_window_view returns every window of the requested shape as a view into the original array, without copying the data [2]. For a variable whose elements are themselves vectors, such as the 24 hourly river readings, the window shape is (days, 24). Each window is then flattened into one row. The compile step stacks the rows of all input variables side by side into x, and all target variables into y.

Note the + 1 in slider. A window includes day x itself, so take 30 before x produces 31 days, and take 7 before x and 6 after x produces 14.

The logged numbers can be derived by hand

I reran PlotData.makefile on the data in the repository. The log starts like this:

> Changing field error to attempt_recovery
> Changing field proportions to [75, 25]
> Generated 1709 data points from input configuration
> Non target variables:['NormalizedRes', 'PNove'], Target variables: ['NormalizedResT']

Here is where 1709 comes from. The overlap runs from 23 March 2018 to 31 December 2022, which is 1745 days. The 30-day river history removes 30 samples at the start and the 6 days of future rain remove 6 at the end: 1745 − 30 − 6 = 1709.

The widths follow from the windows. The river history contributes 31 days of 24 hours, 744 values. The rain window contributes 14 daily values. That makes 758 columns in x. The target is 7 days of 24 hours, 168 columns in y. The saved files match: savefileN.npz holds x with shape (1283, 758), testN.npz holds (426, 758).

The log also draws the layout of a sample. This is the part of the log that helps later, when you have to build an input for the trained model by hand:

> DATA FED INTO THE MODEL:

> | NormalizedRes              | PNove
> | 30 before x           X    | 7 before x           X    6 after x
> <---->  ..28.. <----> <-^^-> <---->  ..5.. <----> <-^^-> <---->  ..4.. <---->

> EXPECTED OUTPUT OF THE MODEL:

> | NormalizedResT
> | Y    6 after y
> <-^^-> <---->  ..4.. <---->

The drawing makes one property of this plan visible. The target window starts at Y, the same day as X. Because windows include the current day, the river levels of day x appear in both the input and the target: 24 of the 168 target values are copies of input values. A model can learn to repeat them, which makes the error on the first day look better than the forecast is. The syntax has no offset keyword, so a target that starts on the next day needs a separately shifted variable prepared in .act.

The split is sequential

split NewPlan into NewFile, TestFile as 75, 25 does not shuffle. It cuts the arrays at a computed index:

props: list = self.specs['proportions']

secs = np.ceil(np.cumsum(np.array(props[:-1]) / 100.0) * np.size(x_data, axis=0) + 1).astype(np.int32)
x_s, y_s = np.split(x_data, secs), np.split(y_data, secs)

For the river plan, the cut is at ⌈0.75 × 1709 + 1⌉ = 1283, which gives the 1283 training and 426 test samples above. The + 1 moves every cut by one or two samples, so the proportions are approximate. For 1709 samples this does not matter.

For a time series, a sequential split is the correct choice. Hyndman and Athanasopoulos describe the standard practice: the test set is the most recent part of the series, so that the evaluation looks like forecasting data the model has not seen [3]. Here, the 426 test samples are the last 426 days of the aligned range.

For tabular data the same rule causes trouble. The Iris data is sorted by species. When I ran the Iris makefile on the unshuffled file, the log reported 150 samples and the split gave 106 for training and 44 for test. The 44 test samples were all Iris virginica, and the training set held only 6 of them. This is why the repository runs shuffle_iris.makefile first, as shown at the end of part 1.

Windows on rows that are not time

The breast cancer example shows what a window means when the rows are not ordered in time. It reads the Wisconsin breast cancer data [4] [5] after a shuffle, aligns everything on a row index, and asks for a window on one feature:

begin plan CancerPred expecting attempt_recovery
{
    align Thickness, CSize, CShape, MAdh against Indexes,Indexes,Indexes,Indexes as index
    align SECS, BNuclei, BChrom, NorNucleoli, Mitoses against Indexes, Indexes, Indexes, Indexes, Indexes as index
    align OHotTruth against Indexes as index

    consider x
    take x from Thickness
    ...
    take 4 before x from SECS
    ...
    make OHotTruth the target and take y from OHotTruth
    pair x with target
}
end plan

The mechanics work as before. The shuffled file has 697 rows. The 4-row window on SECS has no budget before the first index, so it costs 4 samples, and the log reports 693. The inputs are 8 single values plus a window of 5, 13 columns. The 30/70 split cuts at ⌈0.3 × 693 + 1⌉ = 209, and the saved arrays have shapes (209, 13) and (484, 13).

The meaning is another matter. After the shuffle, "4 before x" is the SECS value of four other, randomly chosen patients. The example demonstrates the syntax correctly, but the feature carries no information about patient x. A window is only meaningful when the alignment column orders the rows by something real, usually time.

Normalization happens before the split

The order of the sections has a consequence for the statistics. z_score runs in .act, before .make exists, so the mean and standard deviation come from the full series. In the river example, the scaling of the training data uses values from the test period. The scikit-learn documentation lists this as a standard form of data leakage: preprocessing must be fitted on the training data only [6]. The effect on a long hydrological series is probably small, but it is not zero, and the language makes it easy to overlook because the split is written on the last line of the file.

Fixing this properly means the split has to be known before the normalization, for example by letting .act primitives read a training range from the plan. That is a change to the language, not to one makefile.

A plan also carries settings

Besides the statements above, a plan has a small dictionary of fields. set changes a field by name, with the value in double quotes. No example in the repository uses it, so these two lines only show the form the parser accepts:

set RNND x name = "inputs"
set CancerPred fallback file = "fallback.npz"

x name and y name set the keys in the saved .npz. The expecting clause after begin plan sets the error mode. With attempt_recovery, a failed split writes all the data to the fallback file instead of stopping. Only attempt_recovery changes behaviour in the current code; full_recovery, used by the Iris example, is stored and logged but never read.

Reproducing the numbers

All runs in this post used NumPy 2.4.4. In that version, np.in1d, which the aligner calls, no longer exists, and np.isin is its replacement [7]. I patched it at runtime with np.in1d = np.isin before loading the interpreter. The river run reproduced the log stored in the repository exactly. The breast cancer log stored next to the data reports 691 samples, two fewer than my run, so the counts above are from my run. For Iris I replaced the placeholder paths with the unshuffled iris.data file from the repository.

The aligner is small because of one decision: after it finds the first shared key, it trusts row numbers. That works only if every source has exactly one row per day and every source writes its keys in the same way. Part 3 looks at what the repository does, outside the language, to make that true.