Loren on the Art of MATLAB

Contents

Introduction

You saw in my previous blog that high order polynomials can have some problems. Why not go to the opposite extreme? Use a piecewise version of linear interpolation? I like to call it connect-the-dots, after the child's game of that name. This is really the simplest interpolation of all.

In MATLAB, given a list of points, sampled from some functional relationship in one dimension, how would we perform piecewise linear interpolation? There are really two steps.

For any point u, given a set of (x,y) pairs with a monotonic vector x (by monotonic, I mean that x(k) < x(k+1) ), first find the index k, such that

Second, perform the linear interpolation to predict the value of y at x=u, between the pair of points (x(k),y(k)) and (x(k+1),y(k+1)).

Each data point in the list of points becomes a point where the slope of the piecewise linear interpolant changes to a new value. However, the function is still continuous across those locations. So one might call these locations "knots" because at those points consecutive polynomial segments are tied together. Knots is a common term applied to splines for these locations; breaks is a common alternative name.

Create Some Data to Interpolate

x = linspace(0,2*pi,10)';
y = sin(x);
plot(x,y,'o')
title 'A simple trig function'
xlabel X
ylabel Y

Suppose I want to interpolate this function at some intermediate point, perhaps u == 1.3?

u = 1.3;

histc Solves the Binning Problem

The first step I describe above is what I call binning. histc solves that problem efficiently.

[k,k] = histc(u,x);
k

k =
     2

As an aside, look at the way I took the output from histc. Since I only need the second output from histc but not the first output, rather than clutter up my workspace with a variable that I did not need, the output style [k,k] returns only the information I need.

Next, it seems instructive to dive a little more deeply into binning, so let me offer a few alternatives to histc.

Binning - A Loop With An Explicit Test

Just a test inside a loop suffices.

for k = 1:(length(x)-1)
  if (x(k) <= u) && (u < x(k+1))
    break
  end
end
x
k

x =
            0
      0.69813
       1.3963
       2.0944
       2.7925
       3.4907
       4.1888
       4.8869
       5.5851
       6.2832
k =
     2

Binning - A Semi-vectorized Test

Do the binning with a single vectorized test. This works reasonably as long as I have only a scalar value for u, so I call this only a semi-vectorized solution.

k = sum(u>=x)

k =
     2

When I want to place multiple points into their respective bins, this sum fails.

There are other ways to place points in bins in MATLAB, including a sort, or with hash tables. You can find several different binning methods in bindex on the file exchange. It is useful mainly to those with older MATLAB releases, because histc became available with version 5.3 and later of MATLAB.

Fully Vectorized Binning

Next, I may, and often do, have a list of points to interpolate. This is a common event, where I wish to more finely resample a curve that is sampled only at some short list of points. This time, let me generate 1000 points at which to interpolate the sampled function.

u = linspace(x(1),x(end),1000)';

[k,k] = histc(u,x);

This next line handles the very last point. Recall the definition of a bin as

Note the strict inequality on the left. So that very last point (at x(end)) might be technically said to fall into the n|th bin when I have |n break points.

On the other hand, it is more convenient to put anything that lies exactly at the very last break point into bin (n-1).

By the way, any piecewise interpolation should worry about points that fall fully above or below the domain of the data. Will the code extrapolate or not? Should you extrapolate?

n = length(x);
k(k == n) = n - 1;

Interpolation as a Linear Combination

The final step is to interpolate between two points. Long ago, I recall from high school what was called a point-slope form for a line. If you know a pair of points that a line passes through, as (x(k),y(k)) and (x(k+1),y(k+1)), then the slope of the line is simple to compute. An equation for our line as a function of the parameter u is just:

I can also rewrite this in a way that I like, by defining a parameter t as

Then the interpolant is a simple linear combination of the function values at each end of the interval.

Do the Interpolation and Plot the Result

See how nicely this all worked, in a fully vectorized coding style.

t = (u - x(k))./(x(k+1) - x(k));
yu = (1-t).*y(k) + t.*y(k+1);

