I'm pleased to introduce Peter Webb as this week's guest blogger. Peter has worked on the MATLAB Compiler since 1995. This is the first in a series of posts about how to make your M-files compile more easily.

Contents

The MATLAB Compiler and the deployment tools create deployable applications by "freezing" a collection of MATLAB functions and state information into a portable package. Deployed applications run using the MATLAB Component Runtime (MCR), which "thaws" the functions and state back into executable form. The Compiler creates either standalone executables or shared libraries (DLLs). Many of the differences between MATLAB and the MCR arise because a single user-written application may make use of multiple Compiler-generated shared libraries; in this case the shared libraries share the MCR's global data.

The MCR differs from MATLAB in five ways:

Because the MCR differs from MATLAB, M-Files written to run using the MCR may sometimes differ from M-Files that run only in MATLAB. Generally speaking, the differences between MATLAB and the MCR fall into four categories:

This article describes some general guidelines and then focuses on the problem of finding functions at compile time; I will address the other three types of differences in subsequent posts.

General Guidelines

The process of compiling and deploying an application is more efficient if the application follows these guidelines.

    global data;
    y = data + x;
    data = x;

    fp = fopen('Compute.m', 'w');
    fprintf(fp, 'function x = Compute(y)\n');
    fprintf(fp, '    x = y + 17;\n');
    fclose(fp);
    rehash;
    z = Compute(10);

    cd mathFcns
    z = add(x, y);
    cd ../stringFcns
    s = add(s1, s2);

    if ~isdeployed
        print
    else
        deployprint
    end

    switch (userInput)
        % There's no editor available to deployed applications;
        % print a warning message.
        case 'EditMFile'
            if isdeployed
                warning('MyApp:NonDeployable', ...
                    ['The editor is not available ', ...
                    'in compiled applications.']);
            else
                edit(mfile);
            end
        % Deployed applications can't call DOC. But they can
        %  redirect a help query to The MathWorks help site.
        case 'Help'
            if ~isdeployed
                doc(mfile);
            else
                web('http://www.mathworks.com/support');
            end
    end

Finding Functions at Compile Time

In order to produce a deployable package that is smaller than a MATLAB installation, the MATLAB Compiler attempts to narrow the set of deployed functions to only those required by the application. The Compiler calls the MATLAB function depfun to compute the dependent functions of each M-File it compiles. This process differs from the mechanism by which MATLAB determines which functions to call.

MATLAB searches for the best match for a function at runtime, when the types of the inputs and the contents of the MATLAB path are known precisely. depfun analyzes M-files at compile-time, when much less information is available; as a result it may omit functions that your application requires. Functions may be omitted because they are non-deployable or invoked only as callbacks or by eval.

Non-deployable Functions

Licensing restrictions prevent most design time functions from being compiled or deployed. A design time function is one that changes or augments the basic functionality of a program. Design time functions include the MATLAB editor, most toolbox GUI-based tools, like the Image Processing Toolbox function imtool, and functions that create M-Files, like the Fuzzy Logic Toolbox function anfis. This code will compile, but won't run:

    [fis,error,stepsize] = anfis(trainingData);

Attempting to execute non-deployable functions generates MATLAB:UndefinedFunction errors at runtime. Check your application's code to make sure it uses only deployable functions. mcc generates a list of excluded functions into a file called mccExcludedFunctions.log. Search this file for the string "because of compilability rules" to find the non-deployable functions depfun excluded from compilation. Because of the conservative nature of depfun's compile-time analysis, mccExcludedFiles.log may be quite large, and may list functions that your application does not require.

Detecting Dependencies on Callback Functions

The MATLAB Compiler's dependency analysis statically analyzes M-file functions to determine what other M-files they depend on. This analysis cannot determine the types of any variables or examine the contents of any strings. As a result, callback functions, which are often mentioned only in strings, can be overlooked; the Compiler will produce an application that appears to work but which will generate runtime errors when it attempts to invoke the missing callbacks.

For example:

set(gca, 'ButtonDownFcn', 'LogButtonPress');

If LogButtonPress is not included in the deployed package, this program will fail. Include LogButtonPress by explicitly specifying it on the mcc command line, or inserting a %#function pragma in the M-file that uses LogButtonPress.

Specifying a callback using mcc:

mcc [other options] -a dir1/dir2/LogButtonPress.m

