bugGNU Octave - Bugs: bug #48774, Functions for moving statistics

 
 

bug #48774: Functions for moving statistics

Submitter:  Nir Krakauer <nir_krakauer>
Submitted:  Sat 13 Aug 2016 07:06:16 PM UTC
   
 
Category:  Octave Function Severity:  4 - Important
Priority:  5 - Normal Item Group:  Feature Request
Status:  Fixed Assigned to:  juanpi
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 21 Dec 2018 04:33:51 PM UTC, comment #31: 

@Juan: I find that is very hard to speculate about performance and that one is usually wrong.  The only way to be certain is to benchmark the two alternatives.  I have been testing with a matrix that I created with randi (1e3, 1000, 1000) and then saved to a file so the results are reproducible.  Here is an example benchmark which shows that "discard" is currently the very fastest of the endpoint options.


octave:1> load movmin.var
octave:2> tic; y = movfun (@mean, x, 51, "Endpoints", "shrink"); toc
Elapsed time is 6.17102 seconds.
octave:3> tic; y = movfun (@mean, x, 51, "Endpoints", "shrink"); toc
Elapsed time is 6.07471 seconds.
octave:4> tic; y = movfun (@mean, x, 51, "Endpoints", "fill"); toc
Elapsed time is 0.600617 seconds.
octave:5> tic; y = movfun (@mean, x, 51, "Endpoints", "fill"); toc
Elapsed time is 0.628939 seconds.
octave:6> tic; y = movfun (@mean, x, 51, "Endpoints", "discard"); toc
Elapsed time is 0.25172 seconds.
octave:7> tic; y = movfun (@mean, x, 51, "Endpoints", "discard"); toc
Elapsed time is 0.268769 seconds.


I would set up your own dataset and then check the timing pre/post-conversion to a dispatch structure.

Regarding copyright, I don't believe using a single line, in this case the call to movfun, qualifies the entire file as a derivative work.  Particularly, as there is only one real way to call the movfun function.  But, I'm not a lawyer so I think we need to get the advice of the Free Software Foundation.  I would have you ask jwe on how to proceed.

The author and creation date fields were more useful when we had primitive version control systems, or none at all.  Now, the author and modification date of every change is captured by Mercurial, including the original file creation, and these are less necessary.  I occasionally remove them when I come across them, but I haven't done a full purge.

Regarding _parser_movargs_, I think it depends on whether any further movXXX functions will be added.  I anticipated that there might be more, and so I placed _parser_movargs_.m together with movslice.m and movfun.m.  If there will only ever be statistics versions of sliding window functions then it can be made a private function in the statistics/ directory.

Rik <rik5>
Group administrator
Fri 21 Dec 2018 08:33:00 AM UTC, comment #30: 

Hi Rik,

Thank you for taking the time to answer all my comments.

1. I will merge your changes and go back to the dispatch structure, since I think it adds to maintainability. I do not think that the new information hinders its applicability at all. Aslo the working done in "discard" should be put in a "discard_bc" not to violate the concerns of the functions. Do you foresee a strike in performance by moving the code in the switch to a subfunction?

2. Regarding copyright. Since you were exposed and clearly saw my code on the wrappers I do not think that this is a case for Clean-room design. I will be very happy if my name appears in those wrappers so future users can link to me as well.

I also noticed that you removed author and creation date but many functions in octave core have those fields. Is this a personal criteria or is something the community wants to remove?

3. Isn't it better to have _parser_movargs_ as a private function in the scripts/statistics folder? So far only wrappers there are using it, right?

Regards,

Juan Pablo Carbajal <juanpi>
Group Member
Wed 19 Dec 2018 05:55:04 PM UTC, comment #29: 

@Juan: Regarding programming strategies in core Octave, there is no simple answer as to the best approach.  Octave should strike to have code which is readable and maintainable, but it also needs to balance that with performance.  The two goals are often in opposition.

If the only thing being implemented was a jump table then it would be fine to use the original strategy.  But, the current function does require more checking based on the type of Endpoint chosen.  For example, if the user specifies a numeric fill value then one can't use [arg "_bc"] to construct a valid function name so that is one special case.  Similarly, when the "Endpoint" property is "discard" there is no actual discard_bc necessary, but one does have to adjust various internal variables so that the correct length output is generated.  Once several of the input values required identifying and having different actions the dispatch table no longer made sense, at least too me.