plot(x,y,'bo',u,yu,'r-')
xlabel X
ylabel Y
title 'The connect-the-dots interpolant'

Use interp1 Instead

It is always better to use a built-in tool to solve your problem than to do it yourself, so I might just as well have used interp1 to accomplish this interpolation.

yinterp1 = interp1(x,y,u,'linear');

plot(x,y,'bo',u,yinterp1,'r-')
xlabel X
ylabel Y
title 'The connect-the-dots interpolant, using interp1'

In my next blog I'll begin talking about piecewise cubic interpolants. Until then please tell me your ideas here. Are there some associated topics that I should cover?

Get the MATLAB code

Published with MATLAB® 7.6

Contents

The thing I like most about the MATLAB environment is the flexibility I have as a programmer to choose the function/API that I think is required for the task at hand. And what's neat is that if I don't find or like what comes "in-the-box", I can create my own.

The enhanced object-oriented programming capabilities in the R2008a release of MATLAB add some new tools to our programming toolbox. While I can still program with existing tools such as functions and scripts, here are some situations where creating a class makes good design sense.

Create a New Data Type

MATLAB, like most programming languages, has a set of primitive data types with a set of operations defined on them. However, there are situations when we may need to create our own data type. For example, I might want to use structure arrays to track employee information but want to restrict the addition of new fields or modification of the values associated with existing fields.

employee.age = 23;
employee.name = 'Jane';
employee.ID = 17;
employee(2) = struct('age',39,'name','Jerry','ID',45);
employee

employee = 
1x2 struct array with fields:
    age
    name
    ID

The problem I face is that the fields of structure arrays can be modified easily by assigning new values to them just like I can modify any variables that store primitive data types.

On the other hand, if the variable employee is an object, only functions I choose (as creator of the class) can modify the fields.

The ability to control access and modification of data is referred to as encapsulation. Designing a class facilitates encapsulation thereby improving the quality of my code since objects are robust to accidental modification.

At this point, I'd like to clarify the distinction between classes and objects. Classes are a mechanism to describe interfaces to the functionality with which users will interact whereas objects are concrete instances of a class with their own data and related operations. Classes are designed and objects are used.

Another reason I may want to create a new type is to extend or redefine operations on an existing type. For example, I may require integer arithmetic to wrap as opposed to saturate (which is what MATLAB does) on overflow. I may, in addition, require a sqrt function that works on integers which isn't implemented in MATLAB. I can address these requirements by subclassing from the appropriate integer class in MATLAB and overwriting the arithmetic operations such as +, -, .*, ./, .\ and .^ to wrap on overflow. I can also add a method to my new integer class that computes the sqrt for integer values. Subclassing facilitates code reuse thereby eliminating the need to create new capabilities from scratch.

NOTE: If you do choose to override an arithmetic operator like plus, be aware that you may have to overload other functions that may use plus to ensure they give the right answer.

Using MATLAB's external interface APIs and interoperability capabilities, I can connect to external systems be they external devices such as data acquisition tools or external libraries created using Java or Component Object Model(COM). A convenient way to interact with these systems is to design a class that manages data access and modification as well as provides sufficient control over the functionality in the external system that I'd like to expose in MATLAB. By designing a class, I will be able to add new behaviors as the libraries I interface with evolve without compromising compatibility. Also, if the external system is implicitly object-oriented then designing a class is the more natural way to present it in MATLAB.

An example in MATLAB of a class which interfaces with an external library is the timer functionality which uses Java's timer object (java.util.Timer).

% Create a timer object and set its
% execution mode to schedule timer events.
t = timer('TimerFcn','figure;imagesc(magic(5))','Period',1);
t.ExecutionMode = 'fixedrate';
% Set name of the timer
t.Name = 'MyTimer';
% Set the number of times to execute the callback function.
% In this example, the image of magic squares is drawn twice,
% once for each of the two events.
t.TasksToExecute = 2;
% Start the timer.
start(t)
% Stop and release the timer.
stop(t)
delete(t)

In this example, designing a class allowed me to reuse timer functionality available in Java while giving me the flexibility to present APIs in a manner that is familiar to users interacting with objects in MATLAB.

