How to vectorize a piecewise periodic function in MATLAB?

Question

I've noticed that matlab builtin functions can handle either scalar or vector parameters. Example:

sin(pi/2)
ans =
     1

sin([0:pi/5:pi])
ans =
         0    0.5878    0.9511    0.9511    0.5878    0.0000

If I write my own function, for example, a piecewise periodic function:

function v = foo(t)
t = mod( t, 2 ) ;

if ( t < 0.1 )
   v = 0 ;
elseif ( t < 0.2 )
   v = 10 * t - 1 ;
else
   v = 1 ;
end

I can call this on individual values:

[foo(0.1) foo(0.15) foo(0.2)]
ans =
         0    0.5000    1.0000

however, if the input for the function is a vector, it is not auto-vectorized like the builtin function:

foo([0.1:0.05:0.2])
ans =
     1

Is there a syntax that can be used in the definition of the function that indicates that if a vector is provided, a vector should be produced? Or do builtin functions like sin, cos, ... check for the types of their input, and if the input is a vector produce the same result?

Answer 1

You need to change your syntax slightly to be able to handle data of any size. I typically use logical filters to vectorise if-statements, as you're trying to do:

function v = foo(t)

v = zeros(size(t));
t = mod( t, 2 ) ;

filt1 = t<0.1;
filt2 = ~filt1 & t<0.2;
filt3 = ~filt1 & ~filt2;

v(filt1) = 0;
v(filt2) = 10*t(filt2)-1;
v(filt3) = 1;

In this code, we've got three logical filters. The first picks out all elements such that t<0.1. The second picks out all of the elements such that t<0.2 that weren't in the first filter. The final filter gets everything else.

We then use this to set the vector v. We set every element of v that matches the first filter to 0. We set everything in v which matches the second filter to 10*t-1. We set every element of v which matches the third filter to 1.

For a more comprehensive coverage of vectorisation, check the MATLAB help page on it.

Answer 2

A simple approach that minimizes the number of operations is:

function v = foo(t)
t = mod(t, 2);
v = ones(size(t)) .* (t > 0.1);
v(t < 0.2) = 10*t(t < 0.2) - 1;
end

If the vectors are large, it might be faster to do ind = t < 0.2, and use that in the last line. That way you only search through the array once. Also, the multiplication might be substituted by an extra line with logical indices.

Answer 3

I repeatedly hit the same problem, thus I was looking for a more generic solution and came up with this:

%your function definition
c={@(t)(mod(t,2))<0.1,0,...
@(t)(mod(t,2))<0.2,@(t)(10 * t - 1),...
true,1};
%call pw which returns the function
foo=pw(c{:});
%example evaluation
foo([0.1:0.05:0.2])

Now the code for pw

function f=pw(varargin)
for ip=1:numel(varargin)
    switch class(varargin{ip})
        case {'double','logical'}
            varargin{ip}=@(x)(repmat(varargin{ip},size(x)));
        case 'function_handle'
            %do nothing
        otherwise
            error('wrong input class')
    end
end
c=struct('cnd',varargin(1:2:end),'fcn',varargin(2:2:end));
f=@(x)pweval(x,c);

end

function y=pweval(x,p)
todo=true(size(x));
y=x.*0;
for segment=1:numel(p)
    mask=todo;
    mask(mask)=logical(p(segment).cnd(x(mask)));
    y(mask)=p(segment).fcn(x(mask));
    todo(mask)=false;
end
assert(~any(todo));
end