I'm not opposed if you want to re-work the input validation, as long as the user-visible goals of functionality and performance are unchanged.

Rik <rik5>
Group administrator
Wed 19 Dec 2018 05:34:36 PM UTC, comment #28: 

@Juan: Regarding 'omitnan', it does not appear to be implemented so I added a warning.  Here is test code that I used.


x = (1:10)';
x(5) = NaN;
movmean (x, 3, 'omitnan')


Under Matlab, this produces


ans =

   1.5000
   2.0000
   3.0000
   3.5000
   5.0000
   6.5000
   7.0000
   8.0000
   9.0000
   9.5000


Using the original implementation from your repository


octave:8> movfun (@mean, x, 3, 'nancond', 'omitnan')
ans =

   1.5000
   2.0000
   3.0000
      NaN
      NaN
      NaN
   7.0000
   8.0000
   9.0000
   9.5000


For completeness, running with 'includenan' in Matlab


movmean (x, 3, 'includenan')
ans =

   1.5000
   2.0000
   3.0000
      NaN
      NaN
      NaN
   7.0000
   8.0000
   9.0000
   9.5000



Rik <rik5>
Group administrator
Wed 19 Dec 2018 05:17:20 PM UTC, comment #27: 

@Juan: Regarding window length, I specifically added support for even-length windows.  Matlab already has this, and users will want this feature.  Since they are now supported there is no requirement to have input validation which excludes it.

The input validation for wlen > 1 is now done in movslice.m.  The same input validation of wlen was required for both movfun.m and movslice.m.  For performance reasons, it made sense to remove that validation from movfun and delegate it to movslice so that the validation is done just once.

Rik <rik5>
Group administrator
Wed 19 Dec 2018 05:11:16 PM UTC, comment #26: 

@Juan: Regarding copyright, I re-formatted the copyright block to use the same line breaks, and https instead of http in links, as the rest of Octave m-files.  Everything else about the license (GPL3) is the same.

I touch a lot of the Octave code and therefore I usually don't bother to add my name to the copyright holder list.  I agree that I made significant enough changes to movfun.m that I could be listed as an additional copyright holder, but I didn't feel like it.

For the movXXX functions, I did not use the code from your repository.  I wrote a blank function template with copyright block, documentation, nargin checking, call to movfun, and BIST tests for input validation.  I then used perl to change this template in to the many different movXXX functions.  After the base change I went in and edited the documentation.  For example, 'movmean' says "Calculate the moving average over a sliding window"  It does not say "Calculate the moving mean" because that is not the natural phrasing in English.  Since I created these functions separately I assigned myself copyright.  The possible exception to this is movmin.m which did have an initial implementation.  If you would like copyright for that m-file that would be fine.

Rik <rik5>
Group administrator
Wed 19 Dec 2018 01:55:33 PM UTC, comment #25: 

As a concrete example of the problem mentioned in the previous comment, see this snippet:

  ## Obtain function for boundary conditions
  if (isnumeric (bc))
    bcfunc = @replaceval_bc;
    bcfunc (true, bc);  # initialize replaceval function with value
  else
    switch (tolower (bc))
      case "shrink"
        bcfunc = @shrink_bc;

      case "discard"
        bcfunc = [];
        C -= length (Cpre);
        Cpre = Cpos = [];
        N = length (C);
        szx(dperm(1)) = N;


We have branching, function handle, and direct implementation on bc inside the switch. If you are a potential developer coming to this function this is hell, becuse is really hard to track the path of the function (the potential dev might not want to become an expert of movfun). Why do it like this when this can be avoided and become a single line?

Separation of concerns also applies here, the mian function now is not only dispatching but also implementing bc, why?

To add a new bc, in the current form,  one needs to edit at three points: the cell with valid_bc, the switch part (making sure the if branching is not introducing an block), and add the code of the subfunction.
With the dispatch method, only the cell of valid_bc needs to be extended, and the code of the subfunction added.

Juan Pablo Carbajal <juanpi>
Group Member
Wed 19 Dec 2018 01:33:57 PM UTC, comment #24: 

Hi Rik,

Also i see you have removed the dispatch structure for the bc I implemented. This always leads to code that has lengthy argument parsing, and decreases maintainability of functions, besides the original author.
This is one of the main problems in pkg.m and we should stop doing this. I presented this at Octconf, it is a itty you missed that.
Also speaking with other developers the dispatch method seems to be preferred.