Note that when using new classes written in MATLAB 7.6, I don't need to explicitly call the delete method to remove an object from memory; MATLAB does this automatically when the object is no longer used. The ability for an object to clean up, via the destructor, when MATLAB determines that the object is no longer being used is just one of the many enhancements to the object-oriented capabilities in MATLAB 7.6.

Manage the Lifecycle of Variables

Generally MATLAB takes care of managing memory for me. For example, create a 50x50 matrix of doubles.

a = rand(50);
a = 5;

With the second assignment to a, MATLAB takes care of freeing the memory allocated in the previous statement without any user intervention. Whenever variables go out of scope, are modified, or are no longer used, MATLAB automatically releases the memory associated with these variables when these variables happen to be primitive types. However, if the variables are resources such as file handles, fonts, or represent other types, then MATLAB needs to transfer control to that object to release resources when the object is no longer used.

In some situations I may need to take additional action to explicitly manage the lifecycle of the variables that get created. By creating class destructors (the delete method), I can ensure that memory associated with a variable is freed up when the variable is no longer being used. For example, let's look at the memmapfile class in MATLAB. The memmapfile functionality creates an object that lets me map a file to memory whereby any operations that modify the object in memory have the effect of modifying the contents of the file on disk.

Here's how I can map a file named records.dat to a series of unsigned 32-bit integers and set every other value to zero via the data property in memory.

   m = memmapfile('records.dat', 'format', 'uint32', ...
                  'writable', true);
   m.data(1:2:end) = 0;

