CUMSIM: "CUMulative SIMpson" integral

Integral(y(x),x(1),x(end)) by Simpson. Nonequidistantlysampled data regularuzed by PARFIL

Description

CUMSIM integrates sampled vectors or matrices y as functions of vector x, not necessarily uniform. Parabolic interpolation is made by PARFIL.

Contents of CUMSIM.zip

CUMSIM
CUMSPLINE
COMPCUM

CUMSPLINE: similar to CUMSIM, but interpolation is made by interp1(x,y,X,'spline'), where X=parfil(x,x);

COMPCUM: compares results of
 
CUMSIM,
CUMSIMPSON (cf. FEX) and
CUMSPLINE,

integrating sin(x) with randomly generated x.
Repetitive call of COMPCUM generates different versions of the screenshot above.

As a rule, CUMSPLINE wins, but if the data are slightly noisy, CUMSIM can win (cf. PARFIL).

• MATLAB release: MATLAB 7 (R14)
• Other requirements: PARFIL

[John D'Errico]

Please keep these battles off the file exchange.
John

[Vassili Pastushenko]

Hey John!
I have submitted this file on monday, before the PARFIL has appeared on FEX!!!

Now what does your voting for deletion mean? You don't like that my program works better than that of your friend? You are a battler, not me. Sorry, but I do not like to be accused for the crimes done by others.

[John D'Errico]

No. This is not a question of friend or foe. I note that you chose to compare this tool to only one of at least 3 alternatives that I can think of offhand that do exactly the same thing. One is in matlab, the other two are on the file exchange. Only one of the three was written by Duane. Given that this tool seems not to offer any significant enhancement, I can hardly take any other inference. Please behave professionally on the file exchange. Take any feuds offline.

[Vassili Pastushenko]

.

I have seen only CUMSIMPSON in "integration" as a comparable alternative to CUMSIM.

 Another program, CUMSIMPSUM works only at equally spaced data and cannot be directly compared in the case of irregularly spaced data. CUMSPLINE was written by me to show that the results depend on quality of interpolation, submitted together with CUMSIM. I pretend however that CUMSIM is a TOP QUALITY product within the frames of QUADRATIC INTERPOLATION. What is the third alternative? If you vote for deletion of CUMSIM, I expect that you know something better, other possible reasons would be not objectively correct.

[John D'Errico]

After offline discussions with him, I'll accept that Vassili placed this tool on the exchange as an enhancement to what was there, rather than for other motives.

The cumulative spline tool in here is slower than it needs to be, by about 33% due to a lack of vectorization.

Cumsim itself, as it is based on parfil, does appear to be often faster than its alternatives, while retaining accuracy. The claim of higher accuracy in the presence of minor noise is probably due to the more local nature of parfil than a spline. A C2 spline is known to have higher interpolation variance in the presence of noise than other local mehods. The fact that parfil does not solve a system of linear equations explains both the speed gain and the lowered noise component.

I'll argue the subtle point that a Simpson integral of a parfil interpolant makes it not truly a cumulative simpson, but perhaps a near relative.

Finally, when submitting the file, the file exchange ID number of those related files which one must have to run this code should be included in the "Inspired by" edit box. This creates a link to those files with no other effort required. In this case, had the number 10229 been typed intop that box when he submitted this tool or when it is updated, the user would find a direct link to parfil, which is not included in this submission. Likewise, there is also a functin o_o this tool uses. I tried a search on the file exchange but it did not turn up. Its silly to have to edit all of his functions with a find to remove any references to a function that was never necessary anyway. I'd have rated this tool higher but for these flaws and the unnecessary slowness of the cumulative spline code. I'll re-rate it if they are repaired.

[Vassili Pastushenko]

Yes, I recognize that in both cases, in CUMSIM and in CUMSPLINE, a simple vectorization can be done. For instance, in cumspline all commands starting from the first FOR loop, can be substituted by 1-liner
which represents the proper Simpson integration block:
CS(2:ROWS,:)=cumsum((Y(1:2:N-2,:)+4*Y(2:2:N-1,:)+Y(3:2:N,:)).*repmat(diff(x(:)),1,COLS)/6);

Somehow intuitively I did not like it, but now I have checked it.
Indeed, if Y is only one column, this saves about 30% time, as John sais.

However, if I have many columns, such an improvement can act according to the proverb "better is an enemy of good". Take an example
 x=sort(5000*rand(1,10000))';
y=[sin(x) cos(x) exp(-x)];
y=repmat(y,1,20);
[CS,CS2]=cumspline(x,y);
t = 0.172 0.156
The first time here is the calc. time of vectorized version, the second one is the time of the criticized version, which appears
to be better than that of the vectorized one!

 I am not sure, whether such a vectorization is meaningful here. Again, we have relatively big data arrays (10000*60), and we discuss the time differences about 0.02 sec. Does it make sense?

Here my timing in CUMSIM, if any further questions:

just after interp1:

tic
CS2(2:ROWS,:)=cumsum((Y(1:2:N-2,:)+4*Y(2:2:N-1,:)+Y(3:2:N,:)).*repmat(diff(x(:)),1,COLS)/6);
t(1)=toc;

tic
for i=1:3
 IND(:,i)=i:2:N+i-3;
end
FIL=[1 4 1]'/6;
  
CS=zeros(ROWS,COLS);
D=diff(x);

for i=1:COLS
    w=Y(:,i);
CS(2:ROWS,i)=(w(IND)*FIL).*D;
end
CS=cumsum(CS);

t(2)=toc

[Vassili Pastushenko]

Well It is known that preallocation can help to save time. I have mofied the vectorization by adding preallocation:

tic
CS2=zeros(ROWS,COLS);
DX=repmat(diff(x(:))/6,1,COLS);
CS2(2:ROWS,:)=cumsum((Y(1:2:N-2,:)+4*Y(2:2:N-1,:)+Y(3:2:N,:)).*DX);
t(1)=toc;
the rest as before.

I obtain

t = 0.172 0.125

So it does not help too much, and the measured time fluctuates so that I could conclude that my stupid version is by about 30 % better than the vectorized one. I have some arguments concerning ML speed, but I am not sure that here we have a proper place for such a discussion.

In the previous post a misprint: timing was not in cumcim, but in cumspline, as is also clear from the text.