Unless you can argument pro such a change, I would put back the dispatch structure.


Juan Pablo Carbajal <juanpi>
Group Member
Wed 19 Dec 2018 01:30:27 PM UTC, comment #23: 

Hi Rik,

More on the port. th version of Octave has this check which in my versino of the function s unnecesary, omitnan does work


if (omitnan)
    warning ('movfun: "omitnan" is not yet implemented, using "includenan"');
  endif


Also it would be good to add identifiers to the warnings

Juan Pablo Carbajal <juanpi>
Group Member
Wed 19 Dec 2018 01:28:27 PM UTC, comment #22: 

Hi Rik,

I am merging your changes ot the external repo and I can't fin the following checks


   # Check for proper window length
    # TODO: Matlab accepts even windows
    if (mod (wlen, 2) == 0)
      error ('Octave:invalid-input-arg', 'Window must be of odd length');
    endif
    if (wlen == 1)
      error ('Octave:invalid-input-arg', 'Window must be larger than 1');
    endif


Where are they? The algortium won't work correctly for even window length, and the result for length 1 is not correct due to the changing behavior of functions that run on columns (e.g. mean of an 1-by-n array will give a single scalar)

Juan Pablo Carbajal <juanpi>
Group Member
Wed 19 Dec 2018 01:00:26 PM UTC, comment #21: 

Hi Rik,

Regarding comment #17. It is a good tought, hoever, the argout concept relates to functions with multiple output arguments, not functions with a single output argument of multiple dimensions. So I think using "argout" would add to the users confusion between these two type of functions. Maybe "dimargout" is more intuitive, but seems overly long...

Juan Pablo Carbajal <juanpi>
Group Member
Wed 19 Dec 2018 11:16:16 AM UTC, comment #20: 

... and I would like to add that you more than deserve having your name in the copyright of movfun.m. I see the improvements you have done to the doc and the demos. Thanks

Juan Pablo Carbajal <juanpi>
Group Member
Wed 19 Dec 2018 10:46:37 AM UTC, comment #19: 

Hi Rik,

Thanks for the transfer. However I am not sure you are respecting the license in the code I shared. I understand you did a lot of edition of the code, but you used mainly my work. Hence I do no see removing my name from the copyright noticed as justified.

Am I misunderstanding something?

Juan Pablo Carbajal <juanpi>
Group Member
Wed 19 Dec 2018 12:33:09 AM UTC, comment #18: 

Finally, after a lot of work on movfun, I added the moving statistics function.  See https://hg.savannah.gnu.org/hgweb/octave/rev/f73ca7468864.

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Tue 18 Dec 2018 10:44:38 PM UTC, comment #17: 

@Juan: Thanks, I see what is going on now.  Perhaps a better name for this would be "argout" in analogy with the function nthargout()?.  I will have to adjust the documentation.

Rik <rik5>
Group administrator
Tue 18 Dec 2018 09:41:24 PM UTC, comment #16: 

Hi Rik,

I see... the issue is that in the assumptions of the code mean is not a function returning a multidimensional outpout. It applies it self over multiple input dimension, but the calculation itself is a scalar.
'outdim' is used when the function itself goes form R^m --> R^D, it is meant to select one of the D dimensions of the output.

I agree that in the general matlab/octave jargon this is not directly interpreted like this. That is, one would be lead to think that mean has a multidim output because when you apply it to a multidim input it keeps the other dimensions it doesn't act upon, e.g. size(x) == [5 10 2], then size (mean (x)) == [1 10 2] (or that squeezed). But the calculation itself is from R^m (set of scalars) to a R (a scalar), so the function is a scalar function.

This is why movfun is telling you that mean doesn't returns the 3rd dimension, which you are asking for.

If we can discuss this a little further, I will be glad so I can get an idea how to put it in more clear terms in the docs.

On the internal working, movfun puts the desired dimension as the first dimension of a new array (permute) and then join all other dimension as columns (reshape). This is why the function is asked to work on columns, and not on multidimensional arrays.

The idea behind this convention is to allow the user to think only on the atomic behavior of their function (how it behaves when you give a column vector to it), and not on its behavior for multidimensional arrays.

Try this example, as an illustration

x = magic (3);
x(:,:,2) = 3*x
y = movfun (@statistics, x, 3, 'outdim', 6); % returns movmean
y_ = movmean (x, 3);
assert(y,y_)