Since MATLAB does not explicitly allocate the memory for the data contained in the variable m, it cannot release the memory either (i.e., in this case, unmapping the mapped region in MATLAB's memory) when the memmapfile object goes out of scope or is modified via assignment, perhaps in a situation like this.

   m = rand(4);

If I exposed the memory mapping functionality directly, I would risk inadvertent memory abuse if users of this class do not programmatically unmap mapped regions. However, an object-oriented approach mitigates this situation since it provides the class control over the creation and destruction of objects.

Organize Code

Keywords in MATLAB 7.6 allow me to define my class and organize my code in a single M-file. By having properties, methods, and events defined in the same file, I can gain an insight into class data and operations relatively quickly. With methods, I also have the flexibility of separating the function definition from implementation. I do so by defining the methods in the file containing the class definition and implementing the methods in separate files. Depending on the number of methods my class has, keeping method implementation in separate files can reduce clutter in the class definition.

Summary

MATLAB 7.6 has significant enhancements to object-oriented capabilities that simplify the process of designing classes while providing sufficient functionality to handle your most complex class design. When programming, if you find yourself needing to do any one of the above, then you should consider designing a class. I'd love to hear when you've found creating classes to be useful. Let me know here.

Get the MATLAB code

Published with MATLAB® 7.6

Contents

One of the goals of the compilation process is to turn a flexible and easily modifiable MATLAB progam into a robust software component with consistiently predictable behavior: a compiled application should not be able to change which functions it calls or the way those functions work. As a result, certain aspects of the execution environment that are malleable in MATLAB become fixed or constant in a deployed application. Perhaps the most important of these is the mechanism by which MATLAB locates functions and data files: MATLAB's search paths.

A path consists of a list of directories that MATLAB searches to find files. MATLAB typically searches a path in list order, stopping at the first file that matches the search criteria. MATLAB uses paths for two reasons: to determine which functions to execute; and to locate data files in the file system.

MATLAB primarily interacts with four paths, two internal and two external: the MATLAB path; the MATLAB Java class path; the system path; and the system load library path.

Most familar, perhaps, is the MATLAB path, which affects how MATLAB determines which MATLAB function files to run and which MAT-files to open. Though the MATLAB path is the most visible, of nearly equal importance is the MATLAB Java class path, which specifies the directories in which MATLAB searches for Java functions. The two external paths list the locations which the operating system will look for executables and shared libraries, respectively. The Windows operating system uses the same path to look for executables and shared libraries, but most Unix systems search separate paths for these two types of files.

When writing an application intended for deployment via the MATLAB Compiler, keep in mind that MATLAB applications typically interact with paths in at least three independent ways:

After briefly listing some general guidelines for safely interacting with paths in MATLAB, I provide more detailed techniques to manage problems that can occur in each of these areas.

General Guidelines

A compiled application behaves more consistently if it:

These guidelines become more important when the compiled application executes in an envrionment that differs from the one in which it was developed: for example, when installed at a customer site.

In cases where it is convenient to keep path management commands in your M-files when they run in MATLAB, you can use the isdeployed function to skip over calls to path management functions in deployed applications. For example, this M-code adds my "beta"-quality M-files to the path when my applications run in MATLAB:

if ~isdeployed
    addpath /home/pwebb/mfiles/beta
end

MATLAB Commands that Directly Access the Path

All M-file and Java functions rely on the state of MATLAB's internal paths, but only a few MATLAB commands directly access or modify these paths. Most of these functions should be avoided in deployed applications because they have the potential to change the way the application works, either by changing the functions it calls or the data files it loads.

Path management commands often appear in startup or initialization code; when it comes time to deploy your application, you may have forgotten that the application relies on these path settings. Problematic path settings sometimes occur in these oft-forgotten files:

Protect all calls to path management functions in these files by enclosing them in an ~isdeployed if-block.

cd my/data/directory
fp = fopen('data.file', 'r');

fp = fopen(fullfile(ctfroot, 'my', 'data', 'directory', 'data.file'));

p = path;
save savedpath p

As further incentive to avoid these functions in deployed applications, please note the commands that change the Java class path (javaaddpath, javarmpath) also clear the values of global variables.

For more details on how to manage the MATLAB function namespace with MATLAB objects, please see the polynomial example in the MATLAB documentation.

Paths as Arguments to MATLAB Functions

Many MATLAB programs assume that the directory structure of the filesystem cannot change. When this assumption leads to the use of absolute paths in file names, e.g., load('c:/Work/MATLAB/data.mat') or an implicit dependency on MATLAB's ability to locate a file, e.g., load data.mat, it creates problems for deployed applications. A deployed application is very likely going to be installed on a machine that has neither the same directory structure nor the access to the same filesystem as the machine on which it was developed. The use of absolute paths causes errors because of differences between the development and the deployment file system structures. But relying entirely on relative paths establishes a fragile dependency on the location of the current directory. The solution: base all file locations on a known root.

In MATLAB, there are typically two file system roots that matter: the root of the external file system, typically / or c:/ on Unix and Windows systems respectively; and the MATLAB root, the directory in which MATLAB has been installed. In a deployed application there is an additional root: the CTF root, which is the location of the application's MATLAB content (M-files, data files, JAR files -- everything that the MATLAB Compiler put into the application's CTF file). In MATLAB ctfroot and matlabroot refer to the same directory: MATLAB's installation directory. In a deployed application, matlabroot refers to the directory in which the MATLAB Common Runtime (the MCR) is installed and ctfroot to the CTF root.

Deployed applications typically base their file names on the the CTF root. For example, assume an application that reads and writes MAT files. To ensure that the MATLAB Compiler includes the input data files in the CTF archive, specify the data files with mcc's -a (add) flag:

mcc fcn.m -a data.mat

In the application code, use which or ctfroot to construct a full path to the files being loaded or saved.

Use which when a known filename exists on the application's path:

[pathstr,name]=fileparts(which('data.mat'));
% Load x from the MAT file
load(name);
x = x + 17;
% Save x back to the MAT file
save(fullfile(pathstr, 'data.mat'), x);

Use ctfroot when creating a new file:

data = input('Enter data: ');
save(fullfile(ctfroot, 'more_data.mat'), data, '-append');

This use of ctfroot assumes that the application has permission to write to the file system on the deployment machine. This is not always true by default, especially with network installations. Work with your IT department to secure the appropriate file system access.

Relying on the Structure of External Paths