Using a %#function pragma

%#function LogButtonPress

Including Functions Called by eval

MATLAB's eval function poses the same kind of problem for the Compiler's dependency analysis that callbacks do: the Compiler cannot examine text strings to determine if they contain function calls. Applications that use eval fail at runtime if the functions invoked by eval are not included in the deployed package.

For example, this code that uses eval to create dynamically named MATLAB variables requires the function getdata:

for n = 1:count
   eval(sprintf('x%d = getdata(%d);', n, n));
end

The solutions to problems created by eval are the same as those for callbacks: use the %#function pragma, or explicitly deploy the function in question by listing it as an argument to your mcc command.

Keep Reading

I plan to address issues and workarounds surrounding runtime path management, non-supported functions and differences in function behavior in future posts. In the meantime, you can refer to the documentation for the MATLAB Compiler or post follow-up questions here.

Get the MATLAB code

Published with MATLAB® 7.6

I'd like to introduce a new guest blogger - John D'Errico - an applied mathematician, now retired from Eastman Kodak, where he used MATLAB for over 20 years. Since then, MATLAB is still in his blood, so you will often find him answering questions on the newsgroup and writing new utilities to add to MATLAB Central.

Contents

I'll assume you have some data points through which you wish to pass a curve, interpolating your data. (Initially, I will only talk about problems with one independent variable.) In these coming blogs, I'll try to show some ways to do exactly this, i.e., find a curve that passes through your data. Along the way I'll try to give some pointers on curve fitting, interpolation, modeling, approximation, etc.

Polynomial regression

A valid question for some to ask is why start out with a discussion about polynomial regression , when we really wanted to talk about interpolation. Many people mistake the ideas of interpolation with the approximation produced by a regression model, calling both of these things interpolation. So I'm starting out with some discussion about what interpolation is not. Plus, I want to assure an understanding of polynomials, since many of the tools for interpolation are polynomial based in some way.

Let us start by creating some data. An exponential is a good place to start, a simple curve shape that is easy to fit.

x = -1:.1:1;
y = exp(x);

Plot your data

It is always a good idea to plot your data. In fact, I'll suggest that you should plot everything. Plots are useful, since your eye and brain are splendid at things like pattern recognition. Only you know your data, as the scientist, engineer, or analyst. You will always benefit if you can employ your knowledge of a system as part of the modeling process.

Next, always think in advance about your goals for any model.

For example, I once had a problem where I knew I had some significant noise in our process, but I chose not to do any smoothing anyway. Any such smoothing would also have smoothed out some potentially important features of our process. Since I could survive with the noise in my interpolant, I chose the lesser evil.

Always know your goals for any such task.

plot(x,y,'bo')
xlabel X
ylabel Y
grid on
title 'Exponential data'

This is a nice, well-behaved function. It is of the form y = f(x). I'll pretend for the moment that I have no idea what was in the underlying functional relationship.

One thing I learned in some long past calculus course was that a Taylor series will provide an approximation for many functions. Polynomials are useful things. They are simple to use, simple to build, simple to work with. And a truncated Taylor series is basically a nice, simple polynomial.

So we will start with a linear polynomial approximation for this curve, built using polyfit . This is a utility provided in MATLAB to estimate a polynomial model using linear regression techniques. We could also use many other tools to build our polynomial model, but polyfit is a useful one, and easy to use.

A linear, or first degree polynomial (many use the words "order" and "degree" interchangeably), might be written mathematically as y(x) = a1*x + a2. In MATLAB we will merely store the coefficients, as a vector [a1,a0]. Note that a polynomial in MATLAB has it's coefficients stored with the highest order term first.

P1 = polyfit(x,y,1)

P1 =
    1.1140    1.1937

We can evaluate the polynomial with polyval.

yhat = polyval(P1,x);

plot(x,y,'bo',x,yhat,'r-')
xlabel X
ylabel Y
grid on
title 'Linear polynomial fit'

Look at the residuals

In fact, I'll claim the relationship we are modeling is not terribly well represented by a linear model. Depending on your needs for this model, you might have decided differently.

When you build a regression model, look at the residuals.

ALWAYS PLOT EVERYTHING!

Plot your residuals. In general, some good ways to plot the residuals are versus

Look at your plots. Think about what you see there. Compare it to your expectations.

