Problem with piece-wise function with logical indexing.

I created the following piece-wise fuction using logical indexing.

%Piece-wise function
F = @(coef,t) ((t < 0).*...
                     0 + ...
               (t >= 0 & t < 1200).*...
                    (1/coef(1) + t./coef(2) + 1/(coef(3))*(1-exp(-t))) + ...
               (t >= 1200).*...
                    (1200/coef(2) + 1/(coef(3))*(1-exp(-t)) - 1/(coef(3))*(1-exp(-(t-1200)))));
coef0 = [200, 200000, 100];

% Test that function is valid for wide range of t.
x = linespace(-10, 1300, 1310);
plot(x,F(coef0,x))

I must be missing somethign in the syntax. The function is invalid when t<~500. The rest of the graph looks correct. When I poke around, it seems that the very last part of the fuction:

exp(-(t-1200))

is being evaluated outside of it's defined region (t >= 1200). And when t< ~500, it blows up and goes to inf.

[madhan ravi]

Also did you notice the first 10 elements becomes NaN?? illustrate with a picture.

[Bryan Wilson]

Everything below 491 is NaN.

I think it's becuase that last part of the equation, exp(-(490-1200)) = exp(710) is too large to evalate. But I don't understand WHY matlab is evaluating this bit of the fuction when it's outside the specified range (t>=1200).

[Bryan Wilson]

If I comment out the entire last piece...

F = @(coef,t) (t < 0).*...
                    0 + ...
              (t >= 0 & t < 1200).*...
                   (1/coef(1) + t./coef(2) + 1/(coef(3))*(1-exp(-t)));% + ...
%             (t >= 1200).*...
%                  (1200/coef(2) + 1/(coef(3))*(1-exp(-t/coef(4))) - 1/(coef(3))*(1-exp(-(t-1200))));

Then everything else plots correctly, even below zero.

Answer 1

[Guillaume]

But I don't understand WHY matlab is evaluating this bit of the fuction when it's outside the specified range (t>=1200).

That expression:

exp(-(t-1200))

is evaluated for all t. It does not matter what it's multiplied by, even if it's (t >= 1200), it's still evaluated for all t, since that's what is in that particular expression.

As you have noticed it's going to generate -Inf for some t. Now your (t >= 1200) is going to generate 0 for the same values, and when you multiply 0 by +/-Inf you get NaN.

One possible workaround is to make sure that the input to exp is 0 for invalid t, so you could replace that particular expression by

exp(-(t-1200) .* (t >= 1200))

As long as you ensure that no term is +/-Inf then the mulitplication by (t >= 1200) will ensure the result is 0 outside the desired range.

Note that strictly speaking, there's no logical indexing here. You're using logical vectors, but are not indexing anything.

Also note, that your first part of the expression,

(t < 0).* 0 + ...

is a bit pointless. It's always going to result in a vector of 0. And adding 0 to something does not have much effect on that something.

[Bryan Wilson]

Thanks!

As for the first part (t < 0).*0, you're right, it's not necessary. It's just there to make things tidy for myself, to show that all t's are accounted for.

[Stephen23]

"As long as you ensure that no term is +/-Inf then the mulitplication by (t >= 1200) will ensure the result is 0 outside the desired range"

Whereas logical indexing would work for all values.

https://www.mathworks.com/help/matlab/math/array-indexing.html

[Guillaume]

I'm not convinced you could actually use logical indexing to construct that piecewise function in an anonymous function for a generic t.

If you assume that t is monotonically increasing, you could write

F = @(coef,t) [zeros(1, nnz(t<0)), ...%not logical indexing
               (1/coef(1) + t(t >= 0 & t <= 1200)./coef(2) + 1/(coef(3))*(1-exp(-t(t >= 0 & t < 1200)))), ... logical indexing
               (1200/coef(2) + 1/(coef(3))*(1-exp(-t(t >= 1200))) - 1/(coef(3))*(1-exp(-(t(t >= 1200)-1200))))]

If you use a plain function, then yes:

function out = F(coef, t)
    part1 = zeros(size(t));
    part2 = (1/coef(1) + t./coef(2) + 1/(coef(3))*(1-exp(-t)));
    part3 = (1200/coef(2) + 1/(coef(3))*(1-exp(-t)) - 1/(coef(3))*(1-exp(-(t-1200))));
    out = zeros(size(t));
    out(t<0) = part1(t<0);
    out(t>=0 & t<=1200) = part2(t>=0 & t<=1200);
    out(t>1200) = part3(t>1200);
end

I think Bryan's method is actually simpler.