A MATLAB application that uses the system or loadlibrary functions or calls MEX files linked to non-MathWorks shared libraries has a dependency on the external system paths. The system command allows MATLAB programs to directly execute operating system commands and capture their output as text. This can be very convenient, but it of course requires that the command be available in the current file system. loadlibrary allows MATLAB programs to load shared libraries and call their exported APIs; it similarly requires that the library in question be available on the deployment machine's file system.

The techniques for managing external path dependencies are very similar to those used to manage internal path dependencies: use full paths, favor explicit dependencies over implicit ones, and include all required files in the deployed application.

Next Up: Non-supported functions

The next post in this series will explore mechanisms for identifying and managing calls to non-supported functions. 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

Contents

Setup

First I set up the data.

m = 1e5;
n = 100;
A = rand(m,n);

Code I'm Tempted to Write

And now here's the code I'm tempted to write, safely tucked inside the confines of a try statement.

try
    AZeroMean = A - mean(A);
catch ME
    disp(ME.message);
end

Matrix dimensions must agree.

As you can see, MATLAB does not allow binary operators to work on arrays with different sizes (except when one of the inputs is a scalar value). There are at least two ways to remedy this.

Timing repmat

Using the most excellent timeit utility that Steve Eddins posted to the file exchange, I now time the repmat calculations. First I create an anonymous function that does my calculation. Then I pass that function handle to timeit. timeit carefully warms up the function by running it enough so the times are not subject to first-time effects, figuring out how many times to run it to get meaningful results, and more.

frepmat = @() A - repmat(mean(A),size(A,1),1);
timeit(frepmat)

ans =
      0.30964

Indexing with ones

repmat uses a variety of techniques for replicating an array, depending on the details of what's being replicated. One technique is to index into the array with ones in the dimension to replicate. Here's an illustative example with a vector.

q = [17 pi 42 exp(1)];
q5 = repmat(q,5,1)

q5 =
           17       3.1416           42       2.7183
           17       3.1416           42       2.7183
           17       3.1416           42       2.7183
           17       3.1416           42       2.7183
           17       3.1416           42       2.7183

Timing Indexing

One thing I notice with the repmat solution is that I need to create the vector mean(A) for the function. I need to do the same thing without repmat and I want to be able to set up one function call for performing the calculation so I can use timeit. Since I can't index into the results of a function without assigning the output to a variable, I create an intermediate function meanones to help.

type meanones

function y = meanones(A)

mn = mean(A);
y = A - mn(ones(1,size(A,1)),:);

Now I'm ready to do the timing.

findex = @() meanones(A);
timeit(findex)

ans =
      0.31389

Timing bsxfun

Next see the timing calculation done using bsxfun.

fbsxfun = @() bsxfun(@minus,A,mean(A));
timeit(fbsxfun)

ans =
      0.20569

Punchline

In this example, bsxfun performs fastest. Now that you see bsxfun in action, can you think of uses for this function in your work? Let me know here.

Get the MATLAB code

Published with MATLAB® 7.6

Contents

Introduction

The other day I was reading a post about sound card spectral analysis from the excellent new iheartmatlab blog. Beyond being a great post on recording from sound cards, the author ran into an issue involving object lifetime that has tripped me up in the past too.

So, today I would like to talk about two types of objects in MATLAB, scoped objects and user-managed objects. Knowing which kind of object you are working with can avoid all sorts of confusion. Let's start with a small example.

Recording Audio With a Callback

Suppose you want to use the audiorecorder object to record some audio from your sound card. The audiorecorder has a property called TimerFcn, which allows you to specify a callback function. This callback function calls a function you specify while you are recording audio. Here's an example function.

function makeRecorder( timerCallback, timerPeriod )

% Create an audio recorder object
recObject = audiorecorder();

% Install a callback is called every 'timerPeriod' seconds
set( recObject, ...
    'TimerFcn', timerCallback, ...
    'TimerPeriod', timerPeriod );

% start recording
record( recObject );

end

Once recording has started, the function timerCallback is called every timerPeriod seconds. This can be handy for inspecting audio data as it is streaming in.