Since this data was very simply generated, I'll dispense with some of those plots for brevity. Note that the residuals for this linear fit look vaguely like a quadratic polynomial. This is often the case when there is lack of fit in a polynomial. That lack of fit often looks like the first term we truncated from the Taylor series.

res = y - yhat;
plot(x,res,'bo')
xlabel X
ylabel Residuals
grid on
title 'Residuals for the linear fit'

If the residuals looked vaguely parabolic in shape, then it might make sense to use a second order (quadratic) polynomial for our fit.

P2 = polyfit(x,y,2)
yhat = polyval(P2,x);
plot(x,y,'bo',x,yhat,'r-')
xlabel X
ylabel Y
grid on
title 'Quadratic polynomial fit'

P2 =
    0.5402    1.1140    0.9956

Note that the residuals here look vaguely like a cubic polynomial, although they are much smaller in magnitude than the previous fit. Again, at each step as we increase the order of the model, the residuals will often tend to look much like a polynomial of the next higher order.

res = y - yhat;
plot(x,res,'bo')
xlabel X
ylabel Residuals
grid on
title 'Residuals for the quadratic fit'

Higher order polynomials - are more terms always better?

Polynomial modeling with polyfit is indeed simple and easy to do. In fact, we might decide to just use higher and higher order polynomials, always chasing the term we truncated in the previous model. After all, if a quadratic model is better than a linear one, then why not go to a cubic? Ten terms must be better than nine. Look at a tenth order model.

P10 = polyfit(x,y,10);
disp(P10(1:6))
disp(P10(7:11))

    0.0000    0.0000    0.0000    0.0002    0.0014    0.0083
    0.0417    0.1667    0.5000    1.0000    1.0000

We can write the Maclaurin series representation for the exponential function as

We can compare P10 to the coefficients of the known series. How well did we do in the fit?

series = 1./factorial(10:-1:0);
disp(series(1:6))
disp(series(7:11))

    0.0000    0.0000    0.0000    0.0002    0.0014    0.0083
    0.0417    0.1667    0.5000    1.0000    1.0000

Note that the higher order coefficients deviate somewhat from the known series, although the lower order terms appear to be quite accurate.

yhat = polyval(P10,x);
plot(x,y,'bo',x,yhat,'r-')
xlabel X
ylabel Y
grid on
title 'Tenth order polynomial fit'

The residuals oscillate tightly around zero. In fact, they are so small that this last polynomial begins to approach a true interpolating polynomial. Perhaps if we just add a few more terms we may get there. The numerical issues of floating point arithmetic will often preclude true interpolation down to the least significant bit anyway.

res = y - yhat;
plot(x,res,'bo')
xlabel X
ylabel Residuals
grid on
title 'Residuals for the tenth order fit'

What do you do with interpolation?

I'll start talking about true interpolation in my next blog. But remember that interpolation is different from the approximations provided by polyfit or any other regression modeling tool.

Until that time, please give me your comments on this blog, or ideas for future blog topics on interpolation or modeling in general. Do you have some interesting applications of interpolation? Some interesting problems?

Get the MATLAB code

Published with MATLAB® 7.6

I recently heard Greg Wilson speak about a topic dear to him, Beautiful Code. He and Andy Oran edited a book on this topic in which programmers explain their thought processes. The essay Writing Programs for "The Book", by Brian Hayes, really resonates for me, in part because of the actual problem he discusses, in part because of the eureka moment he describes. Brian wrote his code in Lisp. I will write similar code in MATLAB to illustrate the points.

Contents

Problem Description

Given 3 points in a plane, find out if they lie on a single line, i.e., find out if the points are collinear. Here are the first set of points we'll use for trying out algorithms.

p1 = [1 1];
p2 = [3.5 3.5];
p3 = [-7.2 -7.2];

Method 1 - See if Third Point Obeys Same Slope

In this variant, we compute the straight lines connected 2 pairs of points, and then check that the third point fits the same line.

collinear1(p1,p2,p3)

ans =
     1

dbtype collinear1

1     function tf = collinear1(p1,p2,p3)
2     m = slope(p1,p2);
3     b = intercept(p1,p2);
4     tf = (m*p3(1)+b) == p3(2);
5     end
6     function m = slope(p1,p2)
7     m = (p2(2)-p1(2))/(p2(1)-p1(1));
8     end
9     function b = intercept(p1,p2)
10    m = slope(p1,p2);
11    b = -p1(2)+m*p1(1);
12    end

