dear all,
I’d like simply to interpolate data following piecewise constant method, as illustrated in code below. I looked into available interpolation functions but I didn’t found proper configuration to do this. I would like extrapolation option “edgevalue” from interp1() but for interpolation… more or less the same behaviour as in plot2d2() graphic function but as function defintion.
The code is doing what I want but the loop makes it not performant when dealing with huge datasets.
Do you have any suggestion to make it clever and faster ?
Thanks
David
I=1:10
inom= [1 2 5 7 8 10]
innom= setdiff(I,inom)
iref= innom
for k=1:length(innom), iref(k)=inom(max(find(innom(k)>inom))) ; end
I(innom)=I(iref)
Hello David,
A piecewise constant function f on [a,b[ is defined e.g. by a sequence (x_i)_{i=1..n} where a=x_1<x_2<\dots<x_n=b and values (y_i)_{i=1..n-1} such that \forall x\in[x_i,x_{i+1}[, f(x)=y_i. Where are defined x and y in your code?
If you have huge sequences a built-in function should be preferable: you can use interp1 in nearest neighbor mode but in that case you should define abscicae values (u_i)_{i=0..n} such that (u_i+u_{i-1})/2=x_{i}, i=1\dots n, to synchronize the discontinuity (the nearest neighbor interpolant value changes at the middle of intervals).
Since the choice is not unique you have to define arbitrarily u_0, for example u_0=x_1-(x_2-x_1)/2 (it places u_1 right in the middle of [x_1,x_2]).
Then you have u_{i}=2x_{i}-u_{i-1}, i=1\dots n, and you can feed interp1 with (u_i)_{i=0..n} and (y_i)_{i=0..n}, where y_0 is arbitary, i.e. y_0=y_1.
EDIT : the method is not general, i.e. there are cases when there does not exist a nearest neighbor equivalent of an arbitrary piecewise C0 function…
S.
Hello Stephane,
the previous code was very specific to current problem,not general at all, I was focused on the index and somewhat confusing myself.
I tried to generalize a little bit but still far from being robust ! :
function yp= interp0(x,y,xp)
x_=xp
for k=1:length(x_), x_(k)=max(find(xp(k)>=x)) ; end
yp=y(x_)
endfunction
x=[1 2 5 7 8 10]
y= [1 2 5 7 8 10]+1.6
xp= [1:10] +0.1
yp= interp0(x,y,xp)
plot2d2(x,y)
plot2d(xp,yp,-2)
David
Sorry for my first complicated/not working solution. The following one is much clever…
function yp= interp0(x,y,xp)
i = floor(interp1(x,1:length(x),xp,"linear","edgevalue"));
yp = y(i);
endfunction
x = [1 2 5 7 8 10]
y = [1 2 5 7 8 10]+1.6
xp = linspace(-1,12,10000);
plot(xp,interp0(x,y,xp));
1 Like
Hi Stephane,
thank you for this simple and fast solution !
May this implementation of piecewise constant interpolation feature be included in future release of Scilab as an extra option of interp1() ?
best wishes,
David
Hi ,
looking for highest speed, I experimented successfully to use native function linear_interpn instead of interp1 macro, significantly faster on large dataset.
function yp= interp0(x,y,xp)
i = floor(linear_interpn(xp,x,1:length(x),"C0"))
yp = y(i)
endfunction
So the suggestion is to add the piecewise constant interplotation feature rather to linear_interpn .
best,
David
Hello,
I’d like to share use case experience with interp0 function above, in general when you need to regenerate a collection of the indexes where a user function meet some user condition along its definition space: these index may constitute intervals of consecutive indexes, intervals of variable lengths corresponding to alternate true and false value of the user condition. From previous analysis of the user function where one identified the lower bound and upper bound of each interval, we may want to regenerate the collection of indexes that verify the user condition in order to perform further calculation limited to extraction of user function values, and interp0 is useful and computationnally efficient for that.
Example:
V=grand(1,20, "def");
x= 1:20; // space
scf()
subplot(1,3,1)
plot(x,V)
//state function leads to iH and iL, arbitrarly for illustration :
iH= [1 11 18] //index begin of high state
iL= [4 14 19] //index begin of low state
//generate fstate as vector of indexes where 1 when state is high, 0 elsewhere
fstate= interp0([iH ; iL](:),[ones(iH) ; zeros(iL)], x );
subplot(1,3,2)
plot(x, fstate)
subplot(1,3,3)
//any calculation on extraction of indexes according to fstate
plot(tabul(V(fstate==1)))
This situation of retrieving indexes’ intervals (previously saved as [lowerbound;upperbound]) of variable length might be quite common in data processing tasks.
David
Yes nice, I think we already discussed this kind of usage one year ago : Datetime index localisation in matrix with dsearch? - #8 by davcheze
S.