On the face of things the makeRecorder function looks like it is ready to run. But, there is one subtlety that causes this function not to work as expected. When the makeRecorder function ends, the recObject variable goes out of scope, is removed from the workspace, and cleans itself up. For an audiorecorder object, cleaning up means:

In this case the timerCallback is never called because when makeRecorder ends, recording immediately stops and the audiorecorder's timer is shut down.

The reason the object is cleaned up is because the lifetime of an audiorecorder object is bound by its scope. To keep the object alive you would need to copy the object's local variable out of the makeRecorder function.

function recObject = makeRecorderBetter( timerCallback, timerPeriod )

% Create an audio recorder object
recObject = audiorecorder();

% Install a callback is called every 'timerPeriod' seconds
set( recObject, ...
    'TimerFcn', timerCallback, ...
    'TimerPeriod', timerPeriod );

% start recording
record( recObject );

end

Now recObject will be in the workspace of the code that calls makeRecorder, which allows the audiorecorder to stay alive.

Scoped Objects vs. User-Managed Objects

The audiorecorder object is what I like to call a scoped object because it cleans up after itself when it goes out of scope. Scoped objects are the most common type of object in MATLAB. They behave just like regular scalar and matrix variables.

% Create and clear a scalar
z = sin(pi);
clear z;

% Create and clear a 3 x 3 matrix
M = magic(3);
clear M;

% Create and clear a timeseries object
t = timeseries();
clear t;

Clearing the timeseries object above removes the timeseries from memory completely, just like clearing a scalar or matrix variable.

The reason the makeRecorder example above can be confusing is that some objects in MATLAB clean themselves up when they go out of scope, and some don't. If you try the same type of thing with a timer object, the results are very different.

function makeTimer( timerCallback )

timerObject = timer('TimerFcn',timerCallback, ...
                    'Period', 0.25, ...
                    'ExecutionMode', 'fixedSpacing' );
start( timerObject );

end

In this example, the variable timerObject goes out of scope at the end of makeTimer, but the timer does not get cleaned up. So, the timer fires every 0.25 seconds and the callback function timerCallback is called appropriately. Let's test it out.

% Create a timer
makeTimer( @(obj, event) display([obj.Name ' fired']) );

% pause for a second to let the timer fire
pause(1);

timer-1 fired
timer-1 fired
timer-1 fired
timer-1 fired

You can see here that the timer fired 4 times in one second, even though I can no longer get to the timer object I created.

Timer objects are what I like to call user-managed objects, in that they don't clean up after themselves. You are responsible for cleaning up a timer by calling the delete method after you are done with the timer.

But, after makeTimer ends, the timerObject variable is no longer in the workspace, so you will need to find it using the timerfind function.

% find all outstanding timers
timers = timerfind();

% stop the timers you find
stop(timers);

% delete and clear the timers
delete(timers);
clear timers;

Why do User-Managed Objects Exist?

Using a user-managed object seems like extra work, but you get flexibility over the object's lifetime expressely because it can remain around regardless of a variable's scope. This can be very helpful in several situations, like when you are building GUIs. For example, you could create a timer object in the CreateFcn of a figure window which refreshes the figure every second, and then clean up the timer in the figure's DeleteFcn.

Also, if you are working with a user-managed object and would like to treat it as a scoped object, you can use onCleanup, which was introduced in R2008a. See Loren's previous post that covers onCleanup, Keeping Things Tidy.

Scoped or User-Managed?

MATLAB has a lot of different user-managed objects including timer, analoginput, and videoinput to name a few. Each of these has a corresponding find function to help you find objects that have been removed from the workspace but are still hanging around in memory.

If you are struggling to figure out if an object is scoped or user-managed, follow these guidelines:

Regardless of the guidelines above, it is good just to know that scoped objects and user-managed objects exist in MATLAB and that there are differences between the behavior of each.

Has this ever happened to you?

Have you ever run into issues with user-managed objects? Has an object in MATLAB ever been cleaned up when you weren't expecting it too? Share you experiences here.

Get the MATLAB code

Published with MATLAB® 7.6

Contents