Method 2 - Method 1 Modified

There's a flaw in method one in the case where the points are collinear and lie on a vertical line, i.e., where there are repeated x values. Then the denominator we calculate will have the differences of two x values that are the same, resulting in division by 0 and a non-finite slope. First, let's see what happens with the first algorithm and vertically aligned points.

q1 = [0 1];
q2 = [0 3.5];
q3 = [0 -7.2];
collinear1(q1,q2,q3)

ans =
     0

The first algorithm says these points are not collinear. Time to patch the algorithm.

collinear2(q1,q2,q3)

ans =
     1

dbtype collinear2

1     function tf = collinear2(p1,p2,p3)
2     m = slope(p1,p2);
3     b = intercept(p1,p2);
4     % put in check for m being finite in case
5     % points are all on y axis
6     if ~isfinite(m)
7         if (p3(1)==p1(1))
8             tf = true;
9         else
10            tf = false;
11        end
12    else
13        tf = (m*p3(1)+b) == p3(2);
14    end
15    end
16    function m = slope(p1,p2)
17    m = (p2(2)-p1(2))/(p2(1)-p1(1));
18    end
19    function b = intercept(p1,p2)
20    m = slope(p1,p2);
21    b = -p1(2)+m*p1(1);
22    end

Methods 3 and 4 - Compare Two Slopes Only

The next thing Brian noticed was he really didn't need to see if the 3rd point fit the same line. What he really needed to determine is if two lines had the same slope. If so, all three lines have the same slope (convince yourself that this is true, at least for Euclidean geometry).

Let's try this algorithm on our two sets of collinear points.

collinear3(p1,p2,p3)
collinear3(q1,q2,q3)

ans =
     1
ans =
     0

dbtype collinear3

1     function tf = collinear3(p1,p2,p3)
2     m1 = slope(p1,p2);
3     m2 = slope(p1,p3);
4     tf = isequal(m1,m2);
5     end
6     function m = slope(p1,p2)
7     q = p2-p1;
8     m = q(2)/q(1);
9     end

Once again we need to modify the code for infinite slopes.

collinear4(q1,q2,q3)

ans =
     1

dbtype collinear4

1     function tf = collinear3(p1,p2,p3)
2     m1 = slope(p1,p2);
3     m2 = slope(p1,p3);
4     if isfinite(m1)
5         tf = isequal(m1,m2);
6     else
7        if ~isfinite(m2)
8            tf = true;
9        else
10           tf = false;
11       end
12    end
13    end
14    function m = slope(p1,p2)
15    q = p2-p1;
16    m = q(2)/q(1);
17    end

Method 5 and 6 - Summing Lengths of Sides of Triangle

The next idea that Brian explores is exploiting the relationship between the sides of the triangle that 3 points define. The longest leg is shorter than the sum of the lengths of the other 2 sides, except when the points are collinear. Then the length of the longest side is equal to the sum of the 2 other lengths. Let's try it out.

collinear5(p1,p2,p3)
collinear5(q1,q2,q3)

ans =
     0
ans =
     1

Let's use the point values that Brian uses.

bh1 = [0 1];
bh2 = [14 5];
bh3 = [8 4];
collinear5(bh1,bh2,bh3);

dbtype collinear5

1     function tf = collinear5(p1,p2,p3)
2     side(3) = norm(p2-p1);
3     side(2) = norm(p3-p1);
4     side(1) = norm(p3-p2);
5     lengths = sort(side,'descend');
6     tf = isequal(lengths(1), ...
7         (lengths(2)+lengths(3)));
8     end

So far, so good. Points that are not collinear don't claim to be. But the first set was and the answer came out negative. Let's try another batch of points to be sure.

bh4 = [0 0];
bh5 = [3 3];
bh6 = [5 5];
collinear5(bh1,bh2,bh3)
collinear5(bh1*1e5,bh2*1e5,bh3*1e5)

ans =
     0
ans =
     0

Uh-oh! What's happened? Looking back at collinear5, we see we are no longer doing just addition, subtraction, multiplication, and division. Now we're computing the norm which involves computing a sqrt. As a result, we have more opportunity for numerical roundoff issues, even if the points are confined to be rational numbers (including integers).