When statistics is presented with a column vector it returns a vector with 9 elements (a 9-dimensional output), we choose the 6th one. The function mean, when aplied to a column vector returns a scalar, hence you can't choose the 3rd dimension.

I hope this clarifies some thing.

Juan Pablo Carbajal <juanpi>
Group Member
Tue 18 Dec 2018 08:22:50 PM UTC, comment #15: 

@Juan: Can you explain the "outdim" parameter?  The documentation right now says


"outdim"
     A row vector that selects which dimensions of the calculation
     will appear in the output Y.  This is only useful when FCN
     returns an N-dimensional array in Format 1.  The default is to
     return all output dimensions.

Programming Note: The property "outdim" can be used to save memory
when the output of FCN has many dimensions, or when a wrapper to
the base function that selects the desired outputs is too costly.
When memory is not an issue, the easiest way to select output
dimensions is to first calculate the complete result with 'movfun'
and then filter that result with indexing.  If code complexity is
not an issue then a wrapper can be created using anonymous
functions.  For example, if 'basefcn' is a function returning a
K-dimensional row output, and only dimension D is desired, then the
following wrapper could be used.

     FCN = @(x) basefcn (x)(:,size(x,2) * (D-1) + (1:size(x,2)));
     Y = movfun (@fcn, ...);


I tried this as an example, but there was an error.


x = magic (3);
x(:,:,2) = 3*x
x =

ans(:,:,1) =

   8   1   6
   3   5   7
   4   9   2

ans(:,:,2) =

   24    3   18
    9   15   21
   12   27    6

movfun (@mean, x, 3, 'outdim', 3)
error: movfun: output dimension OUTDIM (3) is larger than largest available dimension (1)
error: called from



Rik <rik5>
Group administrator
Tue 18 Dec 2018 06:14:19 AM UTC, comment #14: 

Hi all,

I have automatically generated wrappers to the functions in the list provided: https://bitbucket.org/KaKiLa/movfun/src/default/

There are no tests, demos, and I have actually only tested the functions with the following script


fun ={
'movmad',
'movmax',
'movmean',
'movmedian',
'movmin',
'movprod',
'movstd',
'movsum',
'movvar'
};

x = ones (10, 2);
for i=1:numel(fun)
  cmd = sprintf ('%s(x, 3)', fun{i});
  printf ([cmd '\n'])
  eval (cmd);
endfor


I am going lazy here an I will wait for bug reports to fix the wrappers, if needed.


Juan Pablo Carbajal <juanpi>
Group Member
Wed 12 Dec 2018 08:29:29 AM UTC, comment #13: 

Juan: I did not miss your argument. As I said your functions
apply to ordered data.

In any case Rik is handling integrating your functions
appropriately.

And, thanks very much for your contributions!

Michael Godfrey <godfrey>
Group Member
Wed 12 Dec 2018 02:29:19 AM UTC, comment #12: 

To comment #8:

I am working on a small extension that allows to use functions that return an array instead of a scalar. This implies one could use, for example, polyfit on a moving window.

This doesn't affect the compatibility with matlab, as the wrappers could ignore this option and just use the default arguments (output dimension == 1)

I will do these changes in the external repo, and I can merge them into octave later on, when the locaction and changes to movfun in octave are settled.

To comment #9:
I think you missed my argument. The fact of using a moving window implies that you consider your data ordered into a sequence by a monotonic growing index, be it time or anything else with this property. Hence, applying a moving window filter, is equivalent as stating that your data is a signal.
Statistical data with independence of samples, can be ordered in anyway you like (that is, they are not ordered); that is the core of independent samples. Therefore, results that depend on a arbitrary order are weird, hard to interpret, if not plainly inapposite.


Juan Pablo Carbajal <juanpi>
Group Member
Tue 11 Dec 2018 09:46:53 AM UTC, comment #11: 

Compatibility with Matlab seems right in this case.

Michael Godfrey <godfrey>
Group Member
Mon 10 Dec 2018 09:55:02 PM UTC, comment #10: 

These functions should stay in core, rather than the statistics package, to match Matlab.


movmad     Moving median absolute deviation
movmax     Moving maximum
movmean    Moving mean
movmedian  Moving median
movmin     Moving minimum
movprod    Moving product
movstd     Moving standard deviation
movsum     Moving sum
movvar     Moving variance



