Logic Statements in Function Handle

I am trying to perform a logic statement in a function handle, so that I can evaluate a function at a series of coordinates, then make decisions.

I have the following:

p --> an N x 2 Double array. Contains pairs of coordinates [x,y]

Fi --> a griddedInterpolant function. The coordinates in p fall within range and domain of the grid that was interpolated to obtain this (so they are compatible).

Fi2 --> a second griddedInterpolant function, derived from a second set of data more or less the same as in Fi.

Important point: the values that can be evaluated in Fi and Fi2 are partially mutually exclusive. If Fi is negative, the Fi2 must be positive, and the other way around. However, it is possible for both to be positive (there is a gap between them).

I've used this to make a function handles as follow:

fd=@(p) Fi(p);
fd2=@(p) Fi2(p);

These work just fine for calling later on and give nice sensible values when evaluated at p or sub-arrays of p.

Here's the tricky bit. I want to make a new function handle ("fh") that takes the evaluations of these functions at p, and then makes a step-change based on certain results for each individial p entry. Essentially, I want a function handle that does the work of an If statement. So, I tried to accomplish this by making 2 function handles for evaluating, and then treating them as logic statements, so I am multiplying by 1 or 0.

What I have tried (simplified), is this:

fha=@(p) (feval(fd,p));
fhb=@(p) (feval(fd2,p));
fh=@(p) max((fha<(0))*A, (fhb)<(0)*B, ((fha)>(0)*(fhb)>(0))*((A+B)/2));

This should be the equivalent of:

If fha < 0, then A; elseif fhb < 0, then B; else average of A and B

In other words, if a point within p has coordinates in a region where fd is negative, that point will have an associate value A. If it is in a position where fd2 is negative, it will be B. But if it falls in the regions where both evaluate to positive, its assigned value is the average of A and B.

I've also tried this nesting the feval statements into the fh call, and doing it as a single handle. Neither has worked. Is there a smart way to approach this, keeping in mind that I need to pass this handle to another function and then do lots of array manipulation stuff with it?

[Torsten]

For such complicated operations, you shouldn't use function handles, but proper functions in which you can work with if-statements and all that is needed.

[Geoffrey Rivers]

Ok. I'll fork my project and work on a version where the other points in the code can handle this as a separate function. Hopefully it's easy to implement.

But, in the interest of my education: even though I shouldn't do it this way, is there a way to? If not, why not?

[Torsten]

even though I shouldn't do it this way, is there a way to? If not, why not?

Why care if it is much easier to program and to understand if you use a function ? You can use functions as flexible as you can use handles.

But maybe someone else will invest the effort to give you a handle solution.

[Stephen23]

I second the recommendation to use a simple function in a file. A few reasons to avoid stacks of anonymous functions:

• inefficient: multiple function calls, multiple workspaces, all data gets evaluated even if it is not required.
• obfuscated: a simple function written in file can be clearly laid out, with comments, help, input checking. Also makes it easy to create relevant unit tests etc.
• debugging: much easier with a function file!
• easier to avoid numeric issues caused by 0/NaN/Inf (sometimes occurs if using the multiplication approach to selecting data, unlike actual indexing).

[Geoffrey Rivers]

Ok, I'm happy to try to do this. However, the function (fh) needs to be used in the following manner, as part of the Delaunay algorithm for triangulation:

hbars=feval(fh,(p(bars(:,1),:)+p(bars(:,2),:))/2);

Where Bars is an array of node pairs based on the Triangle mesh list that is generated from p, the point list. Main point: I need to be able to evaluate fh for each p based on bar's contents. And frankly, I'm not clear on how to replace this with the call, or how to write the function to capture the passed variables or return the values.

This is all working fine as is if fh is a simple equation based on p, or is constant. It's just this case-based evaluation that I don't quite understand how to set up the calls or how to manage the data structures.

Would it be something like this?:

In main function:

hbars=fhassess((p(bars(:,1),:)+p(bars(:,2),:))/2, fd, fd2, scaledMat1delta, scaledMat2delta));

And then in fhassess:

function [hassign]=fhassess(pbars, fd, fd2, scaledMat1delta, scaledMat2delta)
holder = zeros(size(pbars));
metricA = feval(fd,pbars);
metricB = feval(fd2,pbars);
ixA = metricA<0;
ixB = metricB<0;
ixC = metricA >=0 && metricB >= 0;
holder(ixA) = scaledMat1delta;
holder(ixB) = scaledMat2delta;
holder(ixC) = (scaledMat1delta+scaledMat1delta)/2;
hassign = holder;
end

Sorry for the basic questions here, I'm normally an experimentalist, I code maybe once per year.

[Torsten]

ixA = metricA<0;
ixB = metricB<0;
ixC = metricA >=0 && metricB >= 0;

What about the case metricA<0 and metricB<0 ? Then you have an overlap in ixA and ixB.

And I wouldn't use feval if fd and fd2 are function handles. I'd simply use

metricA = fd(pbars);
metricB = fd2(pbars);

[Paul]

@Geoffrey Rivers

You may be interested in this blog post that shows how to implement if-else construct via anonymous functions.

[Geoffrey Rivers]

@Torsten Good question regarding the cases. Zone A and Zone B are mutually exclusive (*makes anxious sound* ... in my current workflow... i suppose when I generalise later they won't necessarily be... hmm :/ ) because they're derived from two signed distance functions that are measuring off the same geometry (inverted copies of the same dataset being eroded away from each other). So there's on;y a third case because there's a zone between them where they've both been eroded away from. But yes, I need to give some more thought to make sure I don't leave a logic loophole in this step, some way to catch an error if someone feeds this a zone where both can be < 0, which would be nonsense physically.

Oh, I can just do fd(pbar)? Oh, ok, that's not too bad then, thank you.

@Paul Oooooooohhhhhhhhh ok, ok yes this is very informative, thank you. I'm going to have to read through this a very times.

Thank you all for your help so far. I'm going to look at a few implementations of this and get it working, once I fix whatever it was I broke yesterday: I tried an inverted case yesterday, and discovered apparently I forgot to put a 'wall' around my Delaunay triangulation, so all my entries in p are leaking outside my XY array's dimensions and accelerating away at the speed of sound. Trying to fix that has broken everything else. ...Programming!