collinear6(bh1,bh2,bh3)
collinear6(bh1*1e5,bh2*1e5,bh3*1e5)

ans =
     0
ans =
     0

Version collinear6 no longer does an exact equality comparison but looks for differences on the order of eps for the appropriate lengths.

dbtype collinear6

1     function tf = collinear6(p1,p2,p3)
2     side(3) = norm(p2-p1);
3     side(2) = norm(p3-p1);
4     side(1) = norm(p3-p2);
5     lengths = sort(side,'descend');
6     tf = ...
7         abs(lengths(1)-(lengths(2)+lengths(3))) ...
8                              <= eps(lengths(1));
9     end

Method 7 - Eureka, Compute Triangle Area!

Finally, as Brian describes it, he realized there was a completely different way to see if 3 points were collinear, and that was to compute the area of the triangle they define.

collinear7(bh1,bh2,bh3)
collinear7(bh1*1e5,bh2*1e5,bh3*1e5)

ans =
     0
ans =
     0

To compute the area, you could compute base * height / 2, but again this would involve square roots. The trig approach would involve transcendental functions. Another way to think of this is to view any 2 pairs of the points as defining a parallelogram, and the triangle area is half that area. To compute the area, treat the sides as vectors which you translate to the origin. Then the area computation is straight-forward, boiling down to something proportional to the determinant If the result is 0, the points are collinear.

dbtype collinear7

1     function tf = collinear7(p1,p2,p3)
2     mat = [p1(1)-p3(1) p1(2)-p3(2); ...
3            p2(1)-p3(1) p2(2)-p3(2)];
4     tf = det(mat) == 0;

This algorithm requires an equality test following 6 arithmetic operations. Pretty simple, very elegant, much less prone to numerical issues than other algorithms. Much nicer than all of the other code above, no accounting for singular behavior necessary.

Note: Added 9 June 2008

If you read the comments for this blog, you will find an even more satisfactory solution relying on svd. svd (and the related rank function) have better numerical attributes than det.

What's Your Recent Aha Moment?

Have you had a similar epiphany recently? We'd love to hear about it here.

Get the MATLAB code

Published with MATLAB® 7.6

In an earlier post I discussed the issue of floating point accuracy in calculations. I have another example to share today that has been lurking in the Mapping Toolbox product. The issue is with the function wrapTo180.

Contents

Work-around

There is a work-around for those of you using Mapping Toolbox with releases R2007b and R2008a.

Numerical Issue

In R2007b and R2008a, the implementation for wrapTo180 causes results from wrapTo180(lon) to differ very slightly from lon even for for certain values of lon in the interval [-180 180]. The differences are on the order of 2*eps(lon).

Test Case

The following should evaluate to true.

lon = 115.8323;

but doesn't.

isequal(lon,wrapTo180(lon))

ans =
     0

These quantities aren't equal because of the arithmetic being done: 115.8323 + 180 - 180 differs from 115.8323 by

(lon + 180 - 180) - lon

ans =
  2.8422e-014

which happens to be

2*eps(lon)

ans =
  2.8422e-014

The Code

Here's the code for wrapTo180

type wrapTo180

function lon = wrapTo180(lon)
%wrapTo180 Wrap angle in degrees to [-180 180]
%
%   lonWrapped = wrapTo180(LON) wraps angles in LON, in degrees, to the
%   interval [-180 180] such that 180 maps to 180 and -180 maps to -180.
%   (In general, odd, positive multiples of 180 map to 180 and odd,
%   negative multiples of 180 map to -180.)
%
%   See also wrapTo360, wrapTo2Pi, wrapToPi.

% Copyright 2007 The MathWorks, Inc.
% $Revision: 1.1.6.2 $  $Date: 2007/08/20 16:35:59 $

lon = wrapTo360(lon + 180) - 180;

Simple Fix

We can just skip the arithmetic operation for values,like 115.8323, that already are in the interval [-180 180]. A nice way to do this is with logical indexing and we can replace the line of code with this.

q = (lon < -180) | (180 < lon);
lon(q) = wrapTo360(lon(q) + 180) - 180;

Comments?

Have you ever had a similar issue, where a seemingly innocent code expression causes erroneous results? How have you found this out and solved it? I'd love to hear your thoughts here.