Rik <rik5>
Group administrator
Mon 10 Dec 2018 09:04:26 PM UTC, comment #9: 

These functions compute functions of an ordered sequence.

The signal processing package computes functions of
time sequences.

The statistics package computes functions of variables
which may be ordered in various ways.

Unless it is intended to restrict the "move" functions
to time series, they should go in the statistics package.

Michael Godfrey <godfrey>
Group Member
Mon 10 Dec 2018 08:11:25 PM UTC, comment #8: 

I think we should try to get this in to Octave 5.0.

I will see about pulling the code in on the weekend, maybe sooner.

Rik <rik5>
Group administrator
Mon 10 Dec 2018 11:55:16 AM UTC, comment #7: 

Hi all,

I have been developing  the core of all these functions (octave-native). You can get the code from here: https://bitbucket.org/KaKiLa/movfun/src

Writing all the mov* functions should be a matter of writing the corresponding wrapper.

I hope that in the future that code will be shipped with octave.

I actually believe that it is brain-death to ship moving windows functionalities with a statistical package that considers independent samples. Applying a moving window to analyze data is tantamount of a prior of non-independent samples. Hence these functions fit better in the signal package. This is reflected in the demos.

To illustrate, consider the results


x  = randn (1e3, 1);
y  = movfun (@mean, x, 25);
xp = x(randperm (1e3));
yp = movfun (@mean, xp, 25);
plot(y, 'o', yp , 'x');


What's the correct moving sample result?

Bechmarks against the C++ versions would be interesting.

Juan Pablo Carbajal <juanpi>
Group Member
Thu 12 Jul 2018 07:45:06 PM UTC, comment #6: 

These should go in to core Octave since Matlab users will expect to find them there.  I don't think anyone has had a chance to review the files here.

Rik <rik5>
Group administrator
Thu 12 Jul 2018 07:30:27 PM UTC, comment #5: 

Will these appear in octave, or in the statistics package ?

John Donoghue <lostbard>
Group Member
Mon 02 Apr 2018 03:30:34 PM UTC, comment #4: 

The Octave Maintainers are deep in to the release procedure for 4.4 so we won't be able to review this for several more weeks, but we should have these functions in Octave.

Rik <rik5>
Group administrator
Mon 02 Apr 2018 12:22:41 PM UTC, comment #3: 

Sorry, I found a bug in movmax and movmin and attached the new patch here.

(file #43767)

JunWang <junwang>
Sun 01 Apr 2018 04:24:40 AM UTC, comment #2: 

I am trying to implement a part of functions for moving statistics. Movsum, movprod, movmax and movmin have been implemented.

(file #43755)

JunWang <junwang>
Sun 14 Aug 2016 02:54:36 PM UTC, comment #1: 

Yes, this would be nice.  In the mean time I added the new functions to the _unimplemented_ list (http://hg.savannah.gnu.org/hgweb/octave/rev/ed023556d4fa)

Rik <rik5>
Group administrator
Sat 13 Aug 2016 07:06:16 PM UTC, original submission:  

Matlab 2016a has introduced movmean, movsum, movmedian, movmax, movmin, movvar, movstd, which seem like they would be useful for us to have too. They are not yet listed in _unimplemented_.m

Nir Krakauer <nir_krakauer>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #43767:  movsum_movprod_movmax_movmin.patch added by junwang (63KiB - application/octet-stream)
file #43755:  movsum_movprod_movmax_movmin.patch added by junwang (63KiB - application/octet-stream - movsum, movprod, movmax and movmin)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by juanpi (Posted a comment)
  • -email is unavailable- added by lostbard (Posted a comment)
  • -email is unavailable- added by junwang (Updated the item)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by nir_krakauer (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-12-19 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2018-12-10 rik5 Severity3 - Normal 4 - Important
    2018-12-10 juanpi StatusPatch Submitted In Progress
        Assigned toNone juanpi
    2018-04-02 rik5 StatusConfirmed Patch Submitted
    2018-04-02 junwang Attached File- Added movsum_movprod_movmax_movmin.patch, #43767
    2018-04-01 junwang Attached File- Added movsum_movprod_movmax_movmin.patch, #43755
    2018-03-21 pantxo Dependencies- bugs #53402 is dependent
    2016-08-14 rik5 StatusNone Confirmed
        SummaryMoving statistics functions Functions for moving statistics

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code