Running finite sum

[yves-surrel]

Any idea about using GenericScanKernel to make a running sum of the preceding N terms, i.e. S(k) = Sum(i = k - N, k)input[i] instead of S(k) = Sum(i = 0, k)input[i] (e.g. for a window averaging) ?

[inducer]

Isn't that just

s = scan(input)
result[i] = s[i] - s[i-N]

? (with appropriate handling of the left boundary of course)

[yves-surrel]

The problem is to deal with the loss of resolution you get when accumulating millions of values (all GPU cards do not accept double precision). Of course you are right (simply subtracting two scans once they are done), but it gets more and more noisy towards the end of the scans. My question was: can this be handled within the GenericScanKernel mechanism...

[inducer]

Well, if all you need is a few extra bits in the significand, why not use some sort of double-double trick on float32?

https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic

[kif]

Here is a small OpenCL library which does it for basic operations (with academic references).
https://github.com/silx-kit/silx/blob/master/src/silx/resources/opencl/doubleword.cl

[yves-surrel]

Very interesting reference.
However, I probably misunderstood @inducer's answer. It is possible to handle the fixed length scan within the input statement of a GenericScanKernel, so I will not need to switch to real or emulated double precision. Nevertheless, I keep your reference.

[yves-surrel]

To overcome this, I am trying to use a segmented scan. However, I can see nowhere any simple working example. I have seen issue #411 and your comment, which is hard to understand (what does ' the scan expression must act as if the a input to the expression were not supplied' mean ? If, as you say, 'you can't use [across_seg_boundary] to set a flag at exactly the segment boundary', what could be the use of this expression???

How should I modify the following code ?

import numpy as np
from matplotlib.pyplot import plot, title, imshow
import pyopencl as cl
import pyopencl.array as cla
from pyopencl.scan import GenericScanKernel

ctx = cl.create_some_context()

queue = cl.CommandQueue(ctx)

out_buf = cla.empty(queue, (10000,), dtype=np.int64)

knl = GenericScanKernel(
    queue.context, np.int64,
    arguments=" __global long *outBuf",
    input_expr="1",
    scan_expr="across_seg_boundary ? 0 : a+b", neutral="0",
    is_segment_start_expr="i % 1000 == 0",
    output_statement="outBuf[i] = item")

knl(out_buf)

out_buf = out_buf.get()

plot(out_buf.flatten())

[yves-surrel]

Here is a working example:

"""
Created on Thu Mar 31 11:25:45 2022
@author: yves

Example of use of GenericScanKernel to calculate a running average.

In this example , the average is done over the 'abscissa', i.e. the index

A moving window sweeps across the buffer, and the indices within the window
are averaged.

The initial state is when the center of the window is at the first point of
the buffer, and the final state is when it is at the last point of the buffer

"""

import numpy as np
from matplotlib.pyplot import plot
import pyopencl as cl
import pyopencl.array as cla
from pyopencl.scan import GenericScanKernel

ctx = cl.create_some_context()

queue = cl.CommandQueue(ctx)

width = 50
kernel_width = 11 # Prefer odd numbers

out_buf = cla.empty(queue, (width + kernel_width,), dtype=np.float32)

# Buffer to accumulate Sigma(1), Sigma(i)
stats_buf = cla.empty(queue, ((width + kernel_width) * 2,), dtype=np.float32)


preamble = """
static float2 load(int iR, int width, int kernelWidth)
{
     // Input: iR, 'right' pixel, will accumulate
     
     // Index of 'left' pixel, de-accumulate
     int iL = iR - kernelWidth;
     
     // To be accumulated by right pixel  scan
     float2 resultR;
     // To be de-accumulated by left pixel  scan
     float2 resultL;
     
     // Accumulate if right pixel is within range
     // Return increment and abscissa as a float2
     resultR = iR < width ? (float2)(1, iR) : 0;
     
     // De-accumulate if left pixel is within range
     // Return increment and abscissa as a float2
     resultL = iL >= 0 ? (float2)(1, iL) : 0;
     
     // Output to scan sum
     return resultR - resultL;
}
"""

knl = GenericScanKernel(
    queue.context, cl.cltypes.float2,
    arguments=" __global float2 * statsBuf, __global float *outBuf, const int width, const int kernelWidth",
    input_expr="load(i, width, kernelWidth)",
    scan_expr="a+b", neutral="0",
    preamble=preamble,
    output_statement="outBuf[i] = item.s1 / item.s0")

knl(stats_buf, out_buf, width, kernel_width)

out = out_buf.get()[kernel_width//2:-kernel_width//2]

plot(out)