Get the MATLAB code

Published with MATLAB® 7.6

Using MATLAB, there are several ways to identify elements from an array for which you wish to perform some action. Depending on how you've chosen the elements, you may either have the list of elements to toss or the list if elements to retain. And you might not have much if any control yourself how the list gets presented to you since the list could be passed to you from another calculation. The lists might be indices, subscripts, or logical arrays (often referred to as masks). Let's look at how you might arrive at such a situation and see what the code looks like to perform one particular action, setting the desired element values to 0.

Contents

Note: I am not discussing efficiency in this article. It is highly dependent on the number of elements in the original array and how many will be retained or thrown out. This article focuses on specifying what to keep or replace.

General Setup

Here's the setup for this investigation. I will use a fixed matrix for all the methods and always end up with the same final output. The plan is to show you multiple ways to get the result, since different methods may be appropriate under different circumstances.

A = magic(17);
Result = A;
Result( A < mean(A(:)) ) = 0;

Let's look at the nonzero pattern of Result using spy.

spy(Result)

Method #1 - Using Subscripts of Keepers

Here's a list of the subscripts for the elements to keep unchanged.

[rA,cA] = find(A > (17^2)/2);

Next we convert the subscripts to indices.

Result1 = zeros(size(A));
indices = sub2ind(size(A),rA,cA);
Result1(indices) = A(indices);
isequal(Result, Result1)

ans =
     1

Why did I convert subscripts to indices? Let me illustrate with a very small example.

matrix = [ -1 1 0; 2 0 -2; 0 3 -3]
[rows,cols] = find(matrix==0)

matrix =
    -1     1     0
     2     0    -2
     0     3    -3
rows =
     3
     2
     1
cols =
     1
     2
     3

Now let's see what I get if I use the subscripts to address the selected elements:

matrix(rows,cols)

ans =
     0     3    -3
     2     0    -2
    -1     1     0

I get the full matrix back, even though I selected only 3 elements. This definitely surprised me when I first encountered this. What's happening?

MATLAB matches each row element with each column element. matrix([1 2 3],2) returns the elements from rows 1 through 3 in column 1.

matrix(1:3,2)

ans =
     1
     0
     3

To learn more about indexing in general, you might want to read these posts or search the MATLAB documentation.

Method #2 - Using Indices of Keepers

Here we used the single output form of find which returns indices instead of subscripts.

indA = find(A > (17^2)/2);
Result2 = zeros(size(A));
Result2(indA) = A(indA);
isequal(Result, Result2)

ans =
     1

Method #3 - Using Logical Keepers

We'll try keeping about half of the elements unchanged.

keepA = (A > (17^2)/2);
Result3 = zeros(size(A));
Result3(keepA) = A(keepA);
isequal(Result, Result3)

ans =
     1

keepA is a logical matrix the same size as A. I use logical indexing to populate Result3 with the chosen values from A.

Method #4 - Subscripts for Elements to Set to Zero

If instead we have a list of candidates to set to 0, we have an easier time since we don't need to start off with a matrix of zeros. Instead we start with a copy of A.

Result4 = A;
[rnotA,cnotA] = find(A <= (17^2)/2);

Convert indices to subscripts, as in method #1.

indices = sub2ind(size(A),rnotA,cnotA);

Now zero out the selected matrix elements.

Result4(indices) = 0;
isequal(Result, Result4)

ans =
     1

Method #5 - Indices for Elements to Set to Zero

If we're instead given indices, we simply skip the step of converting subscripts and follow similar logic to that in method #4.

Result5 = A;
indnotA = find(A <= (17^2)/2);
Result5(indnotA) = 0;
isequal(Result, Result5)

ans =
     1

Method #6 - Using Logical Arrays to Specify Zero Elements

Finally, if we have a mask for the values to set to 0, we simply use it to select and set elements.

Result6 = A;
keepnotA = (A <= (17^2)/2);
Result6(keepnotA) = 0;
isequal(Result, Result6)

ans =
     1

Which Method(s) Do You Prefer?

Which method or methods do you naturally find yourself using? Do you ever invert the logic of your algorithm to fit your way of thinking about addressing the data (the ins or the outs)? Please post your thoughts here. I look forward to seeing them.

Get the MATLAB code

Published with MATLAB® 7.6