[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

Martin Cermak <mcermak at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |ASSIGNED
                 CC|                            |mcermak at redhat dot com

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #1 from Martin Cermak <mcermak at redhat dot com> ---
Created attachment 9311
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9311&action=edit
proposed patch

The variance of N data points is V = S / (N - 1) where S is the sum of squares
of the deviations from the mean.  Here is an attempt to implement @variance()
operator using Knuth's algorithm [1]:

=======
def online_variance(data):
    n = 0
    mean = 0.0
    M2 = 0.0

    for x in data:
        n += 1
        delta = x - mean
        mean += delta/n
        M2 += delta*(x - mean)

    if n < 2:
        return float('nan')
    else:
        return M2 / (n - 1)
=======

This patch is based on current systemtap implementation of the aggregation
operators, which first pre-aggregates the data per each CPU (__stp_stat_add()),
and then, when the aggregations are actually being read via e.g. @sum (or
@variance), they are aggregated again, this time across all the CPUs
(_stp_stat_get()) and outputted.  This approach saves shared resources at the
collection time.  So, in this patch, per cpu variances are being collected
first and then they are being aggregated again across all the CPUs to give the
resulting @variance.  The N is assumed to be N >> 1 and so the resulting
@variance() is being computed as a simple mean of per-cpu variances.  Integer
arithmetic is being used.  With this patch, we get something relatively small
for data points closely spread along the mean, and something relatively big for
data points widely spread along the mean.  So it passes a rough sanity test:

=======
# stap -e 'global a probe oneshot { for(i=0; i<1000; i++) { a<<<42 } }  probe
end { printdln(", ", @count(a), @max(a), @variance(a)) }'
1000, 42, 1
# stap -e 'global a probe oneshot { for(i=0; i<1000; i++) { a<<<42 } for(i=0;
i<20; i++) { a<<<99 } }  probe end { printdln(", ", @count(a), @max(a),
@variance(a)) }'
1020, 99, 65
# stap -e 'global a probe oneshot { for(i=0; i<1000; i++) { a<<<i } }  probe
end { printdln(" ", @count(a), @max(a), @variance(a)) }'
1000 999 332833
# 
=======


-------------------------------------------
[1]
https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[fche at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #2 from Frank Ch. Eigler <fche at redhat dot com> ---
Just came across this:

https://gist.githubusercontent.com/kastnerkyle/f5ddbb9b647ba21577e5/raw/62cb0de74e4a8dff2732277e7e547e7715c18a1a/streaming_variance.py

and this:

http://web.cse.ohio-state.edu/~kamatn/variance.pdf

I'll do some more digging; maybe we have some even better algorithms available

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #3 from Martin Cermak <mcermak at redhat dot com> ---
Created attachment 9369
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9369&action=edit
working version of a patch

Attached patch computes the individual per CPU variances using the Knuth's
algorithm from Comment #1. Based on that, the aggregated variance over all the
CPUs is being computed using the "Total Variance" formula from the above paper.

This gives somewhat reasonable results for a "few" "small" integers with normal
distribution, but almost any other set of values makes it behave crazily
because of the integer arithmetic being used for the dividing.  Below I am
going to attach a little python script that helps comparing this stap variance
implementation with python's statistics.variance().

At the first glance, the floating point arithmetic inside the linux kernel
doesn't look like something usual or straightforward.  But an attempt to
implement it might be an interesting one.  Not sure about this though.

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #4 from Martin Cermak <mcermak at redhat dot com> ---
Created attachment 9370
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9370&action=edit
a little testing script

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[fche at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #5 from Frank Ch. Eigler <fche at redhat dot com> ---
Excellent progress!

> This gives somewhat reasonable results for a "few" "small" integers with
> normal distribution, but almost any other set of values makes it behave
> crazily because of the integer arithmetic being used for the dividing. 

I suspect the loss of precision is occurring in the per-cpu sd->variance =
stp_div64() ... code rather than the cross-cpu agg->variance one.  Your test
case (probe oneshot) runs only on one CPU, so a lot of the cross-cpu terms (S1)
should be zero.


> At the first glance, the floating point arithmetic inside the linux kernel
> doesn't look like something usual or straightforward.  But an attempt to
> implement it might be an interesting one.  Not sure about this though.

We'll eventually get -some- FP capabilities for stap code, but let's see how
much farther we can get without.

If you have any more time/interest in this problem, I'd suggest investigating
whether we could track the sd->avg / sd->_M2 / sd->variance values in a scaled
form.  We continuously track max/min, so have a good idea about the dynamic
range of the <<<'d values.  The code could track a "shift" parameter that
scales those numbers (via << shift) during the online update phase so as to
preserve as much precision as possible.


Another thing we should consider afterwards is bug #10234.  For a stap script
not interested in @variance, this code should not be run at all.

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #6 from Martin Cermak <mcermak at redhat dot com> ---
Created attachment 9413
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9413&action=edit
updated working version of a patch

Introduced the "shift" parameter (and fixed some errors).  Now the resulting
values correspond to those of python's statistics.variance() for all my tests
(more than one CPU core is effectively involved). Wheee!

Small shifts like 1, 2, or 3 significantly help.  Higher shifts bring slightly
more precision, but increase the risk of overflow.

Now the "shift" could either be updated on the fly, or it could be passed to
stap as a -DPARAM.  Given that small shifts help so nicely, ITSM the latter
might be a good trade-off, but not sure.

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #7 from Martin Cermak <mcermak at redhat dot com> ---
Created attachment 9414
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9414&action=edit
improved test program

Attached test program distributes the computations across available CPU cores.

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[fche at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #8 from Frank Ch. Eigler <fche at redhat dot com> ---
Comment on attachment 9413
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9413
updated working version of a patch

considered making the shift value a parameter for the @variance() operator,
kind of like @hist_linear()'s?  i.e.,  @variance(var,5) to suggest a shift of
2<<5 to the runtime for that variable?

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #9 from Martin Cermak <mcermak at redhat dot com> ---
Created attachment 9421
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9421&action=edit
updated working version of the patch

Introduced optional parameter for specifying the bit shift, so that both
@variance(var) and @variance(var, N) are valid.  Additionally -DSTAT_BIT_SHIFT
can be used, and has actual effect when N is not specified or is zero.

Now I plan to listen to the possible code review, clean the code up, fix a bug
residing in the bit shift being shared among all the stat operators (rewriting
one another ;-)  Write a dejagnu testcase, and do more testing.

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #10 from Martin Cermak <mcermak at redhat dot com> ---
Created attachment 9490
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9490&action=edit
The attached patch is just being tested.

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[mcermak at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

--- Comment #11 from Martin Cermak <mcermak at redhat dot com> ---
The @variance() part pushed as commit 63ead7fa4ed58e12efdff1b1ab59e6056c277371
.

-- 
You are receiving this mail because:
You are the assignee for the bug.

[Bug runtime/11308] aggregate operations for @variance, @skew, @kurtosis

[fche at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=11308

Frank Ch. Eigler <fche at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |RESOLVED
         Resolution|---                         |FIXED

--- Comment #12 from Frank Ch. Eigler <fche at redhat dot com> ---
Let's call it a day at @variance.

-- 
You are receiving this mail because:
You are the assignee for the bug.