Nick Trefethen has framed the original flip chart that I described in my Lake Arrowhead Coauthor Graph blog earlier this week and has it on the wall of his office at Oxford University. His colleage Nick Hale has taken these photos. (The faded name in the cell at the lower left just northwest of mine is Tony Chan.)

Let's repeat John Gilbert's social network analysis of this data.

load housegraph
r = symrcm(A(2:end,2:end));
prcm = [1 r+1];
drawit

Get the MATLAB code

Published with MATLAB® R2013b

Published with MATLAB® R2013b

Greg's pick this week is Dynamic Mask Field Changing From Selected Pulldown Menu by David Manegold.

Ever wanted to change how your Simulink block mask appears depending on certain parameter conditions of the mask? David has a nice example.

I selected this submission because it demonstrates an advanced Simulink development technique in a straight-forward manner.

Contents

Change Thy Appearance, and Quick!

Masks have the ability to change what parameters are available depending on which parameters have been selected. In this example, if the value of the "Selector" parameters is "A" then the set of parameters associated with "A" appears in the mask. But if you select "C" instead, then a new set of parameters is displayed. Also option "C" has a subset of parameters that can also change depending on what you choose for "C alt field type 3".

What Doest Thou Hide Behind the Mask?

The mask is a dialog box that appears when you double-click on a block in Simulink. It provides a means to enter parameters for a Simulink block or subsystem to affect its functionality. This means you can have two instances of a block or subsystem in the same model that exhibit different behavior depending on the parameter values entered in the mask. You can create your own masks for subsystems or blocks by right-clicking on the block and select

Mask -> Create mask...

For more information on creating masks see Create Mask Documentation.

Mask, Reveal Thy Magic!

David leverages two capabilities to enable this behavior in a block mask.

In the mask editor, if you select a parameter to edit, you can populate the "Dialog callback" field to execute MATLAB functions that get executed when the parameter is changed in the block mask.

Within the callback function ChangeMaskFields, David changes the visibility of specific parameters based on the value of the "Selector" parameter. He defaults the parameter visibilities to "off":

[M{:}]=deal('off');

except for the first parameter (which is the Selector parameter):

M{1}='on';

Then he enables specific parameters based on the value of the "Selector" parameter:

[M{i0:ie}]=deal('on');

Comments

If you would like to leave any comments regarding this post, please click here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

It may look like a dirty sock, but I promise that it’s not.

Paul Kassebaum, who recently told us about the Battlebots competition, made this with a 3D printer. One of the cool things about 3D printing is that you can start with a shape that exists only in an abstract mathematical sense, and then you can extrude it onto the palm of your hand.

So here’s a hint. The shape above is related to another shape that Paul printed.

Maybe you already figured that out. But what is the relationship between the two shapes? Over the next couple of posts, Paul is going to tell us the story of these two shapes: how they were generated as abstract objects, and how they were given physical form in the real world. ]]> Uncategorized http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/zosI0qyMrZ0/ It can be difficult to see a line that is drawn over an image. The line is often lost in the background colors. That is why cursors are colored as they are, so that they are visible on any background. I demonstrate some code that makes a line more visible with the same technique. function [...] http://blogs.mathworks.com/videos/?p=1123 Mon, 10 Jun 2013 19:12:57 +0000

function h = cursorLine(x,y,innerThickness, outerThickness)

if nargin == 2
    innerThickness = 2;
    outerThickness = 4;
end

h.thick = line(x,y);
h.thin  = line(x,y);

set(h.thick, 'color', [1 1 1]);
set(h.thin , 'color', [0 0 0]);

set(h.thick, 'linewidth', outerThickness);
set(h.thin , 'linewidth', innerThickness);

Twenty years ago, during the Householder Symposium at the Lake Arrowhead conference center , John Gilbert carried out one of the world's first computational social network analyses.

Contents

Householder XII

Twenty years ago, in June, 1993, the UCLA Lake Arrowhead Conference Center in the San Bernardino Mountains, 90 miles east of Los Angeles, was the site of the Householder XII Symposium on Numerical Linear Algebra, organized by Gene Golub and Tony Chan.

John Gilbert Remembers

John Gilbert is now a Professor at U. C. Santa Barbara. A couple of years before the Lake Arrowhead meeting, in 1990 and 1991, John, Rob Schreiber, and I had developed the initial implementation of sparse matrices in MATLAB. John's research activities now include sparse matrix techniques for graph and social network problems. He recently wrote to me in an email, "It was at Householder 1993 in Lake Arrowhead that we used MATLAB's sparse matrices to do the first computational social network analysis I ever saw."

John used to be at PARC, the Xerox Palo Alto Research Center. A file containing his work from 1993 is still archived at PARC, <ftp://ftp.parc.xerox.com/pub/gilbert/graph.html>. Most of this blog post is derived from that archived work.

Coauthor Matrix

Nick Trefethen posted a flip chart with Gene Golub's name in the center. He invited everyone present to add their name to the chart and draw lines connecting their name with the names of all their coauthors. The diagram grew denser throughout the week. At the end it was a graph with 104 vertices (or people) and 211 edges. We entered the names and coauthor connections into MATLAB, creating an adjacency matrix A. Let's retrieve the names and matrix as they are stored in the PARC archive.

load housegraph
size(A)

ans =
   104   104

The matrix A is 104-by-104 and symmetric. Elements A(i,j) and A(j,i) are both equal to one if i and j are the indices of coauthors and equal to zero otherwise. Most of the elements are zero because most pairs of attendees are not coauthors. So the matrix is sparse. In fact, the sparsity, the fraction of nonzero elements, is less that five percent.

format short
sparsity = nnz(A)/prod(size(A))

sparsity =
    0.0486

Here is a picture of the matrix, a spy plot that shows the location of the nonzeros, with Golub in the first column and the other authors in the fairly arbitrary order in which they appeared on Trefethen's flip chart.

spy(A);
title('Coauthor matrix');
snapnow

Most Prolific

As the adjacency matrix is loaded from the archive, the most prolific coauthor is already in the first column. It was not a surprise to those of us at the conference to find that the most prolific coauthor is Gene Golub.

m = find(sum(A) == max(sum(A)))
name(m,:)

m =
     1
ans =
Golub               

Reorder the Authors

We want to create a circular plot with Golub in the center, the other authors around the circumference, and edges connecting the coauthors. If we were to place the authors around the circumference in the order we retrieve them from the chart, the edges would cross the circle pretty much randomly. We want to rearrange the authors so that the coauthor connections are as close as possible to the circumference of the circle. This corresponds to a symmetric permutation of the matrix that minimizes, or at least reduces, its bandwidth.

Reverse Cuthill McKee

The Reverse Cuthill McKee algorithm is the only algorithm I know that is usually used backwards. In 1969 Elizabeth Cuthill and J. McKee described a heuristic for reordering the rows and columns of a matrix to reduce its bandwidth. In 1991 J. Alan George pointed out that reversing the Cuthill McKee ordering almost always leads to fewer arithmetic operations in Gaussian elimination. The algorithm applies to symmetric matrices, so MATLAB has a symrcm function, but no symcm function. Here we are reducing bandwidth; we are not doing any arithmetic, so the difference between Cuthill-McKee and Reverse-Cuthill-McKee is purely cosmetic.

Coauthor Graph

r = symrcm(A(2:end,2:end));
prcm = [1 r+1];
spy(A(prcm,prcm))
title('Coauthor matrix with reduced bandwidth');

Now we are were able to actually plot the coauthor graph. Fortunately, plotting text at arbitrary angles had just been provided in MATLAB.

drawit;
snapnow

My Connectivity

You will find my name at about 5 o'clock on the coauthor graph. There is no edge connecting me to Golub because Gene and I were (not yet) coauthors. But how many coauthors do we have in common? Taking the second power of the connectivity matrix gives the paths of length two.

A2 = A^2;
count = A2(Golub,Moler)

count =
   (1,1)        2

And who are those common coauthors?

twos = name( find ( A(:,Golub) .* A(:,Moler) ), : )

twos =
Wilkinson           
VanLoan             

Arrowhead

Finally, we find a minimum degree reordering of the matrix and discover its true nature as, of course, ....

pmd = amd(A);
spy(A(pmd,pmd))
title('The (Lake) Arrowhead Matrix');

Get the MATLAB code

Published with MATLAB® R2013b

Published with MATLAB® R2013b

Once more, I am pleased to introduce guest blogger Kai Gehrs. Kai has been a Software Engineer at MathWorks for the past five years mainly working on features for the Symbolic Math Toolbox. He has a background in mathematics and computer science and already contributed to my blog in the past.

Contents

In a Nutshell: What Is This Article About?

In the article ODEs, from Symbolic to Numeric Code, Loren demonstrated how you can

In this article, we consider a slightly different scenario. Suppose that we have defined an ODE. First, we try to solve it symbolically. But we find that the symbolic ODE solver cannot find a closed form solution (something which is likely to happen, because only particular classes of ODEs can be solved in closed symbolic form).

Since we cannot solve the ODE symbolically, we must switch to a numeric ODE solver. But numeric solvers require a particular form and type of arguments, which often differ from the arguments of a symbolic solver. Symbolic solvers require a scalar symbolic (= "textbook-like") representation of an ODE. Numeric solvers usually require them to be transformed to a coupled first-order system.

This article demonstrates the following workflow of transforming a symbolic representation of an ODE to a form accepted by the MATLAB numeric ODE solver ode45.

Defining a Differential Equation in Symbolic Form

We start with the second-order non-linear ODE with two given initial values:

$$y'' = (y - y^2) \cdot y' - y, \qquad y = y(t), \qquad y(0) = 2, \qquad y'(0) = 0.$$

We define it in a symbolic form using a symbolic function symfun y(t):

syms y(t);
ode = diff(y,2) == (y-y^2)*diff(y) - y;

When trying to compute a solution of this ODE using dsolve, we see that the symbolic ODE solver does not find an explicit closed form solution:

dsolve(ode)

Warning: Explicit solution could not be found. 
ans =
[ empty sym ]

Adding the initial values does not help:

Dy = diff(y); dsolve(ode,y(0) == 2,Dy(0) == 0)

Warning: Explicit solution could not be found. 
ans =
[ empty sym ]

So we must use a numeric ODE solver to be able to plot the solution. But before we can convert the symbolic form of the ODE to a function handle accepted by ode45, we must convert the scalar form of the ODE to a coupled first-order ODE system.

From Scalar ODE to Coupled First-Order System

The Symbolic Math Toolbox function odeToVectorField converts a scalar ODE to a first-order ODE system:

V = odeToVectorField(ode)

V =
                          Y[2]
 - (Y[1]^2 - Y[1])*Y[2] - Y[1]

For details about the algorithm used to convert a general n-th order scalar ODE to a first-order coupled ODE system, see the odeToVectorField documentation page.

The next step is to convert the system representation V of the ODE to a function handle accepted by ode45.

From Coupled First-Order System to Function Handle Representation

To convert the expression V to a MATLAB function handle we use matlabFunction:

F = matlabFunction(V,'vars',{'t','Y'})

F = 
    @(t,Y)[Y(2);-(Y(1).^2-Y(1)).*Y(2)-Y(1)]

We are almost done. We only need to call the numeric ODE solver ode45 for the function handle F, and then plot the result.

From Function Handle Representation to Numeric Solution

Now we solve the differential equation converted to the function handle F:

sol = ode45(F,[0 10],[2 0]);

Here, [0 10] lets us compute the numerical solution on the interval from 0 to 10. Another additional parameter, [2 0], corresponds to the initial values: $y(0) = 2$ and $y'(0) = 0$.

Before plotting the solution, we use linspace to create 100 points in the interval [0,10] and deval to evaluate the solution at each of these points:

x = linspace(0,10,100);
y = deval(sol,x,1);

Now we plot the solution:

plot(x,y)

Have You Tried This Worklfow?

Have you tried using odeToVectorField and matlabFunction or creating symbolic functions with symfun? Let me know here. If you find Kai's article interesting, take a look at his other posts on my blog:

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Jiro's pick this week is Plot (Big) by my co-worker Tucker McClure.

"Big data"

That seems to be the buzz lately. But when people say "big data", it could mean a lot of things.

This entry by Tucker is all about "visualizing" big data. This was a recommendation from another one of my co-workers, Adam.

Visualizing Big Datasets

One of the issues with visualizing big datasets is that you have so many points displayed on the screen, and exploring (zooming and panning) may become sluggish. It's nice that you have all of those points available, but your screen has only so many pixels, and you don't get additional information by putting more points than what the screen is able to display.

So what can you do? Well, you can choose to display fewer points, just enough to give you a sense of what the dataset looks like. But then, you'd like to be able to see the details if you zoom in. It's sort of like Google Maps, for instance; when you're zoomed out, you see a rough outline of the city and some major highways. When you zoom in, you start to see some of the local roads.

Can we do something like this with MATLAB plots? I created something like that for line plots. Doug even wrote a Pick on it. However, that was more of a proof-of-concept that worked for only 2D line plots. Tucker went all the way. He implemented this using MATLAB Classes, and it could be applied to other plot types. You could also call it with additional arguments that you might pass into these plotting commands. Here are some additional things I like about Tucker's entry:

Here's a quick animation of reduce_plot in action, compared to a regular plot. You can see the difference in responsiveness of the plot when I am trying to pan the plot.

Thanks for an awesome entry, Tucker. And thanks for the recommendation, Adam. Some swag going to both of you!

Comments

We'd like to hear from you if you do, or are interested in, big data visualization. What kind of visualization do you use? Do you use any special techniques to deal with your big data? Give this a try and let us know what you think here or leave a comment for Tucker.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

I've mentioned the R2013a release in two previous posts (15-May-2013 and 12-Mar-2013). Today I want to point out that R2013a is a pretty significant release in terms of new features related to image processing and computer vision.

Here's the quick summary of what's new in Image Processing Toolbox, Computer Vision System Toolbox, and Image Acquisition Toolbox.

Image Processing Toolbox version 8.2

Computer Vision System Toolbox version 5.2

Image Acquisition Toolbox version 4.5

I've written extensively here in the past about geometric transformations (the spatial transformations category includes 21 different posts over a three-year period). Now there is much more to say. The new geometric transformation functionality represents a ground-up redesign, based partially on customer response and feedback to the functionality (such as imtransform) added way back in 2001. I hope to recruit a developer on the Image Processing Toolbox team to post here soon about the changes.

In the meantime, here's an example from the documentation for the new function imsharpen.

a = imread('hestain.png');
imshow(a)
title('Original Image');

b = imsharpen(a);
imshow(b)
title('Sharpened Image');

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

This week I received a large model where I wanted to find all the Lookup Table blocks and modify their Index Search Method parameter. Since the model was quite large and I was not familiar with it, manually navigating to each block, opening its dialog and modifying it was not an option.

So my question is: How would you accomplish this task?

Here are two options I considered.

Using functions at the MATLAB command prompt

For a long time, when I needed to find or modify blocks in a large Simulink model, I used the command prompt and functions like find_system, hilite_system and set_param.

For the example described above, I can find the blocks using find_system and see how they are connected in the model using hilite_system:

To change the Index Search Method, I can use set_param on the previous results:

Using the Model Explorer

Recently, I started using the Model Explorer for this type of task and I think it is very efficient.

For the same task as above, I can easily search by block type. One thing I like is that I do not need to remember the exact string for the block type. The drop down looks at my model and offers me a list based on the blocks in my model.

In the list of results, the Path column is a hyperlink I can click to immediately see where the block is in the model.

To edit the Index Search Method for all the blocks found, I:

• Add a column for this property
• Select all the blocks using "shitf+Click"
• Set the value for one instance, and it applies to all selected blocks

Now it's your turn

What is your workflow to find and edit blocks in a large model? Let us know by leaving a comment here.

So, it's been just over a year since I last asked readers to nominate their favorite (formerly unrecognized) File Exchange submissions. (My, how time has flown!)

Since mid-April, when I issued that challenge, MATLABbers around the world have shared approximately 3100 new files--and that doesn't include files that were simply modified in that timeframe. That means that more than 22% of the 14000 files currently on the Exchange didn't exist a year ago. And there's some great new stuff out there. Heck, there's even a "flatulence simulator." (I think I'll send Maxim some swag just for having posted that!)

So...time to reissue the challenge, but this time with a twist:

Nominate your favorite new submission--let's say, with a submission date (not modification date) of April 12, 2012 or later. If we agree, and if we feature your nomination in this blog, we'll send both you and the file's author some cool MATLAB swag.

The same ground rules apply:

I think we'll make this an annual event. Thoughts?

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

A Golden Spiral is simulated by a continuously expanding sequence of golden rectangles and inscribed quarter circles.

Contents

Golden Rectangles

We begin with an animated .gif of an imitation Golden Spiral. You can see that removing a square from a golden rectangle leaves a smaller rectangle with the same shape. Connecting inscribed quarter circles produces a continuously expanding spiral. The aspect ratio of the rectangles is the golden ratio.

$$ \phi = \frac{1 + \sqrt{5}}{2} $$

A MATLAB function for generating these expanding golden rectangles and quarter circles is given at the end of this blog post, but this is not a true Golden Spiral.

True Golden Spiral

A logarithmic spiral is a curve given in polar coordinates by

$$ r = \alpha e^{\lambda \theta} $$

The angular coordinate $\theta$ must be multi-valued as the point circles around the origin multiple times. We get the particular logarithmic spiral known as the Golden Spiral by involving the golden ratio and setting

$$ \lambda = \frac{2}{\pi}{\ln{\phi}} $$

Then the radius is scaled powers of $\phi$, with integer powers as $\theta$ crosses the cartesian axes.

$$ r = \alpha \phi^{\frac{2}{\pi} \theta} $$

Here is the plot for $0 \le \theta \le 2 \pi$.

I do not see any obvious analytic way to specify the scale factor $\alpha$. By taking $\alpha$ = 1.48699214 the numeric values generated by the logarithmic Golden Spiral approach those generated by golden rectangles and inscribed quarter circles as the angle $\theta$ gets large.

golden_spiral.m

function golden_spiral
% GOLDEN_SPIRAL  Explosion of golden rectangles.
%    GOLDEN_SPIRAL  Constructs a continuously expanding sequence
%    of golden rectangles and inscribed quarter circles.

%   Copyright 2013 Cleve Moler
%   Copyright 2013 The MathWorks, Inc.

   % Initialize_variables

   % Golden ratio
   phi = (1+sqrt(5))/2;

   % Control speed of zoom
   n = 256;
   f = phi^(1/n);

   % Scaling
   a = 1;
   s = phi;
   t = 1/(phi+1);

   % Centers
   x = 0;
   y = 0;

   % A square
   us = [-1 1 1 -1 -1];
   vs = [-1 -1 1 1 -1];

   % Four quarter circles
   theta = 0:pi/20:pi/2;
   u1 = 2*cos(theta) - 1;
   v1 = 2*sin(theta) - 1;
   u2 = 2*cos(theta+pi/2) + 1;
   v2 = 2*sin(theta+pi/2) - 1;
   u3 = 2*cos(theta+pi) + 1;
   v3 = 2*sin(theta+pi) + 1;
   u4 = 2*cos(theta-pi/2) - 1;
   v4 = 2*sin(theta-pi/2) + 1;

   initialize_graphics

   % Loop

   k = 0;
   while get(klose,'value') == 0
      if mod(k,n) == 0
         power
         switch mod(k/n,4)
            case 0, right
            case 1, up
            case 2, left
            case 3, down
         end
      end
      zoom
      k = k+1;
   end
   pause(1)
   close(gcf)

% ------------------------------------

   function power
      a = s;
      s = phi*s;
      t = phi*t;
   end % power

% ------------------------------------

   function zoom
      axis(f*axis)
      drawnow
   end % zoom

% ------------------------------------

   function right
      x = x + s;
      y = y + t;
      line(x+a*us,y+a*vs,'color','black')
      line(x+a*u4,y+a*v4)
   end % right

% ------------------------------------

   function up
      y = y + s;
      x = x - t;
      line(x+a*us,y+a*vs,'color','black')
      line(x+a*u1,y+a*v1)
   end % up

% ------------------------------------

   function left
      x = x - s;
      y = y - t;
      line(x+a*us,y+a*vs,'color','black')
      line(x+a*u2,y+a*v2)
   end % left

% ------------------------------------

   function down
      y = y - s;
      x = x + t;
      line(x+a*us,y+a*vs,'color','black')
      line(x+a*u3,y+a*v3)
   end % down

% ------------------------------------

   function initialize_graphics
      clf reset
      set(gcf,'color','white','menubar','none','numbertitle','off', ...
          'name','Golden Spiral')
      shg
      axes('position',[0 0 1 1])
      axis(3.5*[-1 1 -1 1])
      axis square
      axis off
      line(us,vs,'color','black')
      line(u3,v3)
      klose = uicontrol('units','normal','position',[.04 .04 .12 .04], ...
         'style','toggle','string','close','vis','on');
      drawnow
   end % initialize graphics

end % golden_spiral

Get the MATLAB code

Published with MATLAB® R2013b

Published with MATLAB® R2013b

Idin's pick for this week is the Delta Sigma Toolbox by Richard Schreier.

Delta-Sigma (or sigma-delta) modulators are commonly found in electronic components such as analog-to-digital and digital-to-analog converters (ADCs and DACs), and increasingly in [fractional-N] frequency synthesizers (PLLs), and switch-mode power supplies. Generally speaking, a delta-sigma modulator produces a highly over-sampled binary signal at its output (0/1) which can be low-pass filtered to reproduce the input signal. Figure 1 below shows an example input & output signal. To read more about delta-sigma modulators, go here, here, and here. Richard’s textbook is also a good resource: "Understanding Delta-Sigma Data Converters" by Schreier and Temes (ISBN 0-471-46585-2).

There is a large body of literature dedicated to analysis and design of delta-sigma modulators, and the Delta Sigma Toolbox provides a great tool for analyzing these components in MATLAB. The toolbox contains a fairly complete list of well documented functions to construct, analyze, and simulate delta-sigma modulators of arbitrary order. The full documentation for the toolbox is provided in DSToolbox.pdf (part of the download).

To start using the toolbox, I highly recommend studying the provided example and demo files (dsdemo1-8.m and dsexample1-4.m). I like dsdemo2.m as a starting tutorial. It begins by creating a noise transfer function (NTF), and then simulates the modulator with an example sinusoidal input:

OSR = 32;
H = synthesizeNTF(5,OSR,1);
N = 8192;
fB = ceil(N/(2*OSR)); ftest=floor(2/3*fB);
u = 0.5*sin(2*pi*ftest/N*[0:N-1]);  % half-scale sine-wave input
v = simulateDSM(u,H);

The time domain output shows the sinusoidal input (red) and the binary output (green).

The spectrum of the output can also be easily computed using MATLAB's FFT function:

This figure looks like the superposition of two signals: the spike at ftest/N (0.0104 in this case) which represents the input sinusoid, and some high frequency "noise". The high-pass noise is in fact one of the desirable properties of delta-sigma converters; these converters shift the quantization noise to high frequencies, i.e., away from our signal/carrier of interest. It should be self-evident from Figure 2 that a simple low-pass filter of this output signal would yield the original sinusoidal input.

The next few lines in dsdemo2.m compute the SNR on the output signal, and also use Delta Sigma Toolbox functions to compute and display the expected (theoretical) response of the modulator (shown in pink blow).

The next section of the code computes the expected and simulated signal-to-quantization-noise ( SQNR), which is a figure of merit when analyzing analog-to-digital converters.

The remainder of this demo essentially repeats the same process, but for a band-pass signal. That is, instead of using a low frequency sinusoid as the input, it uses a carrier signal, which makes our output spectrum look as in Figure 5 (the carrier here is at Fs/8 or 0.125 on the normalized frequency axis).

I encourage you to download the toolbox and run through the demos. Even if you do not use the toolbox for any simulations of your own, these demos can serve as a great learning tool (along with freely available online resources mentioned above).

Other MathWorks resources:

If you are interested in delta-sigma modulators (and other mixed-signal components), consider downloading the Mixed-Signal Library for Simulink. Also, take a look at the MathWorks Mixed-Signal page, particularly the ADC section.

Delta-Sigma Toolbox Usage notes:

mex simulateDSM.c

But as noted in the comments on File Exchange, this fails on some systems. Simply defining the _STDC_ symbol should resolve the issue:

mex simulateDSM.c –D__STDC__

Note that you will need a C compiler for this step. Run mex –setup at the MATLAB command prompt to setup your C compiler.

addpath('<your_intall_path>\delsig')
savepath

Suggestions for improvements

Comments

As always, your thoughts and comments here are greatly appreciated.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

A long time ago, I covered the topic of the duality between command and function syntax.

Contents

What is Command-Function Duality?

For calling functions with only literal string input values, i.e., strings delimited by the usual single quote ('), you can avoid the overhead of using parentheses and comma separators to invoke the function, especially if you want either no output or only one output. Let me show you an example.

Checking Directory Contents

Suppose I want to see the MATLAB code files in my directory beginning with the letter 'q'. Here are three ways to do such a query.

First way:

dir q*.m

qichen31.m  quadvec.m   quantum.m   

Second Way:

dir('q*.m')

qichen31.m  quadvec.m   quantum.m   

Third Way:

out = dir('q*.m')

out = 
3x1 struct array with fields:
    name
    date
    bytes
    isdir
    datenum

Discussion of Results

From these three methods, you can see that the first two give the same output, a list of the names of matching files. When using the command syntax (the first way), you can't store the output in a named variable. In this case, the decision was to print something out in a useful format, but not to place the information into a variable for further programmatic use.

The second result is equivalent to the first. You can see the duality looking at these two statements. The first and second are equivalent. Simply replace the space and following string (or strings) into a comma-separated list of the same string(s) inside parentheses. This means that the statement myfun A B c is equivalent myfun('A','B','c').

And now for the third statement. Clearly it is different than the first two. First, you can see that something is returned in the output variable out. Second, you can see that out is not simply a list of names, but a struct containing several fields with information, including name, date, etc.

Default MATLAB Behavior

MATLAB functions often return at least one output, even if the user does not supply an output variable. If the function does return an output when the user does not specify one, the result goes into a variable named ans.

Overriding Default Behavior

It is possible to override default behavior for your function outputs. To do so, we take advantage of the function nargout. This function allows us to query how many output variables the function is called with. In the case of the function dir, you can imagine the logic of the code goes something like this:

With the command form for calling a MATLAB function, the value for nargout is 0. Despite this, some functions are designed to return a value, though without a specified output, ans is created or updated.

What About You?

Do you take advantage of command-function duality? In functions you create? At the MATLAB prompt or in application code your write? Let me know here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Why using Linear Analysis Points?

Those doing linear analysis with Simulink for a long time are probably familiar with the linmod function. With linmod, the inputs and outputs of the linearization are the root-level Inport and Outport blocks of your model. This means that your model needs to be built in a specific way to accommodate linear analysis.

With Simulink Control Design, you can use Linear Analysis Points to specify which portion of a model to linearize.

To see all the advantages of Simulink Control Design compared to linmod, I recommend going through this documentation page: Linearization Using Simulink Control Design Versus Simulink.

Specifying Linear Analysis Points in a model

Let's take this model as example.

To design a controller, the first thing I am interested in is linearizing the plant. For that, I can right click on the signals entering and exiting the plant, and mark them as Open-loop Input and Open-loop Output. For those already familiar with Simulink Control Design, notice the new improved menu in R2013a.

Once this is done, you will see markers in your model for the marked signals.

Specifying Multiple sets of Linear Analysis Points

When developing controllers for your model, you will very likely need to accomplish linear analysis multiple times for different parts of your model, open and closed loop, with different inputs and outputs. For that, you can click on Create new linearization I/Os... in the Linear Analysis Tool.

This will launch a window where you can define multiple sets of points for your different tasks.

Once you are done, you can select the set of points to be used for linearization.

Now it's your turn

For those familiar with the earlier version of this menu, I recommend looking at the R2013a release notes to find a table describing the mapping between the previous and current Linear Analysis Points menu. Hopefully, this new menu will help avoid confusion when specifying linear analysis points, especially with the open and closed loop concept.

If you prefer, all the above is also available from the command line with functions like linio, getlinio, and setlinio.

Give that a try and let us know what you think by leaving a comment here.

A few weeks ago, MathWorks hosted a Student Robot Challenge at our new office building in Cambridge. The event was a great success. The students had fun and learned plenty about modelling, simulation, control design, code generation, and teamwork. I wrote a blog about the event here:

Fun, learning, and drama at the MathWorks Student Robot Challenge

We’ve compiled a brand-new video with the highlights of the day, which explains more about the contest and shows the LEGO MINDSTORM NXT robots in action!

Jiro's pick this week is GPUBench by Ben Tordoff.

Some of you may know of the bench function. It allows you to benchmark your MATLAB on your machine and compare it against other machines. It performs some mathematical computation tests, as well as graphics tests. Note, as mentioned in the documentation, that bench is for comparing a particular version of MATLAB on different machines, not comparing different versions of MATLAB on a single machine.

GPUBench is GPU version of bench. With Parallel Computing Toolbox, you can perform MATLAB computations on NVIDIA CUDA GPUs with Compute Capability of 1.3 or greater. But your mileage may vary depending on the hardware you have. Sometimes, looking at just the Compute Capability may not be enough to see if one card is better than another for MATLAB computations. GPUBench will test your card and compare with other common GPU cards. Here's a sample table that gets generated from the benchmark. It reports in gigaFLOPS, so higher the number, the better the performance.

When I bought a MacBookPro last year and noticed that it had a supported GPU card, I ran the benchmark on it and got the following result. Note that this was done with an older version, so the list of cards is a bit old.

The "Host PC" refers to the CPU on my Mac, and "GeForce GT 650M" is the GPU that's equipped. We can see that the GPU seems to outperform the CPU on single-precision computations but not so much on double-precision computations.

Ben periodically updates the data files to include the newer GPU models that have come out. As of the writing of this post, the version includes benchmark data for K20.

Comments

Do you use GPUs, or are you interested in using GPUs for your work? Let us know about it here. If you do GPU computing in MATLAB, give this App a try and leave a comment for Ben.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

The first MathWorks general product release of the year, R2013a, shipped a couple of months ago. I've already mentioned it once here in my 12-Mar-2013 post about the new MATLAB unit test framework.

With each new release, I peruse the release notes for MATLAB to see what things I find particularly interesting. (This helps me remember which product features have actually been released, as opposed to still being in development. My memory needs all the help it can get.)

The first thing to note is the reappearance of the table of contents for navigating in the Help Browser and in the online Documentation Center. This is a direct result of helpful feedback we received from many of you about the R2012b release.

My favorite "make-it-go-faster-without-sacrificing-accuracy" people (the MATLAB Math Team, that is) have been busy again. People with computers based on Intel or AMD chips using the AVX instruction set should see their calls to fft speed up. Anybody running permute on 3-D or higher-dimensional arrays should also get a nice boost. I've done a lot of development work related to image and scientific format support, so I know that a fast permute can be pretty useful when reading image and scientific data. That's because most of these formats store array elements in the file in a different order than MATLAB uses in memory.

In the small-but-nice category, the MATLAB Math Team also simplified a common programming pattern in my own neck of the woods (image processing). Specifically, it's a bit easier to initial an array of 0s or 1s whose type is based on existing array. Here's an example to illustrate:

clear  % Let's start with a fresh workspace.
rgb = imread('peppers.png');
imshow(rgb)
title('Obligatory image screenshot')

Now I want a 100-by-100 matrix of 0s with the same data type as rgb.

A = zeros(100,100,'like',rgb); % Make a 100-by-100 matrix that's "like" rgb.
whos

  Name        Size                Bytes  Class    Attributes

  A         100x100               10000  uint8              
  rgb       384x512x3            589824  uint8              

My developer friend Tom Bryan really "likes" this (ahem) because it enables much easier solutions to some common programming tasks for users of Fixed-Point Designer.

I have occasionally done a little web scripting in MATLAB, so it's nice to see urlread and urlwrite get a little love. These functions can now handle basic authentication via the 'Authentication', 'Username', and 'Password' parameters.

Do you use a Mac? You can now write MPEG-4 H.264 files using VideoWriter (requires Mac OS 10.7 or later).

A couple of handy new string functions have appeared, strsplit and strjoin. Based on how often users have submitted their own versions to the MATLAB Central File Exchange, I'm sure these will be popular.

out = strsplit(pwd,'\')

out = 

    'B:'    'published'    '2013'

You can now do extrapolation with both scattered and gridded interpolation. For extrapolation with scattered interpolation, use the new scatteredInterpolant. Here's an example I lifted from the doc.

Query the interpolant at a single point outside the convex hull using nearest neighbor extrapolation.

Define a matrix of 200 random points.

P = -2.5 + 5*gallery('uniformdata',[200 2],0);

Sample an exponential function. These are the sample values for the interpolant.

x = P(:,1);
y = P(:,2);
v = x.*exp(-x.^2-y.^2);

Create the interpolant, specifying linear interpolation and nearest neighbor extrapolation.

F = scatteredInterpolant(P,v,'linear','nearest')

F = 

  scatteredInterpolant with properties:

                 Points: [200x2 double]
                 Values: [200x1 double]
                 Method: 'linear'
    ExtrapolationMethod: 'nearest'

Evaluate the interpolant outside the convex hull.

vq = F(3.0,-1.5)

vq =

    0.0031

Disable extrapolation and evaluate F at the same point.

F.ExtrapolationMethod = 'none';
vq = F(3.0,-1.5)

vq =

   NaN

I encourage you to wander over to the R2013a release notes for MATLAB or any other product that you use and see what's new that might be helpful to you.

There are also lots of new things in the image processing and computer vision worlds, of course. I'll look at those next time.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

In my previous blog post I reprinted the Cleve's Corner article from the 1995 issues of MATLAB News and Notes and SIAM News about the Pentium division bug. In today's post I would like to describe some of the effects that affair had on the emerging Internet, the MathWorks, and the Intel Corporation,

Contents

1994 Internet

In the fall of 1994 we did not have anything like the Web as we know it today. There was an Internet, but the connections and protocols were all text based. We did have email and file transfer. News groups, including comp.soft-sys.matlab and comp.sys.intel, existed. The Mosaic graphical browser had been available for just a few months, but most people, including me, had not yet started to use it. Internet Explorer was not yet on the scene and it would be four more years before we would hear of Google.

When I began to get email asking for information about the bug, I would start my reply with "If you have access to the Internet ..." and then give instructions on how to use the MathWorks FTP site. Imagine doing that today.

And, as far as I can remember, there was not yet any spam. All of the email, and all of the postings to the newsgroups, were legitimate. Maybe they weren't all worthwhile, but at least there wasn't yet any of the absolute junk we unfortunately see today.

Initial Activity

Prof. Thomas Nicely started all the activity by sending a memo about an error in the calculation of the reciprocal of a twin prime to the CompuServe network on October 30, 1994. This soon found its way to Terje Mathisen in Norway who confirmed a bug in floating point division on the new Intel Pentium processor and posted an announcement, along with a test program, to the comp.sys.intel newsgroup on November 3.

A couple of MathWorks customers saw the news and called our tech support, asking how the bug might affect MATLAB. At nearly the same time I heard about the bug from a mailing list covering floating point arithmetic.

On November 14 Tim Coe posted what turned out to be an instance of the worst possible error. The next day I made my first posting, to both comp.soft-sys.matlab and comp.sys.intel, summarizing what was known up to then. I said that as far as I was concerned the worst aspect of the situation was that we didn't know how bad it was and so we had to be concerned about it.

Intel's Response

Here is Intel's initial response, in the form of the FAXBACK document produced by their customer support system.

There has been a lot of communication recently on the Internet about a floating point flaw on the Pentium processor. For almost all users, this is not a problem.

Here are the facts. Intel detected a subtle flaw in the precision of the divide operation for the Pentium processor. For rare cases (one in nine billion divides), the precision of the result is reduced. Intel discovered this subtle flaw during on going testing after several trillions of floating point operations in our continuing testing of the Pentium processor. Intel immediately tested the most stringent technical applications that use the floating point unit over the course of months and we have been unable to detect any error. In fact, after extensive testing and shipping millions of Pentium processor-based systems there has only been one reported instance of this flaw affecting a user to our knowledge, In this case, a mathematician doing theoretical analysis of prime numbers and reciprocals saw reduced precision at the 9th place to the right of the decimal.

In fact, extensive engineering tests demonstrated that an average spreadsheet user could encounter this subtle flaw of reduced precision once in every 27,000 years of use. Based on these empirical observations and our extensive testing, the user of standard off-the-shelf software will not be impacted. If you have this kind of prime number generation or other complex mathematics, call 1 800 628-8686 (International) 916 356-3551). If you don't, you won't encounter any problems with your Pentium processor-based system. If ever in the life of the computer this becomes a problem, Intel will work with the customer to resolve the issue.

Needless to say, that kind of response did not satisfy activists on the net.

Internet Emboldened

Intel had only recently begun to appeal directly to consumers with its "Intel Inside" campaign. Intel's traditional customers were computer manufacturers, and they were not concerned with "subtle flaws" that had already been fixed in the latest stepping of the part.

But now posts criticizing Intel started to appear on the net. A small, but important and vocal portion of the customer base was protesting.

On November 22, two engineers at the Jet Propulsion Laboratory suggested to their purchasing department that the laboratory stop ordering computers with Pentium chips. Steve Young, a reporter with CNN, heard about JPL's decision, found my posting on the newsgroup, and called me. After talking to me for a few minutes, he sent a video crew to MathWorks in Massachusetts and interviewed me over the phone from wherever he was in California. That evening, CNN's Moneyline used Young's news about JPL and his interview with me to make the Pentium Division Bug mainstream news. Two days later -- it happened to be Thanksgiving -- stories appeared in The New York Times, The Boston Globe, The San Jose Mercury News, and elsewhere. Hundreds of articles followed in the next several weeks.

On November 27, Andy Grove, CEO of Intel, tried to post his view of the situation on comp.sys.intel. But I'm afraid he only made matters worse. For one thing, he didn't even have his own logon, so he had a colleague, Richard Wirt, post the message on his behalf. Many netters were suspicious that it was a fake. The general tone was similar to the FAXBACK I reproduced above, aggressively defensive. And the posting was too long. I have included the entire posting in the collection of Pentium bug documents available from MATLAB Central.

In the three weeks from Mathisen's first posting until the end of November the traffic on comp.sys.intel increased from a trickle to hundreds messages a day. Most of it was rants about Intel and rants about the rants. I had to stop reading it. Some of it was cross posted to comp.soft-sys.matlab, so the traffic was heavy there too for a while.

The net had been heard. It had gone from being a plaything for a few academic geeks to something with economic clout. The Internet was beginning to grow up.

Effect on MathWorks

In the fall of 1994 the MathWorks was celebrating its 10th anniversary. We had not yet moved to Apple Hill. We were still located on Prime Parkway, off of Speen Street. We had less than 250 employees.

Our product name, MATLAB, was known to our customers, but our company name, MathWorks, was not widely known. On November 23, we announced that we were preparing a version of MATLAB that could detect and correct the division bug. Our public relations firm issued a press release with the headline,

THE MATHWORKS DEVELOPS FIX FOR THE INTEL PENTIUM(tm) FLOATING POINT ERROR

So on the day the story appears in The New York Times, this message showed up in the fax machines of media outlets all over the country. It turned out to be a very successful press campaign.

MathWorks has been at the forefront of companies using the Internet. We were 73rd company to register as a ".com" URL. So with the Mosaic browser becoming available about this time, we put up one of the first company Home pages.

I collected as many of the documents associated with this affair as I could and made them available in a .zip file from an FTP site. We put a notice about that collection on our home page, even though most people learned about it from newsgroups and email. Here is that first MathWorks Home Page.

Reaction at IBM

On December 12, IBM issued a report entitled "Pentium Study" that claimed financial spreadsheet calculations did not produce floating point numbers with uniformly distributed random bit patterns, but rather frequently produced numbers with long strings of 1's in their binary representation. These are the numbers that are "at risk" for encountering the FDIV bug. As a result, they predicted division errors would be much more frequent than Intel claimed.

Citing their report, IBM announced they were halting shipment of Pentium-based personal computers. Intel's stock dropped 10 points within a few hours of the IBM press release.

Intel Response

A week after the IBM announcement, on December 20, Intel finally acknowledged the situation. They announced a no-questions-asked return policy. They would replace the Pentium chip for anyone that requested it. They set up a number of replacement centers around the country (and around the world, I assume). They set aside $475-million to pay for it all. Their announcement said:

Our previous policy was to talk with users to determine whether their needs required replacement of the processor. To some people, this policy seemed arrogant and uncaring. We apologize. We were motivated by a belief that replacement is simply unnecessary for most people. We still feel that way, but we are changing our policy because we want there to be no doubt that we stand behind this product.

When I first thought about writing a blog post recalling the Pentium affair, I ordered the recent biography of Andy Grove listed in the References. It reads like a ghost-written autobiography. There is a substantial chapter entitled "The Pentium Launch: Intel Meets the Internet" which provides Grove's view of these events and their consequences. He still doesn't quite get it -- he continues to refer to newsgroup traffic as email. But the chapter does reveal that the affair had a deep affect on him. I think his message is well summarized by this quote on a key chain containing a surplus defective Pentium.

Bad companies are destroyed by crises;
Good companies survive them;
Great companies are improved by them.;
--Andy Grove, December 1994

Pentium key chain. Image thanks to Thomas Johansson, thomas@cpucollection.se.

Public Reaction

The entire affair fascinated the public. In addition to the key chains, there were necklaces, pendants, and other jewelry. There were cartoons in newspapers and popular magazines, and jokes in late-night TV monologues. The interest lasted for a couple of months. Some of the jokes, unfortunately some of the worst ones, are preserved in Web sites today. Just Google "Pentium Jokes".

One of my favorite stories is personal. My daughter came to visit for Christmas. She uses a Macintosh and I accidentally dropped it while unloading the car. I took it into a shop in Natick that I had never been in before to see about getting it repaired. My picture had recently been in the local newspaper and a guy behind the counter recognized me. He said, "Hey, you're Cleve Moler. I've got something to show you." He clicked on an icon on a nearby Mac. The icon said "Pentium Powered Calculator". Sure enough, the results of divisions had errors in the fifth decimal place. This is on a Mac! The splash screen said the program came from an outfit called "Sarcastic Software". Who are those guys? I'd love to meet them.

Postmortem

I'm not sure of the actual number, but I understand that relatively few people actually asked to replace their Pentiums after all. Intel did not have to spend very much of the $475-million they had allocated. Most people were satisfied just to be assured that they could get a replacement if they ever felt they really needed it.

Only a few dozen people requested the Pentium-aware version of MATLAB. And none of them reported setting off the floating point division warning message without deliberately trying to trigger it.

So, as far as I know, Thomas Nicely was the only person to ever actually encounter the FDIV bug without actually looking for it. We had other reports of suspected division error sightings, but they were all ultimately traced to arithmetic errors in other places. In fact, the probability of the occurrence of the FDIV error is extremely small. But we didn't know that at the beginning of this affair, so that was why we had to be concerned.

References

A collection of original documents associated with the Pentium division bug is available from MATLAB Central. Pentium Papers

Richard S. Tedlow, Andy Grove, The Life and Times of an American Business Icon, Portfolio, Penguin Books, 2006, 568 pp.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

For example, when I see things like this, it makes my eyes hurt a little:

Here are a few tips to format the Sum block to make your models easier to understand.

Round Sum block

When you drag the Sum block from the Simulink Library Browser, its Icon Shape is set to Round, and it has one input port on the left and one on the bottom:

When the Sum block Icon Shape is set to Round, the ports are spread evenly from top to bottom and the vertical bar "|" can be used to skip one position. For example, if we want to add a port on top of the block, we can modify the list of signs:

Rectangular Sum block

If you change the Icon Shape property of the default Sum block to Rectangular, the block will look like the following:

Personally, I can not think of a good reason for skipping a port in Rectangular shape. So I almost always remove the "|" when I change the shape to Rectangular.

Sum of Elements

One way I like to use the Sum block is to sum all the elements of a vector or matrix:

For matrices, the Sum block can also be configured to Sum only over one specific dimension:

Round or Rectangular?

One question remains: When should the Sum block be Round and when should it be Rectangular?

There is no absolute rule, but personally I like to use the Rectangular shape when implementing equations that flow from left to right, without obvious feedback. For example:

I keep the Round shape for when I want to make it obvious that there is a feedback involved:

Now it's your turn

Please share with us if you have rules or guidelines on the format of your blocks by leaving a comment here.

Brett's Pick this week is "Colorspace Transformations", by Pascal Getreuer.

For inspiration for this week's Pick, I went back to my post from April of last year, in which I asked readers to suggest files to feature. I started to write up Pascal's submission long ago, but somehow got sidetracked, and I wanted to circle back around to it.

File Exchange champ Oliver Woodford casually suggested "a couple of useful ones." (His other selection was featured last August by Jiro.) Both of his "Picks" are quite useful.

So, what's so great about this colorspace converter? In particular, I really like the fact that Pascal's code extends the colorspace conversion capabilities beyond those offered by the Image Processing Toolbox, and that Pascal has included C-code (and instructions for calling it via the MEX interface) to make them faster.

Colorspace Transformations facilitates conversion to any of these colorspaces:

'RGB'              sRGB IEC 61966-2-1
'YCbCr'            Luma + Chroma ("digitized" version of Y'PbPr)
'JPEG-YCbCr'       Luma + Chroma space used in JFIF JPEG
'YDbDr'            SECAM Y'DbDr Luma + Chroma
'YPbPr'            Luma (ITU-R BT.601) + Chroma
'YUV'              NTSC PAL Y'UV Luma + Chroma
'YIQ'              NTSC Y'IQ Luma + Chroma
'HSV' or 'HSB'     Hue Saturation Value/Brightness
'HSL' or 'HLS'     Hue Saturation Luminance
'HSI'              Hue Saturation Intensity
'XYZ'              CIE 1931 XYZ
'Lab'              CIE 1976 L*a*b* (CIELAB)
'Luv'              CIE L*u*v* (CIELUV)
'LCH'              CIE L*C*H* (CIELCH)
'CAT02 LMS'        CIE CAT02 LMS

That's a lot of different ways to represent a color image!

Whenever I am asked to segment a color image--and trust me, that's pretty often!--my initial thought is to look at the R,G, and B colorplanes individually to understand where the information I'm after lies. I do this so often that I wrote and shared a utility that facilitates this exploration process. (In "Advanced Mode," ExploreRGB shows--in addition to RGB [and R, G, and B]--HSV, YCbCr, and L*a*b* versions of the image. I can often find a single-plane version of the original image in one of these transformed colorspaces that makes the segmentation easy. For instance, if I needed to isolate the central yellow pepper in the "peppers.png" image that ships with the Image Processing Toolbox:

I might consider (ironically) starting from the blue chrominance image (i.e., the second colorplane of the YCbCr representation of the image):

Starting with that image plane:

mask = ~im2bw(img,0.24);
cc = bwconncomp(mask);
stats = regionprops(cc,'Area');
A = [stats.Area];
[~,biggest] = max(A);
mask(labelmatrix(cc)~=biggest) = 0;
mask = imfill(mask,'holes');
imshow('peppers.png')
showMaskAsOverlay(0.7,mask,'c');

Working in a non-standard colorspace made this segmentation problem easier than it might otherwise have been, and certainly easier than doing it in RGB-space would be.

(showMaskAsOverlay is a utility I shared for overlaying a transparent mask on an image.)

Okay, so back to Pascal's file.

With a little bit of code, I can use "Colorspace Transforms" to examine quickly many additional representations of the original image:

cs = {'RGB','YCbCr','JPEG-YCbCr','YDbDr','YPbPr','YUV',...
'YIQ','HSV','HSL','HSI',...
'XYZ','Lab','Luv','LCH','CAT02 LMS'};
peppers = imread('peppers.png');
figure('color','w');
ax = tight_subplot(numel(cs),4);
ind = 1;
for ii = 1:numel(cs)
    newim = colorspace(['->' cs{ii}],peppers);
    if max(newim(:)) > 1
        newim = newim/max(newim(:));
    end
    axes(ax(ind));
    imshow(newim,[]);title(cs{ii});
    for jj = 1:3
        ind = ind+1;
        axes(ax(ind));
        imshow(newim(:,:,jj),[]);title(sprintf('%s(:,:,%i)',cs{ii},jj));
    end
    ind = ind+1;
end
expandAxes(ax)

(tight_subplot was a previous Pick-of-the-Week.)

Now I have a lot more options at my disposal for segmenting images!

P.S. I know you can't really see what's happening in those tiny thumbnails*, but then that's the beauty of expandAxes; if you were to run the code snippet above, you would be able to click on any of those tiny axes and expand them to full-frame. And you could right-click on any expanded axes to export the image to the Workspace!

*P.P.S. If your image is large, doing what I just did could be a very bad idea. These aren't really "thumbnails"; even if the visualization is small, each of those image objects has in its "cdata" container a full copy of the image it reflects.

Thank you, Oliver, for the suggestion. And thank you, Pascal, for the very useful submission. Swag on the way to both of you!

As always, I welcome your thoughts and comments. Or leave feedback for Pascal here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Recently someone asked me to explain the speed behavior doing a calculation using a loop and array indexing vs. getting the subarray first.

Contents

Example

Suppose I have a function of two inputs, the first input being the column (of a square array), the second, a scalar, and the output, a vector.

myfun = @(x,z) x'*x+z;

And even though this may be calculated in a fully vectorized manner, let's explore what happens when we work on subarrays from the array input.

I am now creating the input array x and the results output arrays for doing the calculation two ways, with an additional intermediate step in one of the methods.

n = 500;
x = randn(n,n);
result1 = zeros(n,1);
result2 = zeros(n,1);

First Method

Here we see and time the first method. In this one, we create a temporary array for x(:,k) n times through the outer loop.

tic
for k = 1:n
    for z = 1:n
        result1(z) = myfun(x(:,k), z);
    end
    result1 = result1+x(:,k);
end
runtime(1) = toc;

Second Method

In this method, we extract the column of interest first in the outer loop, and reuse that temporary array each time through the inner loop. Again we see and time the results.

tic
for k = 1:n
    xt = x(:,k);
    for z = 1:n
        result2(z) = myfun(xt, z);
    end
    result2 = result2+x(:,k);
end
runtime(2) = toc;

Same Results?

First, let's make sure we get the same answer both ways. You can see that we do.

theSame = isequal(result1,result2)

theSame =
     1

Compare Runtime

Next, let's compare the times. I want to remind you that doing timing from a script generally has more overhead than when the same code is run inside a function. We just want to see the relative behavior so we should get some insight from this exercise.

disp(['Run times are: ',num2str(runtime)])

Run times are: 2.3936      1.9558

What's Happening?

Here's what's going on. In the first method, we create a temporary variable n times through the outer loop, even though that array is a constant for a fixed column. In the second method, we extract the relevant column once, and reuse it n times through the inner loop.

Be thoughtful if you do play around with this. Depending on the details of your function, if the calculations you do each time are large compared to the time to extract a column vector, you may not see much difference between the two methods. However, if the calculations are sufficiently short in duration, then the repeated creation of the temporary variable could add a tremendous amount of overhead to the calculation. In general, you should not be worse off always capturing the temporary array the fewest number of times possible.

Your Results?

Have you noticed similar timing "puzzles" when analyzing one of your algorithms? I'd love to hear more here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Simulation Accuracy

I recently had to validate system behavior and analyze numerical accuracy of a model. Making small modifications in one part of the model was leading to unexpectedly large differences at the outputs. Since the model was very large, I thought it would be painful to insert Scopes or To Workspace blocks in the model to try identifying the source of the problem.

I found out that the Simulation Data Inspector can be very useful to compare simulation results without adding blocks to your model and without writing any scripts. Let's see how this work!

Logging Signals

The first step is to configure signals for logging. I like to use the Model Explorer for that. It is significantly faster than right-clicking on the signal lines one by one. Here are a few tips to configure the Model Explorer with the goal of controlling signal logging.

Also, ensure that signal logging is enabled in the model configuration.

Record & Inspect Simulation Output

In your model, click on the Record button.

When the simulation is completed, open the Simulation Data Inspector using the link that pops up.

Compare Runs

Make some changes to your model and simulate it again. In the Simulation Data Inspector, go to the Compare Runs Tab, select your 2 runs and click Compare.

Using the absolute tolerance I specified, the Simulation data Inspector helped me to quickly identify where the divergence appear in my model.

Now it's your turn

Give a look and the Validate System Behavior section of the documentation to learn more about the Simulation Data Inspector and let us know if it helps your workflow by leaving a comment here. ]]> http://feedproxy.google.com/~r/mathworks/pick/~3/vtFfv1dlqTg/ Doug's pick this week is Raspberry Pi DC Motor H-Bridge Driver Block by Joshua Hurst. As my previous posts may suggest, I've been a long-time fan of the Arduino platform. Among other things I use it as the basis of a demonstration to showcase a variety of MathWorks' tools that includes driving a motor to [...] http://blogs.mathworks.com/pick/?p=4536 Fri, 03 May 2013 13:00:27 +0000 Doug's pick this week is Raspberry Pi DC Motor H-Bridge Driver Block by Joshua Hurst.

As my previous posts may suggest, I've been a long-time fan of the Arduino platform. Among other things I use it as the basis of a demonstration to showcase a variety of MathWorks' tools that includes driving a motor to rotate a camera to track an object. The Arduino works great for doing closed-loop position control of the camera, but falls short when it comes to actually processing the camera video to detect the object. So currently the video processing is relegated to my laptop which prevents me from creating a completely embedded solution.

Enter R2013a which introduced support for a number of new target platforms which Jiro pointed out earlier this year. I was drawn to Raspberry Pi in particular because it is a similar price and form factor to the Arduino, but has built-in support for USB webcams:



Looking at the blocks that are included with the target, however, you'll notice that there is nothing included for sending voltages to a DC Motor like you can do with the Arduino PWM Output blocks. Because of that I had initially planned to use both the Raspberry Pi and Arduino for my demo until I ran across Joshua's submission the other day. His model uses S-Function Builder blocks to take advantage of the new WiringPi support in the Raspberry Pi image that is used in the support package. Here's the top-level of Joshua's model:



After some breadboard work, I built and downloaded the model and I was controlling my motor in minutes by simply changing the Input Volts value in External Mode. Now I just need to decide on the best way to read the motor's potentiometer then the whole system can run on the Pi. Or perhaps I should bite the bullet and switch to an encoder. If anyone has any suggestions, leave a comment!

I was happy to see that Joshua was inspired by Giampiero's Device Driver guide that explains how to create custom drivers for the hardware support packages. Joshua has also recently uploaded several other related submissions. I'm looking forward to seeing what Joshua and others come up with in the future to increase the I/O options for the different hardware platforms.

One minor suggestion I would have for Joshua is to include the pin numbers that correspond to the PWM outputs for this model. I had to do a bit of hunting to track that down (it's 7 and 11 on the P1 header for those that are interested). But other than that, this submission has been a big help, thanks Joshua!

Comments
Let us know what you think here or leave a comment for Joshua.]]> Simulink http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/2mFpr8eP0sw/ Sometimes when you are working on a User Interface you want to use the profiler to speed up the code. However, with a UI so much time can be spent by the user rather than the code of interest that the timing is not as precise as would be preferred. I like this technique so [...] http://blogs.mathworks.com/videos/?p=1101 Fri, 03 May 2013 12:51:07 +0000

I like this technique so that I can do ad hock timings as needed in a complicated User Interface in MATLAB.

I was reviewing enhancement requests recently when I came across this one: Add support for a 'Quality' parameter when using imwrite to write a PNG file, just like you can currently do when writing a JPEG file.

Well, there is actually a pretty good reason why there is no 'Quality' parameter for writing PNG files, but it's not an obvious reason unless you know more about the differences between various image compression methods.

A very important characteristic of a compression method is whether it is lossless or lossy. With a lossless compression method, the original data can be recovered exactly. When you make a "zip" file on your computer, this is what you certainly expect.

Here's an example: I can use gzip to compress the MATLAB script file I'm using for this blog post.

gzip('lossless_lossy.m');
dir('lossless_lossy.*')

lossless_lossy.m     lossless_lossy.m.gz  

Let's compare their sizes to see if the output of gzip is really compressed.

d1 = dir('lossless_lossy.m');
d1.bytes

ans =

        4098

d2 = dir('lossless_lossy.m.gz');
d2.bytes

ans =

        1785

Indeed, the compressed does actually have fewer bytes. Can we get the original file back exactly?

gunzip('lossless_lossy.m.gz','./tmp')
isequal(fileread('lossless_lossy.m'),fileread('./tmp/lossless_lossy.m'))

ans =

     1

Yes.

Let's switch to image formats. The PNG image format uses lossless compression. When you save image data to a PNG file, you can read the file back in and get back the original pixels, unchanged. For a sample image I'll use my imzoneplate function on the MATLAB Central File Exchange.

I = im2uint8(imzoneplate);
imshow(I)

Let's write I out to a PNG file, read it back in, and see if the pixels are the same.

imwrite(I,'zoneplate.png');
I2 = imread('zoneplate.png');
isequal(I,I2)

ans =

     1

OK, now let's try the same experiment using JPEG.

imwrite(I,'zoneplate.jpg')
I2j = imread('zoneplate.jpg');
isequal(I,I2j)

ans =

     0

No, the pixels are not equal! It turns out the JPEG is a lossy image compression format. (Full disclosure: there is a lossless variant of JPEG, but it is rarely used.)

Why in the world would we use a compression format that doesn't preserve the original data? Because by giving up on exact data recovery and by taking advantage of properties of human visual perception, we can make the stored file a lot smaller. Let's compare the file sizes of the PNG file with the JPEG file.

z1 = dir('zoneplate.png');
num_bytes_png = z1.bytes

num_bytes_png =

      218864

z2 = dir('zoneplate.jpg');
num_bytes_jpeg = z2.bytes

num_bytes_jpeg =

       72660

The JPEG file is only one-third the size of the PNG file! But it looks almost exactly the same.

imshow('zoneplate.jpg')

So, what about that 'Quality' parameter that I mentioned at the top of today's post? It turns out that we can make the JPEG file even smaller if we are willing to put up with some visible distortion in the image. Let's try a quality factor of 25. (The default is 75 on a scale of 0 to 100.)

imwrite(I,'zoneplate_25.jpg','Quality',25)
I2j_25 = imread('zoneplate_25.jpg');
imshow(I2j_25)

You can see some distortion, especially around the high-frequency part of the pattern. How about the file size?

z3 = dir('zoneplate_25.jpg');
z3.bytes

ans =

       39544

That's about 54% of the size of the zoneplate.jpg.

Sometimes people think they can get "lossless JPEG" by using a 'Quality' factor of 100, but unfortunately that isn't the case. Let's check it:

imwrite(I,'zoneplate_100.jpg','Quality',100)
I2j_100 = imread('zoneplate_100.jpg');
isequal(I,I2j_100)

ans =

     0

So there you go -- that's why there's no 'Quality' parameter when writing PNG files. PNG files are always perfect!

When I write blog posts, sometimes I use PNG image files and sometimes I use JPEG. My choice is based on what kind of graphics I have in that particular post. And that's a blog topic for another day.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

The data cursor feature has since then evolved from tap-and-hold invocation to a more prominent ‘Show Data Cursor’ button. Tapping on this button brings up the crosshair. Move the crosshair to any location on the figure to display the coordinates of that location.

We’re happy to announce that the latest release of MATLAB Mobile for Android supports this feature, much like its iOS counterpart.

                     Tap on Show Data Cursor                                    Position crosshair at desired location

Download MATLAB Mobile from Google Play.

If you have used the data cursor functionality on your mobile device, we would love to hear from you. Leave us a comment here with your thoughts.

The Pentium division bug episode in the fall of 1994 was a defining moment for the MathWorks, for the Internet, for Intel Corporation, and for me personally. In this blog I am reprinting the article that I wrote for the Winter 1995 issue of MATLAB News and Notes and for SIAM News. In my next blog I want to discuss the episode's impact.

Contents

A Tale of Two Numbers (Reprinted from MATLAB News and Notes, Winter 1995)

This is the tale of two numbers, and how they found their way over the Internet to the world's front pages on Thanksgiving day, embarrassing the world's premier computer chip manufacturer.

Thomas Nicely

The first number is a twelve digit integer

$$ p = 824633702441 $$

Thomas Nicely of Lynchburg College in Virginia is interested in this number because both p and p+2 are prime. Two consecutive odd numbers that are both prime are called twin primes. Nicely's research involves the distribution of twin primes and, in particular, the sum of their reciprocals,

$$ S = 1/5 + 1/7 + ... + 1/29 + 1/31 + ... + 1/p + 1/(p+2) + ... $$

It is known that this infinite sum has a finite value, but nobody knows what that value is. For over a year, Nicely has been carrying out various computations involving twin, triple and quadruple primes. One of his objectives has been to demonstrate that today's PCs are just as effective as supercomputers for this kind of research. He was using half a dozen different computers in his work until March, when he added a machine based on Intel Corporation's Pentium chip to this collection.

In June, Nicely noticed inconsistencies among several different quantities he was computing. He traced it to the values he was getting for $1/p$ and $1/(p+2)$. Although he was using double precision, they were accurate only to roughly single precision. He blamed at first on his own programs, then the compiler and operating system, then the bus and motherboard in his machine. By late October though, he was convinced that the difficulty was in the chip itself.

Alex Wolfe

On October 30, in an e-mail message to several other people who had access to Pentium systems, Nicely said that he thought there was a bug in the processor's floating point arithmetic unit. One of the recipients of Nicely's memo posted it on the CompuServe network. Alex Wolfe, a reporter for the EE Times, spotted the posting and forwarded it to Terje Mathisen.

Terje Mathisen

Terje Mathisen, a computer jock at Norsk Hydro in Oslo, Norway, has published articles on programming the Pentium, and posted notes on the Internet about the accuracy of its transcendental functions. Within hours of receiving Wolfe's query, Mathisen confirmed Nicely's example, wrote a short, assembly language test program, and, on November 3, initiated a chain of Internet postings in newsgroup comp.sys.intel about the FDIV bug. (FDIV is the floating point divide instruction on the Pentium.) A day later, Germany's Andreas Kaiser posted a list of two dozen numbers whose reciprocals are computed to only single precision accuracy on the Pentium.

Tim Coe

Tim Coe is an engineer at Vitesse Semiconductor in southern California. He has designed floating point arithmetic units himself and saw in Kaiser's list of erroneous reciprocals clues to how other designers had tackled the same task. He surmised, correctly it turns out, that the Pentium's division instruction employs a technique known to chip designers as a radix 4 SLT algorithm. This produces two bits of quotient per machine cycle and thereby makes the Pentium twice as fast at division as previous Intel chips running at the same clock rate.

Coe created a model that explained all the errors Kaiser had reported in reciprocals. He also realized that division operations involving numerators other than one potentially had even larger relative errors. His model lead him to a pair of seven digit integers whose ratio

$$ c = 4195835 / 3145727 $$

appeared to be an instance of the worst case error. Coe did his analysis without actually using a Pentium -- he doesn't own one. To verify his prediction, he had to bundle his year-and-a-half old daughter into his car, drive down to a local computer store, and borrow a demonstration machine.

This true value of Coe's ratio is

$$ c = 1.33382044... $$

But computed on a Pentium, it is

$$ c = 1.33373906... $$

The computed quotient is accurate to only 14 bits. The relative error is $6.1 \cdot 10^{-5}$. This is worse than single precision roundoff error and over ten orders of magnitude worse than what we expect from double precision computation. The error doesn't occur very often, but when it does, it can be very significant. We are faced with a very small probability of making a very big error.

Press Coverage

I first learned of the FDIV bug a few days before Coe's revelation from an electronic mailing list maintained by David Hough; at that point I began to follow comp.sys.intel. At first, I was curious but not alarmed. It seemed to be some kind of single- versus double-precision problem which, although annoying, is not all that uncommon. But Coe's discovery, together with a couple of customer calls to The MathWorks tech support, raised my level of interest considerably.

On November 15, I posted a summary of what I knew then to both the Intel and the MATLAB newsgroups, using Nicely's prime and Coe's ratio as examples. I also pointed out that the divisor in both cases is a little less than 3 times a power of two:

$$ 824633702441 = 3 \cdot 2^{38} - 18391 $$

and

$$ 3145727 = 3 \cdot 2^{20} - 1 $$

By this time, the Net had become hyperactive and my posting was redistributed widely. A week later, reporters for major newspapers and news services had Xeroxed copies of faxed copies of printouts of my posting.

The story began to get popular press coverage when Net activists called their local newspapers and TV stations. I was interviewed by CNN on Tuesday, the 22nd, and spent most of the next day on the telephone with other reporters. My picture was in the New York Times and the San Jose Mercury News on Thursday which happened to be Thanksgiving. The Times story included a sidebar, titled "Close, But Not Close Enough", which used Coe's ratio to illustrate the problem. The story was front-page news in the Boston Globe, with the headline "Sorry, Wrong Number." The Globe's sidebar demonstrated the error by evaluating Coe's ratio in a spreadsheet.

In the month since I learned of Nicely's prime and Coe's ratio they have certainly earned their place in the Great Numbers Hall of Fame.

One challenging aspect of this has been how to explaing to the press how big the error is. The focus of attention has been on what we numerical analysts call the residual,

$$ r = 4195835 - (4195835 / 3145727) (3145727) $$

With exact computation, r would be zero, but on the Pentium it is 256. To most people, 256 seems like a pretty big error. But when compared to the input data, of course, it doesn't seem so large. This gets us into relative error and significant digits -- terms that are not encountered in everyday journalistic prose. There was even confusion on the part of some reporters about where to starting counting "decimal digits". Not everybody got the details right. The New York Times was the only publication I saw that thought to show the actual value of the erroneous quotient. The British publication, The Economist, had the good sense to describe the error as 0.006%. In retrospect, perhaps the best description of the error is "about 61 parts per million."

Since double precision floating point computation is vital to MATLAB, and since Pentium-based machines account for a significant portion of new MATLAB usage, we decided early on to provide a release of our software that works around the bug. My first thought was to modify our source code so that every division operation

   |z = x/y|

was followed by a residual calculation

   |r = x - y*z|

When the residual is too large, the division could be redone in software using Newton's method. We abandoned both aspects of this approach -- we don't compute the residual, and we don't do software emulation.

Peter Tang

Coe, Mathisen and I are now working with Peter Tang of Argonne Laboratory (on sabbatical in Hong Kong) and a group of hardware and software engineers from Intel, who are in California and Oregon. We have never met face-to-face. Our collaboration is all by e-mail and telephone. (11 am in Massachusetts is 8 am in California and Oregon, 5 pm in Oslo and 1 am in Hong Kong.) We are developing, implementing, testing and proving the correctness of software that will work around the FDIV bug and related bugs in the Pentium's on-chip tangent, arctangent and remainder instructions. We will offer the software to compiler vendors, commercial software developers and individuals folks via the Net. We hope this software will be used by anyone doing floating point arithmetic on a Pentium who is unable to replace the chip.

SRT Algorithm

The Pentium's division woes can be traced to five missing entries in a lookup table that is actually part of the chip's circuitry. The radix-4 SRT algorithm is a variation of the familiar long division technique we all earned in school. Each step of the algorithm takes one machine cycle and appends another two-bit digit to the quotient. Both the quotient and the remainder are represented in a redundant fashion as the difference between two numbers. In this way, the next digit of the quotient can be obtained by table lookup, with several high-order bits from both the divisor and remainder used as indices into the table. The table is not rectangular; a triangular portion of it is eliminated by constraints on the possible indices.

When the algorithm was implemented in silicon, five entries on the boundary of this triangular portion were dropped. These missing entries (each of which should have a value of +2) effectively mean that, for certain combinations of bits developed during the division process, the chips make an error while updating the remainder.

Our Workaround

The key to our workaround is the fact the chip does a perfectly good job with division almost all the time. It is simply a question of avoiding the operands that don't work very well, which our software does with a quick test of the divisor before each attempted FDIV. The absence of dangerous bit patterns indicates that the FDIV can be done safely. The presence of one of the patterns does not guarantee that the error will occur, just a signal that it might. In this case, scaling both numerator and denominator by 15/16 takes the divisor out of the unsafe region and insures that the subsequent division will be fully accurate.

With this approach, it is not necessary to test the magnitude of the residual resulting from a division. It is known a priori that all divisions will produce fully accurate results. If desired, an additional test can compare the result of scaled and unscaled divisions and count the number of errors that would occur on an unmodified Pentium. We will offer this test in MATLAB, but it can be turned off for maximum execution speed.

The regions to be avoided can be characterized by examining the high order bits of the fraction in the IEEE floating point representation. Floating point numbers in any of the three available precisions can be thought of as real numbers of the form

$$ \pm (1 + f ) \cdot 2^e $$

where e is an integer in the range determined by under- and overflow and f is a rational number in the half open interval [0, 1). The eight high order bits of f divide this interval into $2^8$ = 256 subintervals. Five of these subintervals contain all the at-risk divisors. The hexadecimal representations of the leading bits of the five subintervals are

  1F, 4F, 7F, AF, and DF

(As I write this article, Coe, Tang and some of the Intel people are still working on a proof of correctness. The ultimate version of the algorithm may well use more than eight bits, but the idea will be the same.)

A picture of the floating point numbers between any two powers-of-two looks something like this

  [____()______()______()______()______()____]
           1F      4F      7F      AF      DF

The length of each of the five subintervals is 1/256 of the overall length. The result of any division involving a divisor outside these subintervals will be accurate. Only a small fraction of divisions involving divisors inside the subintervals produce erroneous results, but it is easiest, and quickest, to avoid the subintervals altogether. Scaling the numerator and denominator by 15/16 shifts it to a safe region. (I had originally proposed scaling by 3/4, but this takes the 7Fthird subinterval into the 1F subinterval.)

For example, the denominator in Coe's ratio is

$$ 3145727 = 1.49999952316284 \cdot 2^{21} $$

The hexadecimal representation is printed as 4147FFFF80000000. So e= 21 and f = $2^{-1} - 2^{-21}$, and the number is close to the right end point of the 7F subinterval. The 17 consecutive high order 1 bits in f conspire with the bits in Coe's numerator to produce an instance of worst-case error.

The Pentium, like most other microprocessors, saves floating point numbers in memory with either a 32-bit single precision format, or a 64-bit double precision format. In the floating point arithmetic unit itself, however, an 80-bit extended-precision format is used for the calculations. Thus, if an at-risk denominator is encountered in either single or double precision, it can be scaled by 15/16 and extended precision can be used for the division. Then, when the resulting quotient is rounded back to the working precision, it will yield exactly the same result, down to the last bit, as would be produced by a chip without the bug.

Pentium-Safe Division

The entire algorithm can be summarized in a few lines of MATLAB pseudo-code. This doesn't actually work, because MATLAB doesn't do bit-level operations on floating point numbers. Moreover, we do not make the distinction between the 64- and 80-bit representations that ensure full accuracy when scaling is required. But I hope this explanation has given readers an idea of the thinking that went into the workaround.

The algorithm is:

  function z = fdiv(x,y)
    if at_risk(y)
       x = (15/16)*x;
       y = (15/16)*y;
    end
    z = x/y;
  end

The safety filter is:

  function a = at_risk(y)
     f = and(hex(y),'000FF00000000000');
     a = any(f == ['1F'; '4F'; '7F'; 'AF'; 'DF'])
  end

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Last week I received the following question and example model:

In the attached model, I implemented a basic bang-bang controller for a plow. Using a Sign block, I generate a 1 for the plow to go up, -1 to go down and 0 to stop. The problem is that the plow never seems to stabilize around zero. Can you help me understand why?

Here is an image of the model:

Nonadaptive zero-crossing detection

With the default settings, this model errors out as soon as we try to give a non-zero command with the following error:

At time 1.0000000000236877, simulation hits (1000) consecutive zero crossings. Consecutive zero crossings will slow down the simulation or cause the simulation to hang. To continue the simulation, you may 1) Try using Adaptive zero-crossing detection algorithm or 2) Disable the zero crossing of the blocks shown in the following table.
--------------------------------------------------------------------------------
Number of consecutive zero-crossings : 1000
Zero-crossing signal name : Input
Block type : Signum
Block path : 'myPlowController/Sign'
--------------------------------------------------------------------------------
You can turn off this message by using the MATLAB command:
set_param('myPlowController','MaxConsecutiveZCsMsg','none');

This error is relatively simple to understand by looking at the ideal plant model (an Integrator block) and the switching logic. If the plow position is slightly negative, the controller makes it move upward, leading to a positive error. The controller sees this positive error, makes the plow go down, bringing us back to our original state .

With zero-crossing detection enabled for the Sign block, Simulink will detect that we crossed zero and try to reduce its step-size to capture this transition accurately. But here, reducing the step size will not help. The step size can be a small as possible, a small negative error will always lead to a jump on the positive side, and vice-versa.

Let's try the suggestions offered by the error message. What suggestions? If you didn't read the error message, go back and read it now.

Suggestion 1: Adaptive zero-crossing detection

In the solver section of the model configuration, we can try selecting adaptive zero-crossing detection. This will dynamically activate and deactivate zero-crossing bracketing, helping with models that exhibit strong chattering.

With this setting the simulation completes, but the variable step solver takes way too many steps when it is not needed, and the position error is too large.

The reason is that the adaptive algorithm encountered multiple consecutive zero-crossings for the Sign block and decided to disable zero-crossing detection for this block. Once this is done, the solver takes larger steps and gives inaccurate results.

Let's try suggestion 2 from the error message.

Suggestion 2: Disabling zero-crossing

The second suggestion from the error message is to disable zero-crossing detection for the problematic block.

In this case, the tracking is significantly better because the solver uses the state of the integrator to limit its step size, but we see a lot of switching happening when the plow desired height is constant.

The main problem here is that it is impossible for the Simulink engine to stop exactly on the zero value of the Sign block with this system. With the equations involved, Simulink can take very small steps around zero, but will not stop exactly on it. This not what the Sign block is designed to do.

Possible solutions

If we want a system that can be simulated efficiently using a variable step solver, we have to modify the equations. If we want to keep the plant as it is, we need a controller with a logic that will stop when the error is within a certain range around zero. For example, we could use a Stateflow chart or a block with hysteresis like the Relay block:

In this case, the plow stops moving when the command is fixed and the zero-crossing detection helps ensure that the solver takes appropriate steps when needed.

Now it's your turn

If you have some interesting plowing or zero-crossing stories to share, leave a comment here

Jiro's pick this week is togglefig by our own Brett Shoelson.

Brett and I, along with Bob, started writing for this blog about 5 years ago, back when Doug was still the owner of the blog. Since the beginning, we typically highlighted submissions that were written by the general public and not by one of us. But ultimately, we like to highlight files that can be of some value to our readers, and I must say that there's no shortage of such files in Brett's entries.

So what does togglefig do? Have you ever wanted to better manage your figure windows? I've often had a need to reuse my figure windows. Once, I was developing some piece of code that created multiple figures, and as I was iterating on my code, I kept creating the same number of windows each time, and I would end up with tens of figure windows open after a while. Of course I could "close all" of them at the beginning of my script, but I had some windows I needed open for other reasons. As a solution, I would create/select figures using some numbers, e.g.

figure(123)

This would always use the same figure window, or create one if it doesn't exist. But I would want to give a meaningful name to it, so then I would need to set the name property after getting the handle to the figure.

h = figure(123);
set(h, 'Name', 'Important Window');

It works, but it's just some extra lines of code that are not adding any value to the real purpose of my program.

Brett's togglefig addresses this in a very elegant manner. It would raise the figure with the specified name, or create one if it doesn't exist.

togglefig('More Important Window')

You can also pass in an optional argument to clear the figure window when it's given focus.

Thanks, Brett, for this extremely useful function that makes me more efficient in using MATLAB!

Comments

Let us know what you think here or leave a comment for Brett.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Once more I am pleased to introduce guest blogger Kai Gehrs. Kai has been a Software Engineer at MathWorks for the past five years working on the Symbolic Math Toolbox. He has a background in mathematics and computer science. He already contributed to my BLOG in the past writing about Using Symbolic Equations And Symbolic Functions In MATLAB as well as on approaches for Simplifying Symbolic Results.

Contents

In a Nutshell: What Is This Article About?

If you are interested in using MATLAB and the Symbolic Math Toolbox in teaching some basics in mechanical engineering, this might be of interest to you.

Computing area moments of inertia is an important task in mechanics. For example, area moments of inertia play a critical role in stress, dynamic, and stability analysis of structures. In this article, we use capabilities of the Symbolic Math Toolbox to compute area moments for cross sections of elliptical tubes.

We start with a basic case involving only numeric parameters, and then make the computations more general by introducing symbolic parameters.

All plots used in this article have been created in the MuPAD Notebook app. You can find the source code for these plots at the end of the article.

Basic Example: Cross Section of an Elliptical Tube

The following picture shows the cross section of an elliptical tube.

The following ellipses define outer and inner contours of the section:

We will compute the area moment of inertia of this section with respect to the $x$-axis.

Area Moment Of Inertia

The moment of inertia of an area $A$ with respect to the $x$-axis is defined in terms of the double integral

$$I_x = {\int \int}_A y^2\, \mathrm{d}A.$$

You can find this definition on Wikipedia.

Math Behind This Example

In our example, the area $A$ is the hatched area shown in the previous plot. Taking advantage of the symmetry of the section with respect to both the $x$- and $y$-axes, we can restrict our computations to the first quadrant of the $x-y$-plane.

To compute the necessary double integral over the hatched area, we divide this area into two separate areas $A1$ and $A2$.

To compute the final area moment of inertia about the x-axis, we multiply the sum of the double integrals over $A_1$ and $A_2$ by four.

$$I_x = {\int\int}_A y^2 \, \mathrm dA = 4 \cdot \Big({\int\int}_{A_1} y^2 \, \mathrm dA_1 + {\int\int}_{A_2} y^2 \, \mathrm dA_2\Big).$$

Now, we only need to compute these two double integrals. The mathematical theory behind multi-dimensional integration tells us that we can rewrite each of these double integrals in terms of two integrals:

$$I_1 = {\int\int}_{A_1} y^2 \, \mathrm dA_1 = \int_0^2 \int_{\sqrt{1 - \frac{x^2}{4}}}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2 \, \mathrm dy \, \mathrm dx,$$

$$I_2 = {\int\int}_{A_2} y^2 \, \mathrm dA_2 = \int_2^3 \int_{0}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2\, \mathrm dy\, \mathrm dx.$$

The int function from the Symbolic Math Toolbox can compute these integrals.

Symbolic Integration

We define $x$ and $y$ as symbolic variables:

syms x y;

We start computing the inner integral $\int_{\sqrt{1 - \frac{x^2}{4}}}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2 \, \mathrm dy$ of $I_1$:

innerI1 = int(y^2, y, sqrt(1-x^2/4), 2*sqrt(1-x^2/9))

innerI1 =
(8*(9 - x^2)^(3/2))/81 - (4 - x^2)^(3/2)/24

Next, we integrate innerI1 with respect to $x$ from $0$ to $2$. This provides $I_1$:

I1 = int(innerI1, x, 0, 2)

I1 =
3*asin(2/3) - pi/8 + (74*5^(1/2))/81

We use the same strategy to compute $I_2$. First, we compute the inner integral $\int_{0}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2\, \mathrm dy$. Then, we integrate the resulting expression with respect to $x$ from $2$ to $3$:

I2 = int(int(y^2, y, 0, 2*sqrt(1-x^2/9)), x, 2, 3);

Hence, the area moment of inertia of the elliptical tube with respect to the $x$-axis is

Ix = 4 * (I1 + I2)

Ix =
(11*pi)/2

We can approximate the result numerically by using the vpa function. For example, we approximate the result using 5 significant (nonzero) digits:

vpa(Ix, 5)

ans =
17.279

Advanced Example: Cross Section of an Elliptical Tube Defined by Symbolic Parameters

We can use a numerical integrator, such as MATLAB's integral2, to compute the area moment of inertia in the previous example. A numerical integrator might return slightly less accurate results, but other than that there is not much benefit from using symbolic integration there.

But what if we consider the cross section of a more general elliptical tube whose shape is defined by arbitrary symbolic parameters?

For example, we want to compute the area moment of inertia of this elliptical tube.

Its contour lines are still ellipses, but they are defined by the symbolic parameters $r_1$, $r_2$, $R_1$ and $R_2$:

Now we need to compute these integrals:

$$I_1 = {\int\int}_{A_1} y^2 \, \mathrm dA_1 = \int_0^{r_1} \int_{r_2 \, \sqrt{1 - \frac{x^2}{r_1^2}}}^{R_2 \, \sqrt{1 - \frac{x^2}{R_1^2}}} y^2 \, \mathrm dy \, \mathrm dx,$$

$$I_2 = {\int\int}_{A_2} y^2 \, \mathrm dA_2 = \int_{r_1}^{R_1} \int_{0}^{R_2 \, \sqrt{1 - \frac{x^2}{R_1^2}}} y^2\, \mathrm dy\, \mathrm dx.$$

Although the situation is more complicated now, we can still use the same strategy as above, using symbolic variables instead of numbers. Nevertheless, we need to add one more step: specify relationships between symbolic parameters. This is because the variables $r_1$, $r_2$, $R_1$, and $R_2$ can be any complex number, unless we explicitly restrict their values. For example, the system does not know if these variables are positive or negative, if $r_1 < R_1$ or otherwise. In this example, we want to specify that $r_1>0$, $r_2>0$, $R_1>r_1$, and $R_2>r_2$. In the Symbolic Math Toolbox, such restrictions are called assumptions on variables. They can be set by using the assume and assumeAlso functions:

syms x y r1 r2 R1 R2;
assume(r1 > 0);
assume(r2 > 0);
assumeAlso(r1 < R1);
assumeAlso(r2 < R2);
I1 = int(int(y^2, y, r2*sqrt(1-x^2/r1^2), R2*sqrt(1-x^2/R1^2)), x, 0, r1);
I2 = int(int(y^2, y, 0, R2*sqrt(1-x^2/R1^2)), x, r1, R1);
Ixgeneral = 4 * (I1 + I2);
pretty(Ixgeneral)

 
        /          4     / r1 \ \         /                     4     / r1 \ \ 
        |      3 R1  asin| -- | |         |             4   3 R1  asin| -- | | 
      3 |                \ R1 / |       3 |      3 pi R1              \ R1 / | 
  4 R2  | #1 + ---------------- |   4 R2  | #1 - -------- + ---------------- | 
        \             8         /         \         16             8         / 
  ------------------------------- - ------------------------------------------ 
                   3                                      3 
               3 R1                                   3 R1 
   
                3 
        pi r1 r2 
      - --------- 
            4 
   
  where 
   
           /     2        3 \ 
           | 5 R1  r1   r1  |    2     2 1/2 
     #1 == | -------- - --- | (R1  - r1 ) 
           \    8        4  /

Now, let's substitute our symbolic parameters with the same values that we had in the first example. Do we get the same result? Let's double-check: when substituting $r_1=2,r_2,=1,R_1=3,R_2=2$ in Ixgeneral, do we get the same result for Ix as obtained for the section initially considered?

As expected, the general solution Ixgeneral reduces to the specific solution Ix when we substitute our original dimensions:

vpa(subs(Ixgeneral, [r1, R1, r2, R2],[2, 3, 1, 2]), 5)

ans =
17.279

Note that it is essential to set the assumptions on $r_1$, $r_2$, $R_1$ and $R_2$. By default, the Symbolic Math Toolbox assumes that all variables including symbolic parameters are arbitrary complex numbers. This makes computing integrals much more complicated and time-consuming. For example, try computing Ixgeneral without setting the assumptions on $r_1$, $r_2$, $R_1$ and $R_2$. In general, it can be impossible to even get the result without any additional assumptions.

How to Create MuPAD Graphics Shown in This Article

To generate the graphics shown in this article, use MuPAD Notebook app. Note that the code in this section does not run in the MATLAB Command Window.

To open the MuPAD Notebook app:

Alternatively, type mupad in the MATLAB Command Window.

Paste the following commands to a MuPAD notebook to obtain the four graphics used above:

Ellipse:= (x,a,b) -> sqrt(b^2*(1-x^2/a^2)):

f1 := plot::Function2d(Ellipse(x,2,1), x = -3..3):
f2 := plot::Function2d(Ellipse(x,3,2), x = -3..3):
f3 := plot::Function2d(Ellipse(x,3,2), x = -3..-2):
f4 := plot::Function2d(Ellipse(x,3,2), x = 2..3):

f5 := plot::Function2d(-Ellipse(x,2,1), x = -3..3):
f6 := plot::Function2d(-Ellipse(x,3,2), x = -3..3):
f7 := plot::Function2d(-Ellipse(x,3,2), x = -3..-2):
f8 := plot::Function2d(-Ellipse(x,3,2), x = 2..3):

f9 := plot::Function2d(Ellipse(x,2,1), x = 0..2, VerticalAsymptotesVisible = FALSE):
f10 := plot::Function2d(Ellipse(x,3,2), x = 0..3, VerticalAsymptotesVisible = FALSE):
f11 := plot::Function2d(Ellipse(x,3,2), x = 2..3, VerticalAsymptotesVisible = FALSE):

plot(plot::Hatch(f1,f2),plot::Hatch(f3),plot::Hatch(f4),f1,f2,
     plot::Hatch(f5,f6),plot::Hatch(f7),plot::Hatch(f8),f5,f6,
     Scaling = Constrained,
     GridVisible = TRUE,
     VerticalAsymptotesVisible = FALSE,
     Height = 120, Width = 200,
     Header = "Cross section of elliptical tube (x-y-plane)"):

plot(plot::Hatch(f9,f10),plot::Hatch(f11),f9,f10,
     Scaling = Constrained,
     GridVisible = TRUE,
     Height = 120, Width = 200,
     Header = "Restriction to first quadrant (x-y-plane)"):

plot(plot::Hatch(f9,f10,Color=RGB::Magenta),plot::Hatch(f11,Color = RGB::Green),f9,f10,
     Scaling=Constrained,
     GridVisible = TRUE,
     plot::Text2d("A1",[1.0,1.25]),
     plot::Text2d("A2",[2.5,0.25]),
     Height = 120, Width = 200,
     Header = "Integration areas"):

plot(plot::Hatch(f1,f2),plot::Hatch(f3),plot::Hatch(f4),f1,f2,
     plot::Hatch(f5,f6),plot::Hatch(f7),plot::Hatch(f8),f5,f6,
     plot::Text2d("r1",[1.8,0.05]),
     plot::Text2d("r2",[0.05,0.85]),
     plot::Text2d("R1",[2.74,0.05]),
     plot::Text2d("R2",[0.05,1.84]),
     Scaling = Constrained,
     GridVisible = TRUE,
     VerticalAsymptotesVisible = FALSE,
     XTicksLabelsVisible = FALSE,
     YTicksLabelsVisible = FALSE,
     ViewingBox = [-3.5..3.5,-2.5..2.5],
     Height = 120, Width = 200,
     Header = "Cross section of elliptical tube (more general situation)"):

Have You Tried Symbolic Math Toolbox?

Have you used the Symbolic Math Toolbox in computations related to mechanics? Let me know here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

In my previous post, we saw how to estimate continuous transfer functions with System Identification Toolbox. We estimated the following transfer function for a simple DC Motor using tfest:

For this transfer function, we designed the following controller using pidtune:

We will now implement the controller on the Arduino Uno and see how the DC motor fares with this controller. To deploy the controller on the hardware, we will use Simulink’s capability to generate an executable and run it on selected hardware.

Deploying controller to the Arduino board

You probably noticed that the controller shown above is in a continuous form. To use it on our target, the first thing to do is to discretize it using the c2d function:

Then we grab the PID block from the Simulink Library and configure it.

To keep the PID controller’s output within the limits of the hardware, we go to the PID Advanced tab and enable output saturation along with anti-windup protection.

To test the controller on the hardware, we created a Simulink model using blocks from the Arduino Support Package.

As you can see, we receive the desired motor position from the serial port and compare it to the measured position from the Analog Input. The position error goes through the PID block which generates a voltage to be sent to the motor. We also send the measured position through the serial port.

Hardware response

We ran the model on the target, sent it some commands and logged the data transmitted through the serial port. Here is what it looks like:

So, it does a pretty good job of tracking the reference signal. When we compare the response from the hardware to that of the simulation, we observe that they are very close! The System Identification Toolbox model is quite good.

Now it is your turn

How are you designing controllers when a system is difficult to model? Have you tried the Run on Target Hardware capability in Simulink to run your models on the supported boards? Let us know by leaving a comment here.

Brett's Pick this week is "setname", by Andrew Bliss.

Some of you may be familiar with my own togglefig function; it allows one to recall and reuse a figure by specifying its name. I find it extremely practical, and I use it almost on a daily basis. (It's particularly useful when you're iterating on a block of code separated into "sections." Andrew's setname makes a very nice companion to togglefig.)

As of this writing, Andrew has 10 files on the Exchange, which have collectively been downloaded more than 500 times in the past 30 days. (Not too shabby!) Nonetheless, I don't believe any of his submissions have previously been recognized as a Pick of the Week. Until now.

Ironically, the file of Andrew's that I am most enamored of is his least-downloaded.

So what does Andrew's function do? If you've previously indicated a title for the axes in your figure, setname automatically renames the figure to match axes name. If you haven't titled your graphic, the function extracts a string from the x- and y- labels of the axes to construct a figure name.

Consider, for example, this snippet of code:

t = 0:pi/64:8*pi;
plot(t,sin(t)+rand(size(t)))
xlabel('Time')
ylabel('Noisy Sine')
title('Noisy Sinusoid')

This creates a figure not-so-usefully named "Untitled" (with a number appended to differentiate it from the other n "Untitled" figures you've previously created.)

Now issue the setname command (no arguments needed--it works by default on the current axes, though you may specify a non-default axes.) Now the figure is renamed to "Noisy Sinusoid," and you can readily pick it out from all the other open figures.

Alternatively, if you didn't explicitly provide a title, setname would automatically name the figure "Time vs Noisy Sine." Now not only can I find the figure, I can easily reuse it:

% Recall the previous figure...
togglefig('Noisy Sinusoid')
plot(t,cos(t)+rand(size(t)))
title('Noisy Cosine')
% ...and then rename the figure!
setname

This is brilliant, Andrew...thanks! One suggestion: it would be really nice if I could pass to setname an argument ('all'?) that tells the function to rename all open figures. (For those times when I've lazily created 20 different 'Untitled' figures.)

As always, I welcome your thoughts and comments. Or leave feedback for Andrew here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

One of Jim Wilkinson's outstanding skills was his ability to pick great examples. A previous post featured his signature polynomial. This post features his signature matrices.

Contents

Tridiagonal Matrices

When working with eigenvalues, we always have to be concerned about multiplicity and separation. Multiple eigenvalues and clusters of close eigenvalues require careful attention. Symmetric tridiagonal matrices are easier to handle than the general case. If none of the off diagonal elements are zero, there are no multiple eigenvalues. But if none of the off diagonal elements are "small", can we conclude that the eigenvalues are reasonably well separated? One of Wilkinson's matrices answers this question with a resounding "no".

Wilkinson's Matrices

Wilkinson introduces the matrices $W_{2n+1}^-$ and $W_{2n+1}^+$ on page 308 of The Algebraic Eigenvalue Problem and then employs them throughout the rest of the book. These matrices are the families of symmetric tridiagonal matrices generated by the following MATLAB function. ($W_n^+$ is also generated by the function wilkinson in toolbox\matlab\elmat that has been part of MATLAB for many years.)

type generate_wilkinson_matrices

function [Wm,Wp] = generate_wilkinson_matrices(n)

   D = diag(ones(2*n,1),1);
   Wm = diag(-n:n) + D + D';
   Wp = abs(diag(-n:n)) + D + D';

Here, for example, is $n = 3, 2n+1 = 7$.

format compact
[Wm,Wp] = generate_wilkinson_matrices(3)

Wm =
    -3     1     0     0     0     0     0
     1    -2     1     0     0     0     0
     0     1    -1     1     0     0     0
     0     0     1     0     1     0     0
     0     0     0     1     1     1     0
     0     0     0     0     1     2     1
     0     0     0     0     0     1     3
Wp =
     3     1     0     0     0     0     0
     1     2     1     0     0     0     0
     0     1     1     1     0     0     0
     0     0     1     0     1     0     0
     0     0     0     1     1     1     0
     0     0     0     0     1     2     1
     0     0     0     0     0     1     3

Eigenvalues

Wilkinson usually has $2n+1 = 21$. Here are the eigenvalues of $W_{21}^-$ and $W_{21}^+$.

format long
[Wm,Wp] = generate_wilkinson_matrices(10);
disp('       eig(Wm(21))         eig(Wp(21))')
disp(flipud([eig(Wm) eig(Wp)]))

       eig(Wm(21))         eig(Wp(21))
  10.746194182903357  10.746194182903393
   9.210678647333035  10.746194182903322
   8.038941119306445   9.210678647361332
   7.003952002665359   9.210678647304919
   6.000225680185169   8.038941122829023
   5.000008158672943   8.038941115814275
   4.000000205070442   7.003952209528674
   3.000000003808129   7.003951798616375
   2.000000000054486   6.000234031584166
   1.000000000000619   6.000217522257097
   0.000000000000003   5.000244425001915
  -1.000000000000618   4.999782477742903
  -2.000000000054490   4.004354023440857
  -3.000000003808127   3.996048201383625
  -4.000000205070437   3.043099292578824
  -5.000008158672947   2.961058884185726
  -6.000225680185168   2.130209219362506
  -7.003952002665363   1.789321352695084
  -8.038941119306442   0.947534367529292
  -9.210678647333047   0.253805817096678
 -10.746194182903357  -1.125441522119985

Properties

Here are some important properties of these eigenvalues.

Remarkable Closeness

All of the off diagonal elements in $W_{21}^+$ equal one, so none are "small", and yet the first two eigenvalues of $W_{21}^+$ agree to fourteen significant figures. Wilkinson calls this "remarkable" and "pathological". The relative distance is

format short e
e = flipud(eig(Wp));
delta = (e(1)-e(2))/e(1)

delta =
   6.6120e-15

Eigenvectors

The behavior of the eigenvalues can be understood by looking at the eigenvectors. Here is the vector $u_1$ associated with the first eigenvalue of $W_{21}^-$ and the vectors $v_1$ and $v_2$ associated with the leading pair of eigenvalues of $W_{21}^+$, Each vector has been normalized so that its first component is one.

[U,~] = eig(Wm);
[V,~] = eig(Wp);
u1 = flipud(U(:,end)/U(end,end));
v1 = flipud(V(:,end)/V(end,end));
v2 = flipud(V(:,end-1)/V(end,end-1));

format long
[u1 v1 v2]

ans =
   1.000000000000000   1.000000000000000   1.000000000000000
   0.746194182903358   0.746194182903392   0.746194182903321
   0.302999941502167   0.302999941502254   0.302999941502075
   0.085902493869951   0.085902493870164   0.085902493869724
   0.018807481330335   0.018807481331049   0.018807481329571
   0.003361464615150   0.003361464618322   0.003361464611748
   0.000508147087273   0.000508147104788   0.000508147068490
   0.000066594309069   0.000066594424058   0.000066594185758
   0.000007705362252   0.000007706235465   0.000007704425844
   0.000000798285440   0.000000805807738   0.000000790218740
   0.000000074882661   0.000000147323229  -0.000000002800555
   0.000000006418173   0.000000777356287  -0.000000820314052
   0.000000000506441   0.000007428942095  -0.000007992139484
   0.000000000037026   0.000064197613845  -0.000069080489810
   0.000000000002522   0.000489858240829  -0.000527118748833
   0.000000000000161   0.003240481200881  -0.003486964947267
   0.000000000000010   0.018130575985481  -0.019509658947141
   0.000000000000001   0.082810753074099  -0.089109664858083
   0.000000000000000   0.292094585462557  -0.314312449184672
   0.000000000000000   0.719337698380754  -0.774053354706959
   0.000000000000000   0.964008718993031  -1.037335016061421

And here is a plot of these three eigenvectors. The components of the first half of each of the three vectors decay rapidly and all three agree with each other to many decimal places.

n = 10;
plot_wilkinson_eigenvectors

For a tridiagonal matrix $A$, the eigenvalue/vector equation

$$ Ax = \lambda x $$

is a three-term recurrence relation for the components of $x$. For the first half of each of the three vectors, we have the same positive elements on the diagonal of $A$, ones on the off diagonal, and nearly the same positive value of $\lambda$, so the recurrence relation is nearly the same and leads to nearly the same rapidly decreasing components.

Once the recurrence relation gets half way through the vector and the components get small, this argument breaks down. For $u_1$ the diagonal elements are negative and so the components of the vector remain small. Perron-Frobenius implies that $v_1$ goes positive. Symmetry implies that $v_2$ goes negative to stay orthogonal to $v_1$.

Perturbation Theory

A closer examination of the elements of $v_1$ reveals that they decay exponentially until they reach a minimum half way through the vector at index $n+1$. It is possible to estimate a priori that this minimum value is of order $O(1/n!)$

format short e
minimum = v1(n+1)
estimate = 1/factorial(n)

minimum =
   1.4732e-07
estimate =
   2.7557e-07

This information about the shape of the dominant eigenvector makes it possible to prove that the corresponding eigenvalue, $\lambda_1$, is a double eigenvalue with separation of order the square of the minimum in the eigenvector, which is $O(1/(n!)^2)$.

separation = e(1)-e(2)
estimate = 1/factorial(n)^2

separation =
   7.1054e-14
estimate =
   7.5941e-14

Acknowledgement

Thanks to Pete Stewart, University of Maryland, for some helpful email while I was working on this post.

Reference

J. H. Wilkinson, The Algebraic Eigenvalue Problem, Claredon Press, Oxford, 1965, 662 pp.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

This is the fourth and last part of my plan (my evil plan?) to rewrite an Image Processing Toolbox example from 20 years ago using more modern MATLAB language features. I got the idea from Dave Garrison's recent article on writing MATLAB apps.

Here's the old app I'm trying to reinvent:

And here's what I had working last time:

There are two things I want to wrap up before calling it good enough:

First let's get the DCT coefficient mask display working. Recall that last time I added a function to compute the reconstructed image, given the desired number of DCT coefficients. I'll add a second output argument to that function in order to return the DCT coefficient mask. Here's the code. The only changes are on the first line (to define the additional output argument) and the last three lines (to compute the mask).

function [I2,mask] = reconstructImage(I,n)
% Reconstruct the image from n of the DCT coefficients in each 8-by-8
% block. Select the n coefficients with the largest variance across the
% image. Second output argument is the 8-by-8 DCT coefficient mask.

% Compute 8-by-8 block DCTs.
f = @(block) dct2(block.data);
A = blockproc(I,[8 8],f);

% Compute DCT coefficient variances and decide
% which to keep.
B = im2col(A,[8 8],'distinct')';
vars = var(B);
[~,idx] = sort(vars,'descend');
keep = idx(1:n);

% Zero out the DCT coefficients we are not keeping.
B2 = zeros(size(B));
B2(:,keep) = B(:,keep);

% Reconstruct image using 8-by-8 block inverse
% DCTs.
C = col2im(B2',[8 8],size(I),'distinct');
finv = @(block) idct2(block.data);
I2 = blockproc(C,[8 8],finv);

mask = false(8,8);
mask(keep) = true;
end

Next I need some code to visualize the coefficient mask. I want to display it as image with gray lines drawn between the mask pixels. So I added a local function called displayCoefficientMask:

function displayCoefficientMask(mask,ax)
imshow(mask,'Parent',ax)
for k = 0.5:1.0:8.5
    line('XData',[0.5 8.5], ...
        'YData',[k k], ...
        'Color',[0.6 0.6 0.6], ...
        'LineWidth',2, ...
        'Clipping','off', ...
        'Parent',ax);
    line('XData',[k k],...
        'YData',[0.5 8.5],...
        'Color',[0.6 0.6 0.6], ...
        'LineWidth',2,...
        'Clipping','off', ...
        'Parent',ax);
end
title(ax,'DCT Coefficient Mask')
end

The last step is to call displayCoefficientMask from the update method (which gets called whenever the slider is moved). In the code below, I have modified the call to reconstructImage to use two output arguments in order to get the mask; I have assigned the various app properties; and I have added the call to displayCoefficientMask at the end.

function update(app)
% Update the computation
[recon_image,mask] = reconstructImage(app.OriginalImage, ...
    app.NumDCTCoefficients);

diff_image = imabsdiff(app.OriginalImage, recon_image);

% Update the app properties
app.ReconstructedImage = recon_image;
app.ErrorImage = diff_image;
app.DCTCoefficientMask = mask;

% Update the display
imshow(app.OriginalImage,'Parent',app.OriginalImageAxes);
title(app.OriginalImageAxes,'Original Image');

imshow(recon_image,'Parent',app.ReconstructedImageAxes);
title(app.ReconstructedImageAxes,'Reconstructed Image');

imshow(diff_image,[],'Parent',app.ErrorImageAxes);
title(app.ErrorImageAxes,'Error Image');

displayCoefficientMask(mask,app.MaskAxes);

drawnow;
end

Here's the result with the DCT coefficient mask visualization included:

The last thing I want to do with this little app is to allow users to set the app's NumDCTCoefficients property from the command line and to have the app automatically update. To do this, I'll make a couple of changes to the NumDCTCoefficients property. First, I'll make it a dependent property. Instead of being stored independently, this property will be computed on-the-fly from slider setting whenever it is queried. That requires that I define a property get method that computes the property's value on demand. And last I'll need a property set method that defines what actions should be taken whenever the user sets the property.

Here's the modified property block that indicates that NumDCTCoefficients is a dependent property.

properties (Dependent)
    NumDCTCoefficients
end

And here are the get and set functions for NumDCTCoefficients. The get function computes the number on-the-fly based on the current slider position. The set function modifies the slider position and then calls the update method to recompute and redisplay everything.

function value = get.NumDCTCoefficients(app)
    value = round(get(app.Slider,'Value') * 64);
end

function set.NumDCTCoefficients(app,num_coefficients)
    set(app.Slider,'Value',num_coefficients/64)
    update(app);
end

Here's an example of this interaction.

app = dctCompressionExample_v5

app = 

  dctCompressionExample_v5 with properties:

         OriginalImage: [240x240 double]
    ReconstructedImage: [240x240 double]
            ErrorImage: [240x240 double]
    DCTCoefficientMask: [8x8 logical]
    NumDCTCoefficients: 3

I can see the DCT coefficient mask on the app, but I can also look at the DCTCoefficientMask property.

app.DCTCoefficientMask

ans =

     1     1     0     0     0     0     0     0
     1     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0

And I can set the NumDCTCoefficients property, which causes the app to update.

app.NumDCTCoefficients = 1;

OK, I think that's enough to get all the basic ideas. If you want to play around with the final version of the code for this blog post, you can download it from here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Jiro's pick this week is Sim.I.am by Jean-Pierre de la Croix.

This week, I had a trip to University of Cincinnati with my colleague (and fellow blogger) Sean de Wolski. There, I did a session titled "Enabling Project-Based Learning with MATLAB, Simulink, and Hardware". Project-based learning is a "collaborative education style facilitated by teachers, aimed at increasing students' retention of content in a way that is directly engaging, through projects applicable to life outside of the classroom." There are various tools that would serve the purpose of engaging students, and MATLAB and Simulink are certainly among them.

With MATLAB and Simulink, you can have software-based projects (simulations) and hardware-based projects. I wrote about some of the hardware support that is included with Simulink, and this allows you to easily take what you have in simulation onto low cost hardware. However, when you are trying to implement a project-based learning approach for a class of hundreds of students, purchasing hardware devices for every student becomes impractical. In that case, simulation becomes even more important.

Jean-Pierre is a Ph.D. student in GRITSLab at Georgia Tech (my alma mater!). Sim.I.am is a simulator that can be used to learn and apply control theory for mobile robots. The robot used in the simulator is based on the Khepera III (K3) mobile robot. This File Exchange submission comes with a detailed manual for the simulator, although the App is user-friendly and self-explanatory. What I like even more is the collection of programming exercises included with the package. It's ready to be used in a robotics course right away!

In this animation, you can see a target-seeking algorithm in action. The green dot indicates the target that the robot seeks.

At my session at University of Cincinnati, I also talked about a similar MATLAB-based Simulator for the iRobot® Create®. This simulator works alongside the MATLAB-iRobot Create interface, and this allows students to test their control algorithms they developed in simulation on the actual robot.

You can find other project-based learning materials created by people from around the world in the Classroom Resources page. Check out the Hardware for Project-Based Learning page to see what kind of hardware platforms interface with MATLAB and Simulink.

Comments

Are you involved with project-based learning? Let us know about it here and leave a comment for Jean-Pierre.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Even if we add more targets and support more features for each target with every new release, it is still possible that you will need a driver that is not included in the Simulink support package. If you desperately need a driver for your hardware, you can always build it yourself.

Let's see how this works using an example from a LEGO NXT project I worked on.

The LEGO Light Sensor

The LEGO NXT kit ships with a Light Sensor including a red LED that can be turned on or off. In the Simulink support package, you can control the light from a checkbox parameter in the block dialog.

For our project, we used this sensor to make a line following robot. To make our algorithm more robust, we thought it would be interesting to turn the light on and off while tracking the line, to actively filter out the ambient light.

Step 1: Determine the code to be generated for your driver block

In the LEGO MINDSTORMS NXT Support from Simulink, you can find many small examples showing how the robot can be programmed in C. By default, they are located in C:\MATLAB\SupportPackages\R2013a\nxtOSEK\samples_c. In one of those example, I found a function initializing the robot that looked like:

The function ecrobot_set_light_sensor_active is exactly what I need to generate from my driver block. I also figured that a corresponding function ecrobot_set_light_sensor_inactive exists to turn the light off, and I found that those functions are declared in a file named ecrobot_interface.h.

Step 2 - Create an S-Function

Here you have 2 options to create the S-Function

Step 2 - Option 1: S-Function Builder Block

My colleague Giampiero Campa published a very good submission on MATLAB Central titled Device Drivers showing how to use the S-Function Builder block to include the code you found during step 1 in your model.

His submission contains detailed procedures and screen captures to guide you through the process step-by-step. If you are intimidated by writing an S-function, I recommend using his S-Function Builder technique.

Step 2 - Option 2: Writing an S-Function and a TLC file

If you are like me and want to understand the magic that is happening when you click build in the S-Function Builder block, this second option is for you.

First, we need to realize that unless you model the interaction of the sensor/actuator with the environment, in simulation, driver blocks typically do nothing. All blocks must specify the number of ports and parameters, and their dimensions, even if they do nothing. In my case, I created an s-function with 1 input port of dimension 1 to specify if the light should be on or off. Here is the entire code for my s-function.

To specify the code generated for your block, you need to use the Target Language Compiler.

Concretely, this means that you need to write a TLC file for your block. For that, I recommend starting with examples from sfundemos.

For this example we need our TLC to do two things: Tell the compiler to include ecrobot_interface.h, and call ecrobot_set_light_sensor_active and ecrobot_set_light_sensor_inactive. Here is what it looks like.

Note that I used BlockTypeSetup to include the header file, and Outputs to define what the block output method should be.

Now it's your turn

Download this LEGO example or a similar example for the Arduino target and begin creating your own driver blocks!

If you develop custom drivers for the Simulink Target Hardware, share them on MATLAB Central and let us know by leaving a comment here.

It's my pleasure to introduce guest blogger Kiran Kintali. Kiran is the product development lead for HDL Coder at MathWorks. In this post, Kiran introduces a new capability in HDL Coder™ that generates synthesizable VHDL/Verilog code directly from MATLAB and highlights some of the key features of this new MATLAB based workflow.

Contents

Introduction to HDL Code Generation from MATLAB

If you are using MATLAB to model digital signal processing (DSP) or video and image processing algorithms that eventually end up in FPGAs or ASICs, read on...

FPGAs provide a good compromise between general purpose processors (GPPs) and application specific integrated circuits (ASICs). GPPs are fully programmable but are less efficient in terms of power and performance; ASICs implement dedicated functionality and show the best power and performance characteristics, but require extremely expensive design validation and implementation cycles. FPGAs are also used for prototyping in ASIC workflows for hardware verification and early software development.

Due to the order of magnitude performance improvement when running high-throughput, high-performance applications, algorithm designers are increasingly using FPGAs to prototype and validate their innovations instead of using traditional processors. However, many of the algorithms are implemented in MATLAB due to the simple-to-use programming model and rich analysis and visualization capabilities. When targeting FPGAs or ASICs these MATLAB algorithms have to be manually translated to HDL.

For many algorithm developers who are well-versed with software programming paradigms, mastering the FPGA design workflow is a challenge. Unlike software algorithm development, hardware development requires them to think parallel. Other obstacles include: learning the VHDL or Verilog language, mastering IDEs from FPGA vendors, and understanding esoteric terms like "multi-cycle path" and "delay balancing".

In this post, I describe an easier path from MATLAB to FPGAs. I will show how you can automatically generate HDL code from your MATLAB algorithm, implement the HDL code on an FPGA, and use MATLAB to verify your HDL code.

MATLAB to Hardware Workflow

The process of translating MATLAB designs to hardware consists of the following steps:

There are some unique challenges in translating MATLAB to hardware. MATLAB code is procedural and can be highly abstract; it can use floating-point data and has no notion of time. Complex loops can be inferred from matrix operations and toolbox functions.

Implementing MATLAB code in hardware involves:

HDL Coder™ simplifies the above tasks though workflow automation.

Example MATLAB Algorithm

Let’s take a MATLAB function implementing histogram equalization and go through this workflow. This algorithm, implemented in MATLAB, enhances image contrast by transforming the values in an intensity image so that the histogram of the output image is approximately flat.

type mlhdlc_heq.m

% Histogram Equalization Algorithm
function [pixel_out] = mlhdlc_heq(x_in, y_in, pixel_in, width, height)

persistent histogram
persistent transferFunc
persistent histInd
persistent cumSum

if isempty(histogram)
    histogram = zeros(1, 2^8);
    transferFunc = zeros(1, 2^8);
    histInd = 0;
    cumSum = 0;
end

% Figure out indices based on where we are in the frame
if y_in < height && x_in < width % valid pixel data
    histInd = pixel_in + 1;
elseif y_in == height && x_in == 0 % first column of height+1
    histInd = 1;
elseif y_in >= height % vertical blanking period
    histInd = min(histInd + 1, 2^8);
elseif y_in < height % horizontal blanking - do nothing
    histInd = 1;
end

%Read histogram
histValRead = histogram(histInd);

%Read transfer function
transValRead = transferFunc(histInd);

%If valid part of frame add one to pixel bin and keep transfer func val
if y_in < height && x_in < width
    histValWrite = histValRead + 1; %Add pixel to bin
    transValWrite = transValRead; %Write back same value
    cumSum = 0;
elseif y_in >= height %In blanking time index through all bins and reset to zero
    histValWrite = 0;
    transValWrite = cumSum + histValRead;
    cumSum = transValWrite;
else
    histValWrite = histValRead;
    transValWrite = transValRead;
end

%Write histogram
histogram(histInd) = histValWrite;

%Write transfer function
transferFunc(histInd) = transValWrite;

pixel_out = transValRead;

Example MATLAB Test Bench

Here is the test bench that verifies that the algorithm works with an example image. (Note that this testbench uses Image Processing Toolbox functions for reading the original image and plotting the transformed image after equalization.)

type mlhdlc_heq_tb.m

%% Test bench for Histogram Equalization Algorithm
clear mlhdlc_heq;
testFile = 'office.png';
RGB = imread(testFile);

% Get intensity part of color image
YCBCR = rgb2ycbcr(RGB);
imgOrig = YCBCR(:,:,1);

[height, width] = size(imgOrig);
imgOut = zeros(height,width);
hBlank = 20;
% make sure we have enough vertical blanking to filter the histogram
vBlank = ceil(2^14/(width+hBlank));

for frame = 1:2
    disp(['working on frame: ', num2str(frame)]);
    for y_in = 0:height+vBlank-1
        %disp(['frame: ', num2str(frame), ' of 2, row: ', num2str(y_in)]);
        for x_in = 0:width+hBlank-1
            if x_in < width && y_in < height
                pixel_in = double(imgOrig(y_in+1, x_in+1));
            else
                pixel_in = 0;
            end

             [pixel_out] = mlhdlc_heq(x_in, y_in, pixel_in, width, height);

             if x_in < width && y_in < height
                 imgOut(y_in+1,x_in+1) = pixel_out;
             end
         end
     end
 end

% Make color image from equalized intensity image
% Rescale image
imgOut = double(imgOut);
imgOut(:) = imgOut/max(imgOut(:));
imgOut = uint8(imgOut*255);

YCBCR(:,:,1) = imgOut;
RGBOut = ycbcr2rgb(YCBCR);

figure(1)
subplot(2,2,1); imshow(RGB, []);
title('Original Image');
subplot(2,2,2); imshow(RGBOut, []);
title('Equalized Image');
subplot(2,2,3); hist(double(imgOrig(:)),2^14-1);
title('Histogram of original Image');
subplot(2,2,4); hist(double(imgOut(:)),2^14-1);
title('Histogram of equalized Image');

Let's simulate this algorithm to see the results.

mlhdlc_heq_tb

working on frame: 1
working on frame: 2

HDL Workflow Advisor

The HDL Workflow Advisor (see the snapshot below) helps automate the steps and provides a guided path from MATLAB to hardware. You can see the following key steps of the workflow in the left pane of the workflow advisor:

Let's look at each workflow step in detail.

Fixed-Point Conversion

Signal processing applications are typically implemented using floating-point operations in MATLAB. However, for power, cost, and performance reasons, these algorithms need to be converted to use fixed-point operations when targeting hardware. Fixed-point conversion can be very challenging and time-consuming, typically demanding 25 to 50 percent of the total design and implementation time. The automatic floating-point to fixed-point conversion workflow in HDL Coder™ can greatly simplify and accelerate this conversion process.

The floating-point to fixed-point conversion workflow consists of the following steps:

Note that this step is optional. You can skip this step if your MATLAB design is already implemented in fixed-point.

HDL Code Generation

The HDL Code Generation step generates HDL code from the fixed-point MATLAB code. You can generate either VHDL or Verilog code that implements your MATLAB design. In addition to generating synthesizable HDL code, HDL Coder™ also generates various reports, including a traceability report that helps you navigate between your MATLAB code and the generated HDL code, and a resource utilization report that shows you, at the algorithm level, approximately what hardware resources are needed to implement the design, in terms of adders, multipliers, and RAMs.

During code generation, you can specify various optimization options to explore the design space without having to modify your algorithm. In the Design Space Exploration and Optimization Options section below, you can see how you can modify code generation options and optimize your design for speed or area.

HDL Verification

Standalone HDL test bench generation:

HDL Coder™ generates VHDL and Verilog test benches from your MATLAB scripts for rapid verification of generated HDL code. You can customize an HDL test bench using a variety of options that apply stimuli to the HDL code. You can also generate script files to automate the process of compiling and simulating your code in HDL simulators. These steps help to ensure the results of MATLAB simulation match the results of HDL simulation.

HDL Coder™ also works with HDL Verifier to automatically generate two types of cosimulation testbenches:

HDL Synthesis

Apart from the language-related challenges, programming for FPGAs requires the use of complex EDA tools. Generating a bitstream from the HDL design and programming the FPGA can be daunting tasks. HDL Coder™ provides automation here, by creating project files for Xilinx® and Altera® that are configured with the generated HDL code. You can use the workflow steps to synthesize the HDL code within the MATLAB environment, see the results of synthesis, and iterate on the MATLAB design to improve synthesis results.

Design Space Exploration and Optimization Options

HDL Coder™ provides the following optimizations to help you explore the design space trade-offs between area and speed. You can use these options to explore various architectures and trade-offs without having to manually rewrite your algorithm.

Speed Optimizations

Area Optimizations

Best Practices

Now, let's look at few best practices related to writing MATLAB code when targeting FPGAs.

When writing a MATLAB design:

When writing a MATLAB test bench:

Conclusion

HDL Coder™ provides a seamless workflow when you want to implement your algorithm in an FPGA. In this post, I have shown you how to take an image processing algorithm written in MATLAB, convert it to fixed-point, generate HDL code, verify the generated HDL code using the test bench, and finally, synthesize the design and implement it in hardware.

See this article about how one of the HDL Coder customers, FLIR has used MATLAB to HDL workflow to achieve good results. You can also learn more about this workflow using the product examples located here.

We hope this brief introduction to the HDL Coder™ and MATLAB-to-HDL code generation, verification framework has shown how you can quickly get started on implementing your MATLAB designs and target FPGAs. Please let us know in the comments for this post how you might use this new functionality. Or, if you've already tried using HDL Coder™, let us know about your experiences here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

The Problem

In the following model, inside the calib function-call subsystem, the Count signal is connected to a Unit Delay block.

When logging data, we can see that the input and output of the Unit Delay are identical, as if the signal was not delayed!

Even though this might seem strange, this behavior is expected. Let's see why.

The simulation loop

To understand what is happening, you first need to be familiar with the order into which the block methods are executed within the simulation loop. To really understand the simulation loop I recommend the following documentation pages:

• First, you need to know that the equations for blocks are implemented in multiple block methods: Outputs, Update, Derivative, etc.
• To get an overview of the order into which those methods are called: How S-Functions Work
• To get the complete details of all the block methods available: Simulink Engine Interaction with C S-Functions

In our case, we have a model that is fixed-step and discrete, so the simulation loop is very simple and looks like this:

The important thing to understand is that Simulink executes the Outputs method of all blocks, and then the Update method of all blocks.

The execution order

The second thing we need to know is the order in which blocks are executed. For that, we display the block sorted order.

When you display the sorted order, you can see red numbers displayed on the block. The first number is the number of the system, so the function-call subsystem in this model is system number 3. The second number is the order of the block in this system. In this case, the Gain is the first (3:0), Unit Delay the second (3:1), Counter the third (3:2), etc.

The effect of the latch

To see the effect of latching, I generated code for this model and added comments to highlight the fact that the feedback loop causes the input signal to change immediately as the output is written:

Such feedback loops break the assumptions on which the Output-Update logic is built... and can make results difficult to interpret. This is why the Latch is available.

With the latch, a copy of the input signal is made at the beginning of the step and this copy is used for the computations:

Conclusion

In conclusion, we can say that introducing the latch is safer, but costs an additional copy. If you know that your function-call subsystem might be involved in a direct feedback loop, use Inport Latching to avoid surprises.

Now it's you turn

Are you latching the inputs of your function-call subsystems? Let us know the reason why by leaving a comment here.

Build Your Own Autonomous Battlebot

A massive maker space. Cutting-edge software. Expert mentors. And just one week to craft the most devastating robot possible.

by Dr. Paul Kassebaum

MathWorks is teaming up with Artisan’s Asylum to wage a robot battle accelerated by Simulink’s run on target hardware capabilities and the computer controlled manufacturing tools at Artisan’s Asylum, the largest maker space on the east coast. The fight will serve as the grand finale of the Cambridge Science Festival on Sunday April 21st at the Center for Arts at the Armory, and will be free and open to the public. The contest begins on Sunday April 14th at Artisan’s Asylum, where teams will design their robots over the course of a week. To sign up and for more info, have a look at designchallengebattlebots.com.

Sergio Biagioni is working at MathWorks to develop a Simulink framework to aid high school and undergraduate contestants in their algorithm designs. Each robot will be equipped with an Arduino Mega 2560, a line sensor, touch sensors, sonar, and a wireless feed akin to GPS sent from an eye in the sky. The overhead camera will use the Computer Vision System Toolbox to locate each robot and its orientation based on colored stickers applied to the chassis.

Rob Masek, facilities manager of Artisan’s Asylum, is organizing the tournament, drawing on his experience as a (three time) contestant on Comedy Central’s BattleBots and as an organizer of several US FIRST Robotics competitions and a home-brew robot competition called Pound of Pain. I’m confident he’ll be able to replicate his previous successes.

Gui Cavalcanti, the founder of Artisan’s Asylum, is whipping up a robot chassis to help Sergio develop the Simulink framework all of the contestants will use. This prototype robot will be made of the same components available to the contestants, and will be offered up to the contestants to help them jump right into software development before their own chassis are fully operational. While the robots in this competition will have wheels, most of Gui’s robots have more than their fair share of legs.

Stay tuned for more details as the story progresses!

And speaking of robot competitions, have you heard about MathWorks’ latest robot competition in the UK?

There is a new datatype we are playing around with that we hope to make available in an upcoming release and we would like your input beforehand.

Contents

New Datatype in Action

Let me show you the new datatype in action so you can first get a feel for it.

inputData = magic(3)

inputData =
     8     1     6
     3     5     7
     4     9     2

outputValues = dis(inputData);

Let's Examine the Output

outputValues

Why are you asking?
     4     2     3
     1     9     6
     8     7     5

Well, that's a bit strange, isn't it? I wonder what the relationship between inputData and outputValues is. What can we learn about outputValues?

whos outputValues

  Name              Size            Bytes  Class    Attributes

  outputValues      1x1               248  dis                

Well, it's a dis array. Let's look at it again.

outputValues

What's it matter to you?
     5     6     9
     4     3     2
     7     8     1

Say what? Let's check it a few more times.

outputValues
outputValues
outputValues

Who wants to know?
     6     3     9
     5     8     7
     4     2     1
Who wants to know?
     5     2     9
     1     4     8
     7     6     3
Who are you to ask me that?
     5     3     7
     9     8     2
     1     4     6

Hoping you get the double meaning here - the dis array not only mixes up the values of the input for display purposes, but also tries to gently *dis*respect you along the way.

Even though this is a silly class, I'll show you the code so you can see how simple it is to make such a class.

type dis

classdef dis
    %dis dis is a class.
    %   In fact, it's a declasse class.
    
    properties
        Data
    end
    properties (Access=protected)
        Original
    end
    
    properties (Constant)
        Answers = {'Why are you asking?' ,...
            'What''s it matter to you?',...
            'Who are you to ask me that?',...
            'Who wants to know?',...
            'What''s the big deal?'}
    end
    
    methods
        function display(obj)
            disp(dis.Answers{randi(length(dis.Answers),1)})
            obj.Data(:) = obj.Data(randperm(numel(obj.Data)));
            disp(obj.Data)
        end
        function obj = dis(in)
            obj.Original = in;
            obj.Data = reshape(in(randperm(numel(in))),size(in));
        end
    end
    
end

Should We Invest More Resources?

Of course, I could also add some numeric functions like plus to dis, but I didn't take the time, in case you didn't find this possible new MATLAB addition useful. So please share your thoughts with us here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Although this is the time of year that I normally mention some of the new features from the "a" release, today I thought I would focus on a feature added several releases ago that I think deserves more attention: parallel cluster-based image display. I can't remember exactly when this was released; I think it might have been R2011c.

Anyway, it turns out image display speed can be dramatically increased if we exploit the capabilities of today's massive parallel computing clusters. You start by opening a MATLAB pool from your desktop session using the Parallel Computing Toolbox function matlabpool. The Parallel Computing Toolbox has a new cluster profile that makes it easy to set things up for fast image display:

matlabpool image-display

Once the pool is configured, then the Image Processing Toolbox function imshow automatically uses the workers in the pool.

Let's try it with my blog author picture:

(I'm writing this blog post at home, so for a cluster I'll be using my home-brew video game console, which runs on a cluster of 128 TRS-80s.)

imshow

Here's what the display looks like on worker #37:

And here is the display on worker #92:

If you look carefully, you can see that the pixel block is not exactly the same size on the two workers. That's because of the automatic worker load balancing feature. Way cool!

Let's get down to some numbers. Specifically, I was interested to see how well image display performance scales with the number of workers. I used my function timeit to measure the display speed as I varied the number of workers. The plot below is "normalized speed-up." That is, 1 is the image display speed for normal MATLAB, with no use of the Parallel Computing Toolbox.

As you can see, we get almost perfectly linear scaling with the number of workers. (I'm not sure what's going on with that blip on the left. I think it might have been caused by a bad tape unit on worker #25.)

So, Parallel Computing Toolbox users, please give this a try and let us know about your experience!

Before I wrap it up for today, I can't resist giving a hint about future development directions. (Long-time readers know how much MathWorks loves to preannounce features.) We recently got our hands on one of the experimental new superfluidic clusters. Early indications are that we'll be able to get image display speed to scale with the square of the number of workers! Woo-hoo!

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

We were musing here about how common it is to want more than two Y axes on a plot. Checking out the File Exchange, there seem to be several candidates, indicating that this is something at least some people find useful.

Contents

Sample Plot

Here's a sample plot using plotyy that comes with MATLAB.

x = 0:0.01:20;
y1 = 200*exp(-0.05*x).*sin(x);
y2 = 0.8*exp(-0.5*x).*sin(10*x);
plotyy(x,y1,x,y2,'plot');

List of Some Possibilities

In addition to plotyy in MATLAB, here's a list of some of the candidates from the File Exchange.

What are You Plotting with More Y Axes?

I am curious to know what kind of data or results you are plotting so that having multiple y-axes makes a compelling presentation. Let us know here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Did you know that a Trigger functionality is now available with the Time Scope in the DSP System Toolbox (R2013a) ? The Time Scope ships both as a System object and also as a Simulink block with the DSP System Toolbox

Triggers are ubiquitous on real-world hardware oscilloscopes. They are useful to analyze signals,specifically to stabilize a repetitive waveform and search for a occurrence of a particular pattern in the signal.

Trigger Modes:

The Trigger in Time Scope supports 3 different modes of operation for a given trigger type. They are:

• Auto: In this mode, the Time Scope stabilizes the display if the trigger criteria is met or continues to display the incoming signal.
• Normal : The time scope displays the stabilized signal only if the trigger criteria is met otherwise nothing is displayed.
• Once: The time scope displays the signal only for the very first time the trigger event is encountered. (Simulation will still keep running)

Here is an animation showing those different modes:

Let's look at a few more examples

Edge Triggering

Let's create a simple sine wave and display it on the Time Scope.

As soon as you click the Trigger button, the triggers panel open allowing you to configure options like level, type, polarity, etc.

Runt Triggering

What are Runt signals? Runt signals are the signals which typically cross a voltage threshold level but fail to cross a second threshold before re-crossing the first threshold level. Runt signals are frequently encountered in digital circuits.

For example if we look at the following waveform:

To look for Runt signals, Let us enable the Trigger, set the Trigger Mode to Normal, set the trigger type as Runt, adjust the voltage threshold and then the Time Scope detects and displays the presence of any Runt signals in the incoming data stream.

Now it's your turn

Go through the trigger documentation, give this a try and let us know what you think by leaving a comment here.

This video shows a technique for catching this problem. It is important to catch these errors from unexpected inputs immediately, because they can have effects on your code much farther down the line. When this happens, it would be hard to trace back to the cause.

The MathWorks Student Robot Challenge

by Dr. Tanya Morton, MathWorks UK

Earlier this year a colleague and I were bouncing around ideas for ways to do more outreach with the local student community. He pointed behind him and said, “Why don’t we make use of the empty space in the new office to host a robot contest?” An idea was born…

Last Friday, this idea came to fruition when MathWorks hosted a Student Robot Challenge. Nine teams from the University of Cambridge participated in the contest. The challenge was to develop a controller for a LEGO MINDSTORM NXT robot to make it visit a set of 12 locations on an arena in the shortest time possible.

After a briefing, the students started working their way through short tasks designed to bring them up to speed with Simulink, the principles of controller design, and the process of generating code to run on the robot. Within five minutes, you could hear the first robot spinning on the spot.

In the afternoon, the students explored a variety of navigation and control strategies, which they tested using a combination of simulation and hardware testing on their robot.

The day ended with a final contest where the teams competed against each other to navigate to the most locations in the shortest time. Team names were picked out of the hat to decide the running order for the final contest. The first to go was the QED-grad team, an experienced team of PhD engineering students. The team successfully hit all 12 locations in a time of 1 minute and 41.3 seconds. Their robot left each location quickly, but tended to spiral around the next target rather than approach it directly.

The LEGO robot with positioning ball used in the final arena

Five of the next seven teams were not able to complete the course within the 3-minute time limit. One of these teams was QED-’12, a team of first-year undergraduate engineers who had been taught to program LEGO by the QED-grad team. Briefly, it seemed that the youngest team in the competition were going to beat their lab demonstrators; however, their robot stopped at the 11th point because they had forgotten to program in the 12th. An important lesson learned! Only the team of PhD mathematicians on Team CCA came close to the time of the QED-grad team with a time of 1 minute and 44.5 seconds.

The last team to go was MatLads, with its leader Richard Peach, a 3rd year undergraduate engineer. MatLads’ robot moved neatly around the course to knock the QED-grad team off the top spot to win by a mere 1.6 seconds!

Paths taken by the fastest two robots (click to see a larger image)

It was an exciting finale to a fun event. The students enjoyed the project-based learning experience and gained some valuable skills to help them in their future careers.

The Final Leaderboard

• If you would like to program your own LEGO robot, then a good starting point is the LEGO MINDSTORMS NXT Support from MATLAB and Simulink.
• If you would like to encourage more project-based learning in your local university or high-school, then encourage your lecturers to learn about Teaching Mechatronics with MATLAB, Simulink, and Low Cost Hardware.
• See the Student Robot Challenge for further information.

This is the third part of my plan to rewrite an Image Processing Toolbox example from 20 years ago using more modern MATLAB language features. I got the idea from Dave Garrison's recent article on writing MATLAB apps.

Here's the old app I'm trying to reinvent:

And here's what I had working the last time:

The code isn't doing any actual DCT computations yet. The idea is to compute 8-by-8 block DCTs, zero out the least significant DCT coefficients of each block, and then reconstruct the image by computing the inverse DCT of each 8-by-8 block.

Let me show how this will work on the sample image I'm using for the app:

pout = imread('pout.tif');
pout2 = pout(1:240,:);
I = im2double(adapthisteq(pout2));
imshow(I)

First, use blockproc to compute the 8-by-8 2-D DCTs.

f = @(block) dct2(block.data);
A = blockproc(I,[8 8],f);
imshow(A,[])
title('DCT coefficients')

Second, I want to compute the variances of each of the 64 coefficients. That is, compute the variance of the (1,1) DCT coefficients from each block, the variance of the (1,2) DCT coefficients from each block, and so on. I'll use the computed variances to determine which DCT coefficients to keep and which to zero out. To do this, I'll rearrange the elements from the previous step so that the corresponding DCT coefficients are together in the columns of a 64-column matrix.

B = im2col(A,[8 8],'distinct')';
vars = var(B);
plot(vars)
title('Variances of the 64 DCT coefficients')

Next, suppose we want to keep the 3 block DCT coefficients with the highest variances. I'll use sort to figure out which ones to keep.

[~,idx] = sort(vars,'descend');
keep = idx(1:3)

keep =

     1     9     2

For an 8-by-8 block of DCT coefficients, index values 1, 9, and 2 correspond to the upper-left DCT coefficient and the two coefficients just the right and down from it. Let's keep just those coefficients.

B2 = zeros(size(B));
B2(:,keep) = B(:,keep);

Now I'll rearrange the columns back into blocks and reconstruct the image by computing the inverse DCT of each block.

C = col2im(B2',[8 8],size(I),'distinct');
finv = @(block) idct2(block.data);
D = blockproc(C,[8 8],finv);
imshow(D)
title('Image reconstructed from 3 DCT coefficients per block')

This algorithm code into a local function to go at the bottom of my app code file, like this:

function I2 = reconstructImage(I,n)
% Reconstruct the image from n of the DCT coefficients in each 8-by-8
% block. Select the n coefficients with the largest variance across the
% image.

% Compute 8-by-8 block DCTs.
f = @(block) dct2(block.data);
A = blockproc(I,[8 8],f);

% Compute DCT coefficient variances and decide
% which to keep.
B = im2col(A,[8 8],'distinct')';
vars = var(B);
[~,idx] = sort(vars,'descend');
keep = idx(1:n);

% Zero out the DCT coefficients we are not keeping.
B2 = zeros(size(B));
B2(:,keep) = B(:,keep);

% Reconstruct image using 8-by-8 block inverse
% DCTs.
C = col2im(B2',[8 8],size(I),'distinct');
finv = @(block) idct2(block.data);
I2 = blockproc(C,[8 8],finv);
end

I already have an update method in my class, so I can add code to that method that will compute the reconstructed image and error image. Here's the revised update method:

function update(app)
    recon_image = reconstructImage(app.OriginalImage, ...
        app.NumDCTCoefficients);

    diff_image = imabsdiff(app.OriginalImage, recon_image);

    imshow(app.OriginalImage,'Parent',app.OriginalImageAxes);
    title(app.OriginalImageAxes,'Original Image');

    imshow(recon_image,'Parent',app.ReconstructedImageAxes);
    title(app.ReconstructedImageAxes,'Reconstructed Image');

    imshow(diff_image,[],'Parent',app.ErrorImageAxes);
    title(app.ErrorImageAxes,'Error Image');

    imshow(zeros(size(app.OriginalImage)),'Parent',app.MaskAxes);
    title(app.MaskAxes,'DCT Coefficient Mask');
    drawnow;
end

Although I don't have all the interactive behavior wired up yet, and don't have the DCT coefficient mask visualization done yet, the app runs and looks a bit more finished than it did before. It now shows the results for reconstructing the image from 3 DCT coefficients per block.

I think it'll take just one more post to finish making the rest of this app work.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Simulink Variants provide increased functionality over Configurable Subsystems:

• They exist in two flavors: Model Variants and Subsystem Variants.
• They can be controlled programmatically via variables in the MATLAB workspace.
• They are in your model instead of a Simulink Library which simplifies model development.

Here are some tips to help you migrate.

Upgrade Advisor

In R2012b, the Upgrade Advisor includes a check titled Identify configurable subsystem blocks for converting to variant subsystem blocks.

If you run the check and a configurable subsystem is found, the result will tell you that "These blocks can be upgraded to variant subsystems as they have better programmatic control and code generation capabilities."

Context Menu

If you right-click on a configurable subsystem (or any subsystem), you will see a new option to convert the subsystem to a variant.

A window will allow you to specify some details needed for the conversion.

and a new model will pop up with your new block converted to a variant subsystem. Notice the different icon in the lower left corner of the block, it indicates that the variant subsystem is configured to override variant conditions and behaves exactly like the original configurable subsystem.

Now it's your turn

Are you ready to move to Variant Subsystems? Let us know by leaving a comment here

1.       Great problems to draw on

There are already over 900 ready-made MATLAB programming problems, in just the first year online.  Search for your subject matter of interest.  Check back often, as the Community area is growing rapidly.

2.       Fits into your workflow

Navigate to Cody problems by hyperlink directly from your Learning Management System (LMS), HTML syllabus page or simply link directly into each Cody problem in email to your students.

3.       Engaging for the learner

As Jeff points out on the Community Blogs, Cody has an addictive quality you can use to deepen student engagement.

4.       Can use online and/or with MATLAB

Copy the Cody test suite into desktop MATLAB.  Let students write their Cody function with the full power of MATLAB, then paste their answer back into Cody to verify correctness.  For example, the Times 2 problem could be solved in MATLAB as follows:

testsuite.m:
%%
assert(isequal(times2(1),2));
%%
assert(isequal(times2(11),22));
%%
assert(isequal(times2(-3),-6));
%%
assert(isequal(times2(29),58));


times2.m:
function y = times2(x)
% Modify the line below so that the output y is twice the incoming value x
y = 2*x;
% After you fix the code, press the "Submit" button, and you're on your way.
end


Running testsuite.m confirms the function passes each test case, hence success.

5.       Create your own problems

Try creating a Cody problem for yourself.

6.       Experiment to see how people try to solve your problems

Draw on the MATLAB Central community to confirm your Cody problem challenges exactly the skills you want students to learn.  We offer you a large population of MATLAB experts eager to try out the latest MATLAB programming puzzle.  Use the Solution Map to see outlier attempts, to better understand how your students could react before offering the practice problem in class.

7.       Identify learning gaps

Using the solution map, keep an eye across multiple problems from the point in time you’ve offered in class.  Compare these across problems to identify subject matter areas meriting more lecture coverage.  For example,

Problem 7: Column Removal

shows a cluster of incorrect attempts worth exploring.  You’ll quickly see which topics need follow-up coverage in class.

I’ve also prepared a five minute video to demonstrate the above points.  Please let me know how you’ve found Cody helpful.

Try Cody for yourself at www.mathworks.com/matlabcentral/cody, and get your students addicted to learning in MATLAB.

Steven Lord, Andy Campbell, and David Hruska are members of the Quality Engineering group at MathWorks who are guest blogging today to introduce a new feature in R2013a, the MATLAB unit testing infrastructure. There are several submissions on the MATLAB Central File Exchange related to unit testing of MATLAB code. Blogger Steve Eddins wrote one highly rated example back in 2009. In release R2013a, MathWorks included in MATLAB itself a MATLAB implementation of the industry-standard xUnit testing framework.

If you're not a software developer, you may be wondering if this feature will be of any use to you. In this post, we will describe one way someone who may not consider themselves a software developer may be able to take advantage of this framework using the example of a professor grading students' homework submissions. That's not to say that the developers in the audience should move on to the next post; you can use these tools to test your own code just like a professor can use them to test code written by his or her students.

There is a great deal of functionality in this feature that we will not show here. For more information we refer you to the MATLAB Unit Testing Framework documentation.

Contents

Background

In order to use this feature, you should be aware of how to define simple MATLAB classes in classdef files, how to define a class that inherits from another, and how to specify attributes for methods and properties of those classes. The object-oriented programming documentation describes these capabilities.

Problem Statement

As a professor in an introductory programming class, you want your students to write a program to compute Fibonacci numbers. The exact problem statement you give the students is:

Create a function "fib" that accepts a nonnegative integer n and returns
the nth Fibonacci number. The Fibonacci numbers are generated by this
relationship:

F(0) = 1
F(1) = 1
F(n) = F(n-1) + F(n-2) for integer n > 1

Your function should throw an error if n is not a nonnegative integer.

Basic Unit Test

The most basic MATLAB unit test is a MATLAB classdef class file that inherits from the matlab.unittest.TestCase class. Throughout the rest of this post we will add additional pieces to this basic framework to increase the capability of this test and will change its name to reflect its increased functionality.

dbtype basicTest.m

1     classdef basicTest < matlab.unittest.TestCase
2         
3     end

test = basicTest

test = 
  basicTest with no properties.

Running a Test

To run the test, we can simply pass test to the run function. There are more advanced ways that make it easier to run a group of tests, but for our purposes (checking one student's answer at a time) this will be sufficient. When you move to checking multiple students' answers at a time, you can use run inside a for loop.

Since basicTest doesn't actually validate the output from the student's function, it doesn't take very long to execute.

results = run(test)

results = 
  0x0 TestResult array with properties:

    Name
    Passed
    Failed
    Incomplete
    Duration
Totals:
   0 Passed, 0 Failed, 0 Incomplete.
   0 seconds testing time.

Let's say that a student named Thomas submitted a function fib.m as his solution to this assignment. Thomas's code is stored in a sub-folder named thomas. To set up our test to check Thomas's answer, we add the folder holding his code to the path.

addpath('thomas');
dbtype fib.m

1     function y = fib(n)
2     if n <= 1
3         y = 1;
4     else
5         y = fib(n-1)+fib(n-2);
6     end

Test that F(0) Equals 1

The basicTest is a valid test class, and we can run it, but it doesn't actually perform any validation of the student's test file. The methods that will perform that validation need to be written in a methods block that has the attribute Test specified.

The matlab.unittest.TestCase class includes qualification methods that you can use to test various qualities of the results returned by the student files. The qualification method that you will likely use most frequently is the verifyEqual method, which passes if the two values you pass into it are equal and reports a test failure if they are not.

The documentation for the matlab.unittest.TestCase class lists many other qualification methods that you can use to perform other types of validation, including testing the data type and size of the results; matching a string result to an expected string; testing that a given section of code throws a specific errors or issues a specific warning; and many more.

This simple test builds upon generalTest by adding a test method that checks that the student's function returns the value 1 when called with the input 0.

dbtype simpleTest.m

1     classdef simpleTest < matlab.unittest.TestCase
2         methods(Test)
3             function fibonacciOfZeroShouldBeOne(testCase)
4                 % Evaluate the student's function for n = 0
5                 result = fib(0);
6                 testCase.verifyEqual(result, 1);
7             end
8         end
9     end

Thomas's solution to the assignment satisfies this basic check. We can use the results returned from run to display the percentage of the tests that pass.

results = run(simpleTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running simpleTest
.
Done simpleTest
__________

results = 
  TestResult with properties:

          Name: 'simpleTest/fibonacciOfZeroShouldBeOne'
        Passed: 1
        Failed: 0
    Incomplete: 0
      Duration: 0.0112
Totals:
   1 Passed, 0 Failed, 0 Incomplete.
   0.011168 seconds testing time.
100% Passed.

Test that F(pi) Throws an Error

Now that we have a basic positive test in place we can add in a test that checks the behavior of the student's function when passed a non-integer value (like n = pi) as input. The assignment stated that when called with a non-integer value, the student's function should error. Since the assignment doesn't require a specific error to be thrown, the test passes as long as fib(pi) throws any exception.

dbtype errorCaseTest.m

1     classdef errorCaseTest < matlab.unittest.TestCase
2         methods(Test)
3             function fibonacciOfZeroShouldBeOne(testCase)
4                 % Evaluate the student's function for n = 0
5                 result = fib(0);
6                 testCase.verifyEqual(result, 1);
7             end
8             function fibonacciOfNonintegerShouldError(testCase)
9                 testCase.verifyError(@()fib(pi), ?MException);
10            end
11        end
12    end

Thomas forgot to include a check for a non-integer valued input in his function, so our test should indicate that by reporting a failure.

results = run(errorCaseTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running errorCaseTest
.
================================================================================
Verification failed in errorCaseTest/fibonacciOfNonintegerShouldError.

    ---------------------
    Framework Diagnostic:
    ---------------------
    verifyError failed.
    --> The function did not throw any exception.
        
        Expected Exception Type:
            MException
    
    Evaluated Function:
            @()fib(pi)

    ------------------
    Stack Information:
    ------------------
    In C:\Program Files\MATLAB\R2013a\toolbox\matlab\testframework\+matlab\+unittest\+qualifications\Verifiable.m (Verifiable.verifyError) at 637
    In H:\Documents\LOREN\MyJob\Art of MATLAB\errorCaseTest.m (errorCaseTest.fibonacciOfNonintegerShouldError) at 9
================================================================================
.
Done errorCaseTest
__________

Failure Summary:

     Name                                            Failed  Incomplete  Reason(s)
    =============================================================================================
     errorCaseTest/fibonacciOfNonintegerShouldError    X                 Failed by verification.
    
results = 
  1x2 TestResult array with properties:

    Name
    Passed
    Failed
    Incomplete
    Duration
Totals:
   1 Passed, 1 Failed, 0 Incomplete.
   0.026224 seconds testing time.
50% Passed.

Another student, Benjamin, checked for a non-integer value in his code as you can see on line 2.

rmpath('thomas');
addpath('benjamin');
dbtype fib.m

1     function y = fib(n)
2     if (n ~= round(n)) || n < 0
3         error('N is not an integer!');
4     elseif n == 0 || n == 1
5         y = 1;
6     else
7         y = fib(n-1)+fib(n-2);
8     end

Benjamin's code passed both the test implemented in the fibonacciOfZeroShouldBeOne method (which we copied into errorCaseTest from simpleTest) and the new test case implemented in the fibonacciOfNonintegerShouldError method.

results = run(errorCaseTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running errorCaseTest
..
Done errorCaseTest
__________

results = 
  1x2 TestResult array with properties:

    Name
    Passed
    Failed
    Incomplete
    Duration
Totals:
   2 Passed, 0 Failed, 0 Incomplete.
   0.010132 seconds testing time.
100% Passed.

Basic Test for Students, Advanced Tests for Instructor

The problem statement given earlier in this post is a plain text description of the homework assignment we assigned to the students. We can also state the problem for the students in code (if they're using release R2013a or later) by giving them a test file they can run just like simpleTest or errorCaseTest. They can directly use this "requirement test" to ensure their functions satisfy the requirements of the assignment.

dbtype studentTest.m

1     classdef studentTest < matlab.unittest.TestCase
2         methods(Test)
3             function fibonacciOfZeroShouldBeOne(testCase)
4                 % Evaluate the student's function for n = 0
5                 result = fib(0);
6                 testCase.verifyEqual(result, 1);
7             end
8             function fibonacciOfNonintegerShouldError(testCase)
9                 testCase.verifyError(@()fib(pi), ?MException);
10            end
11        end
12    end

In order for the student's code to pass the assignment, it will need to pass the test cases given in the studentTest unit test. However, we don't want to use studentTest as the only check of the student's code. If we did, the student could write their function to cover only the test cases in the student test file.

We could solve this problem by having two separate test files, one containing the student test cases and one containing additional test cases the instructor uses in the grading process. Can we avoid having to run both test files manually or duplicating the code from the student test cases in the instructor test? Yes!

To do so, we write an instructor test file to incorporate, through inheritance, the student test file. We can then add additional test cases to the instructor test file. When we run this test it should run three test cases; two inherited from studentTest, fibonacciOfZeroShouldBeOne and fibonacciOfNonintegerShouldError, and one from instructorTest itself, fibonacciOf5.

dbtype instructorTest.m

1     classdef instructorTest < studentTest
2         % Because the student test file is a matlab.unittest.TestCase and
3         % instructorTest inherits from it, instructorTest is also a
4         % matlab.unittest.TestCase.
5         
6         methods(Test)
7             function fibonacciOf5(testCase)
8                 % Evaluate the student's function for n = 5
9                 result = fib(5);
10                testCase.verifyEqual(result, 8, 'Fibonacci(5) should be 8');
11            end
12        end
13    end

Let's look at Eric's test file that passes the studentTestFile test, but in which he completely forgot to implement the F(n) = F(n-1)+F(n-2) recursion step.

rmpath('benjamin');
addpath('eric');
dbtype fib.m

1     function y = fib(n)
2     if (n ~= round(n)) || n < 0
3         error('N is not an integer!');
4     end
5     y = 1;

It should pass the student unit test.

results = run(studentTest);
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running studentTest
..
Done studentTest
__________

100% Passed.

It does NOT pass the instructor unit test because it fails one of the test cases.

results = run(instructorTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running instructorTest
..
================================================================================
Verification failed in instructorTest/fibonacciOf5.

    ----------------
    Test Diagnostic:
    ----------------
    Fibonacci(5) should be 8

    ---------------------
    Framework Diagnostic:
    ---------------------
    verifyEqual failed.
    --> NumericComparator failed.
        --> The values are not equal using "isequaln".
    
    Actual Value:
             1
    Expected Value:
             8

    ------------------
    Stack Information:
    ------------------
    In C:\Program Files\MATLAB\R2013a\toolbox\matlab\testframework\+matlab\+unittest\+qualifications\Verifiable.m (Verifiable.verifyEqual) at 411
    In H:\Documents\LOREN\MyJob\Art of MATLAB\instructorTest.m (instructorTest.fibonacciOf5) at 10
================================================================================
.
Done instructorTest
__________

Failure Summary:

     Name                         Failed  Incomplete  Reason(s)
    ==========================================================================
     instructorTest/fibonacciOf5    X                 Failed by verification.
    
results = 
  1x3 TestResult array with properties:

    Name
    Passed
    Failed
    Incomplete
    Duration
Totals:
   2 Passed, 1 Failed, 0 Incomplete.
   0.028906 seconds testing time.
66.6667% Passed.

Benjamin, whose code we tested above, wrote a correct solution to the homework problem.

rmpath('eric');
addpath('benjamin');

results = run(instructorTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

rmpath('benjamin');

Running instructorTest
...
Done instructorTest
__________

results = 
  1x3 TestResult array with properties:

    Name
    Passed
    Failed
    Incomplete
    Duration
Totals:
   3 Passed, 0 Failed, 0 Incomplete.
   0.015946 seconds testing time.
100% Passed.

Conclusion

In this post, we showed you the basics of using the new MATLAB unit testing infrastructure using homework grading as a use case.

We checked that the student's code worked (by returning the correct answer) for one valid value and worked (by throwing an error) for one invalid value. We also showed how you can use this infrastructure to provide an aid/check for the students that you can also use as part of your grading.

We hope this brief introduction to the unit testing framework has shown you how you can make use of this feature even if you don't consider yourself a software developer. Let us know in the comments for this post how you might use this new functionality. Or, if you've already tried using matlab.unittest, let us know about your experiences here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

R2013a, the latest semi-annual MathWorks product release, just went live. One of the significant new capabilities in the MATLAB release is a new unit test framework (overview video, documentation). I am very happy to see this!

Long-time blog readers might be wondering, though, what about MATLAB xUnit? This is a unit test framework that I created and put on the File Exchange. I first wrote about it in this space in 2009. It is my second most popular File Exchange contribution (over 200 downloads in the last 30 days).

Today I want to recap how MATLAB xUnit came to be, explain what will probably happen to it, and convince its users to give the new "official" unit test framework in R2013a a try.

I learned to write unit tests when I came to MathWorks in 1993 as a new software developer. I've used several generations of test frameworks created for internal use here. Until recently, our internal testing machinery was tightly coupled to the entire system we used to build, test, and release our products, and so it was impractical to use it for small projects. This became an issue for me sometime around 2002-2003, as I was working on the first edition of Digital Image Processing Using MATLAB. As the "software guy" among the three coauthors, I was responsible for overseeing the MATLAB functions in the book, including making sure that everything was working OK. So I wrote a simple test harness that exercised the book's examples and some of the functions.

Later, the MATLAB R2008a release included the new generation of object-oriented language features. I was interested in learning more about it. Coincidentally, I had just read Kent Beck's Test-Driven Development, which included a case study on using test-driven development methods to create an xUnit-style test harness. I decided that would be an excellent project with which to learn about test-driven development as well as the R2008a MATLAB language changes. (With the wisdom of hindsight, I have to say here that I don't actually recommend learning test-driven development by using test-driven development to develop a test harness. It hurt my brain too much.)

Anyway, sometime in the spring on 2008 I had something pretty basic working. At about that time, Greg Wilson (original creator of Software Carpentry) visited MathWorks to give a talk about his edited book, Beautiful Code. I met Greg then, and we talked about the needs of MATLAB scientific users. Greg encouraged me to turn my fledgling tool into something real. He also invited me to write up something for the IEEE magazine, Computing in Science and Engineering. That became the article "Automated Software Testing for MATLAB", published in the November/December 2009 issue.

The MATLAB xUnit package raised the visibility of unit testing tools in MATLAB, and its popularity helped make the case for putting something like it in MATLAB. At the same time, our testing tools team was looking to modernize our testing infrastructure. They wanted to take advantage of the latest test framework design principles, and they wanted to decouple the testing machinery from the rest of our complex systems for building and releasing products. So an idea was born: let's make a test framework that can stand on its own, that can test our own software, that can test customer code, and that we can ship in MATLAB!

And now, here it finally is!

Here are a few examples that I've taken from the documentation to illustrate. The first example shows how to write a couple of tests for a hypothetical function called quadraticSolver.

classdef SolverTest < matlab.unittest.TestCase
    % SolverTest tests solutions to the quadratic equation
    % a*x^2 + b*x + c = 0

    methods (Test)
        function testRealSolution(testCase)
            actSolution = quadraticSolver(1,-3,2);
            expSolution = [2,1];
            testCase.verifyEqual(actSolution,expSolution);
        end
        function testImaginarySolution(testCase)
            actSolution = quadraticSolver(1,2,10);
            expSolution = [-1+3i, -1-3i];
            testCase.verifyEqual(actSolution,expSolution);
        end
    end

end

If you need to define setup and teardown methods for your test cases, here's how you do it. (I've just shown the relevant portion of the entire test file.)

properties
    TestFigure
end

methods(TestMethodSetup)
    function createFigure(testCase)
        %comment
        testCase.TestFigure = figure;
    end
end

methods(TestMethodTeardown)
    function closeFigure(testCase)
        close(testCase.TestFigure);
    end
end

Someone asked me just this week about enhancing MATLAB xUnit to provide for setup and teardown methods that get executed just once for all the methods in a test file, instead of being executed once for every individual test method. MATLAB xUnit does not have this capability, but the new unit test framework in R2013a does. A class-level setup method is identified by using the TestClassSetup method attribute, like this:

methods (TestClassSetup)
    function doSomethingBeforeAllTestMethodsAreExecuted(testCase)
        disp('Hi, Mom!')
    end
end

There are lots of new capabilities in this framework. A little birdie told me that more information about matlab.unittest may appear pretty soon in Loren's blog. For now, though, I'd like to point you to the overview video and to the documentation.

That brings me back to MATLAB xUnit and its fate. Well, this little package is entering its "retirement phase." I imagine it will remain on the File Exchange for quite a while because it can take a long time for the user community to adopt a new release, and also because there's a lot of customer-written code out there that uses it. However, I do not plan to work on it anymore. Our testing tools team is an excellent group of engineers, and they can move matlab.unittest forward better and faster than I could ever do with my little side project.

So, if you're interested in unit testing and you get your hands on MATLAB R2013a, please give it a try.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Perhaps you've read about Chebfun before here or on Cleve''s blog. Chebfun allows you to do numerical computing with functions and not simply numbers. If not, you might want to investigate. If so, you may be very interested in new capabilities available.

Contents

Major New Capability

Just recently, the Chebfun team released a major new capability on their site. While Chebfun is extremely useful, it was restricted to solving problems in one dimension. Now you may solve problems with two independent variables using Chebfun2! As before, their website includes code you can freely download and many interesting examples. The examples are well documented and easy to learn from. I will soon have time to explore Chebfun2 because of an irrational amount of airline time.

I also love how easy it is to install - you can do it from MATLAB. They kindly post the code right on the website so you can paste it into your running session.

Many Examples

They have many examples included in the software, and posted on their web site, using the publish command in MATLAB. The high-level examples are organized into these categories:

I hope with this breadth of examples, you are able to find ones interesting or relevant (or both) to you. Once you've gotten a chance to try out Chebfun2, I'd love to hear about your experiences. Post them here to share with us.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Dancing in the Spotlight

by Joachim Schlosser

I love trade shows and conferences, because you get the opportunity to speak with people creating very different things all on one day. The Embedded World trade show, held each year end of February in Nuremberg, Germany, is such an event. Every time it is amazing to see the breadth of topics that people approach us on, and it is fun to discuss and find out which of our tools around MATLAB & Simulink can make them successful.

Dancing: Enter Anna.

People might stop because they see our hostess Anna making the cute little Nao robot dance, but they come into discussion because they get the vision about what they could achieve themselves. What they see on the booth is just a catalyst for people’s own insight and needs, just a starting point to evaluate which tools serve best their applications.

Photo courtesy Joachim Schlosser.

Movement: Enter John.

Enter booth visitor John (not his real name) from Ynapmoc, Inc. (not his real company), a small to mid size manufacturer with 40 years of experience in the area of medical devices (not his real industry). When he approached me, he said he stopped because he noticed we somehow analyze the real picture data, and he was curious to find out more. He then found yet another demonstration in the booth where a Simulink model made an Arduino board with a webcam track and follow a colored ball. Think of it as a stripped down version of this one, but you get the idea.

Now the reason John was curious is because they are abstracting movements of people filmed in their lab to detect anomalies in their muscle tissue functions and provide a self-assessment tool (not his real application). They had been coding it by hand, with little possibility to do fast prototyping, and with additional effort to adapt the functionality each time they got another camera, each time they have another idea to make the algorithm better. To see that he could do video processing in Simulink was a relief to him, since he knew it from his engineering studies. To see additionally that he could real-time prototype on a low cost Arduino board made John enthusiastic.

Grid stability: Enter Lynn.

Enter Lynn (not her real name), engineer from Sunshine Corp (not her real company), who works on electrical energy generators used in production plants as backup to the regular power grid (not her real application). She had seen our offerings around simulating electrical systems and wondered whether this would not just work for simulating an isolated device but the production plant’s power grid as a whole.

Well, it’s not just possible to simulate an entire plant’s grid, Transpower from New Zealand went a bit further and simulated the entire national grid. They do this to calculate the reserve they need to provide in order to maintain grid stability. So, Lynn has made her next step in solving her challenge, discussing the scenario with one of our Application Engineers to find out the exact needs and a plan for implementing the solution.

Now what if you combine both, image processing and control?

3D Imaging: Enter Nao.

Enter Nao (its real name) from Aldebaran (its real company). Nao is the little robot from the beginning of this story. Think of movement. Think of control theory. Think of feedback control. Think of Simulink. And now attach Microsoft® Kinect® with its support in MATLAB & Simulink, one of the many hardware resources for project-based learning. Kinect provides a 3D image of a person, which is transformed into a mesh in MATLAB & Simulink, allowing it to calculate trajectories for the robot. The rest is normal control system theory with Model-Based Design.

You get: A dancing machine being able to mimic movements of a human. A lot of fun in the realms of engineering, science, commercial and education developments, controls, video, data analysis. And lots of nice discussions with bright people of different disciplines who either are already using MATLAB & Simulink or are exploring their way into doing so. To live and breathe their art.

2012 was an eventful year for MATLAB Mobile, with releases for the iPad, the Android platform, and several enhancements to graphics and usability.

To kick off 2013 on a high note, we are making connecting to the cloud even better.

Introducing… MATLAB Cloud Storage

You can now upload MATLAB files and data to MathWorks Cloud and run them from MATLAB Mobile. You can also download your files from the cloud and save them to your computer.

To upload your files, navigate to the MATLAB Cloud Storage page at https://www.mathworks.com/mobile/cloud-storage.html. To access this page, you will need to have a valid MathWorks Account associated with a license that is current on maintenance. You will also need to have created a user ID on MATLAB Mobile.

This is the second part of my plan to rewrite an Image Processing Toolbox example from 20 years ago using more modern MATLAB language features.

Recall the old app I'm trying to reinvent:

The idea is to experiment with the basic principles of DCT-based image compression. I've decided that I don't want to copy this original layout and functionality exactly. For one thing, computers (and MATLAB!) are a lot faster than they were 20 years ago, and I think we really don't need a separate "Apply" button.

Here's my first attempt at sketching what I'd like to make this time around.

The truth, is, I'm making the code and these posts up as I go along. So far I've only had a chance to get some of the initial pieces up and working. But it's enough to start exploring some of the coding ideas.

Here's what happens when I call the app class, which I've called (for now) dctCompressionExample_v2.

app = dctCompressionExample_v2

app = 

  dctCompressionExample_v2 handle

  Properties:
         OriginalImage: [240x240 double]
    ReconstructedImage: []
            ErrorImage: []
    DCTCoefficientMask: []
    NumDCTCoefficients: 10

You can see I got some output in the Command Window, and I also got a new figure on the screen.

Here's the code so far.

classdef dctCompressionExample_v2 < handle

    properties (SetAccess = private)
        OriginalImage
        ReconstructedImage
        ErrorImage
        DCTCoefficientMask
    end

    properties
        NumDCTCoefficients = 10
    end

    properties (Access = private)
        Slider
        OriginalImageAxes
        ReconstructedImageAxes
        ErrorImageAxes
        MaskAxes
    end

    methods
        function this_app = dctCompressionExample_v2
            this_app.OriginalImage = initialImage;
            layoutApp(this_app);
            update(this_app);
        end

        function layoutApp(app)
            f = figure;
            app.OriginalImageAxes = axes('Parent',f, ...
                'Position',[0.12 0.56 0.28 0.28]);
            app.ReconstructedImageAxes = axes('Parent',f, ...
                'Position',[0.60 0.56 0.28 0.28]);
            app.MaskAxes = axes('Parent',f, ...
                'Position',[0.12 0.16 0.28 0.28]);
            app.ErrorImageAxes = axes('Parent',f, ...
                'Position',[0.60 0.16 0.28 0.28]);
            app.Slider = uicontrol('Style','slider', ...
                'Value',app.NumDCTCoefficients/64, ...
                'Units','normalized', ...
                'Position',[0.12 0.08 0.28 0.04], ...
                'Callback',@(~,~,~) app.reactToSliderChange);

        end

        function reactToSliderChange(app)
            v = get(app.Slider,'Value');
            app.NumDCTCoefficients = round(64*v);
            update(app);
        end

        function update(app)
            imshow(app.OriginalImage,'Parent',app.OriginalImageAxes);
            title(app.OriginalImageAxes,'Original Image');

            imshow(app.OriginalImage,'Parent',app.ReconstructedImageAxes);
            title(app.ReconstructedImageAxes,'Reconstructed Image');

            imshow(zeros(size(app.OriginalImage)),'Parent',app.ErrorImageAxes);
            title(app.ErrorImageAxes,'Error Image');

            imshow(zeros(size(app.OriginalImage)),'Parent',app.MaskAxes);
            title(app.MaskAxes,'DCT Coefficient Mask');
            drawnow;
        end
    end

end

function I = initialImage
pout = imread('pout.tif');
pout2 = pout(1:240,:);
I = im2double(adapthisteq(pout2));
end

At this point I'll remind you of Dave Garrison's article that I mentioned last time. This article teaches about how this kind of code works. I don't intend to repeat all of that article here. I'll just point out a few things.

Here's the function that runs when you launch the app.

        function this_app = dctCompressionExample_v2
            this_app.OriginalImage = initialImage;
            layoutApp(this_app);
            update(this_app);
        end

initialImage is a little function I stuck at the bottom of the file. It's purpose is to create our sample image, which for the purpose of this simple DCT demonstration needs to be a multiple of 8-by-8 in size. (Plus, I threw in a call to adapthisteq because the original image has low contrast.)

function I = initialImage
pout = imread('pout.tif');
pout2 = pout(1:240,:);
I = im2double(adapthisteq(pout2));
end

Then there's the properties block.

    properties (SetAccess = private)
        OriginalImage
        ReconstructedImage
        ErrorImage
        DCTCoefficientMask
    end

    properties
        NumDCTCoefficients = 10
    end

    properties (Access = private)
        Slider
        OriginalImageAxes
        ReconstructedImageAxes
        ErrorImageAxes
        MaskAxes
    end

Some of the properties, such as OriginalImage, I want to make accessible but read-only. Some of them, such as OriginalImageAxes, are private because they are only needed by the guts of the app. Notice that one of them, NumDCTCoefficients, is both readable and settable. I wanted to explore the idea of giving this app a little programmable interface so that you can write code to make it change.

I would like to make it so that the NumDCTCoefficients property changes when you drag the slider. Here's how to accomplish that:

            app.Slider = uicontrol('Style','slider', ...
                'Value',app.NumDCTCoefficients/64, ...
                'Units','normalized', ...
                'Position',[0.12 0.08 0.28 0.04], ...
                'Callback',@(~,~,~) app.reactToSliderChange);

Setting the 'Callback' property of the slider like this causes the app's function reactToSliderChange to get called whenever the slider is dragged. Here's what reactToSliderChange does:

        function reactToSliderChange(app)
            v = get(app.Slider,'Value');
            app.NumDCTCoefficients = round(64*v);
            update(app);
        end

Now we're starting to get a little interactivity going in our app.

That's about as far as I got this week. Next time we'll continue wiring things up and maybe start putting in some algorithmic DCT block processing code.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Go play around with it and tell me what you find. Shown above is one view: author ID vs. problem ID. It shows patterns of how problem authors contribute. The size of the dots is relative to the number of comments on the problem. That horizontal stripe across this middle at around author_id = 3000 is super contributor (and current Cody front runner) Richard Zapor, who has single-handedly contributed 122 problems!

This tool loads a relatively large (200 KB) data set when it starts up, so please be patient as it loads! Also, please note that this is an experimental tool that works with an archived static data set from February 8, 2013. It won’t reflect changes in Cody after that date.

Last month, Dave Garrison wrote a nice article on the MathWorks web site about coding an app (or GUI) by defining a class.

Curiously, this article made me think of a day back in September 1993. I arrived at MathWorks on that day to interview for a job as an image processing software developer. Sometime in the afternoon I found myself in an office with the two developers of version 1.0 of the Image Processing Toolbox. (One of those developers, Loren Shure, was the first employee hired by the MathWorks founders. You can find her over at Loren on the Art of MATLAB.)

Anyway, the Image Processing Toolbox was about to ship, and the developers were trying to think of some demo ideas. I joined in a brainstorming session in front of a whiteboard. One of my ideas that day became the product demo called dctdemo. Here is a screen shot.

(That's a picture of Loren, by the way.)

Prompted by Dave's article, I thought it might be a good time to revisit this oldie-but-goodie demo. It was originally written using MATLAB 4, and the programming techniques available at that time for doing this sort of thing were relatively crude. My plan is to take a fresh look at this demo and reimplement it using the techniques described by Dave.

Here's a start:

classdef dctCompressionDemo < handle

end

That's enough code to actually execute.

dctCompressionDemo

ans = 

  dctCompressionDemo handle with no properties.

And that's about all it does (for now). Next time we'll start filling in the details.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

A couple of weeks ago, Cleve wrote a post about George Forsythe, Cleve's "thesis advisor, colleague, and friend."

In that post, Cleve showed some MATLAB code to load in and display a 1964 picture of the organizing committee of the Gatlinburg conferences on numerical algebra.

load gatlin
image(X)
colormap(map)
axis image
axis off

As Cleve mentioned, this picture is one of the very first images distributed with MATLAB. (I recall getting the first MATLAB version capable of image display sometime around 1990. It was very exciting!)

As it turns out, there is another interesting method in MATLAB to display this picture (or at least a cropped, low-resolution version of it):

default_image = pow2(get(0,'DefaultImageCData'),47);
numbits = 12 - 9 + 1;
b = bitshift(default_image,-9);
b = fix(b);
b = bitand(b,bitcmp(0,4));
b = b/max(b(:));

imshow(b,'InitialMagnification',200)

Forsythe is on the right of the cropped version.

So, what's going on here? Well, it turns out that, if you call image with no input arguments (no data), an image with "default" pixel values gets displayed. Furthermore, that image has many images "hidden" within it. For the full story, see "The Story Behind the Default MATLAB Image."

Cleve knew that I was working on this little hidden-image project in the mid-1990s, and I asked him for suggestions about images to include. He suggested this Gatlinburg Conference picture, as well as the 4-by-4 Durer magic square.

c = bitshift(default_image,-23);
c = fix(c);
c = bitand(c,bitcmp(0,5));
c = c/max(c(:));

imshow(c,'InitialMagnification',200)

Thanks for the suggestions, Cleve!

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

It may feel like cheating at first, but the File Exchange has a mechanism to let you acknowledge the person (or people) that you borrowed code from.

When you submit a file to the File Exchange, you come across a section called “License and Acknowledgments.” This gives you a place to mention the other files that helped you create yours, either because you borrowed code from them or because you just want to publicly acknowledge your appreciation of their influence.

Here’s an example. When Douglas Schwarz wrote his highly regarded Fast and Robust Curve Intersections, he made sure to acknowledge the pre-existing file Curve Intersect 2 by Sebastian Hölz. We’ve had this feature in place for many years, and I decided it was time to take a look at the patterns of attribution that we see.

Of the more than 17,000 files on the File Exchange, there are 2910 separate acknowledgments like the one mentioned above. In one way or another, 2050 files acknowledge their debt to a total of 1870 inspiring files. If you think of each one of these acknowledgments as a pointer from the inspiring file to the inspired file, then the whole site can be thought of as a directed graph. Consider the example above. The inspiration doesn’t stop with Douglas Schwarz’s file. That one went on to inspire others. Here’s a look at the acknowledgment tree.

In this picture, each box is a different submission labeled with the author’s last name. Isn’t that lovely? It’s exciting to see the “footprints” of good ideas as they move through the site.

These patterns of acknowledgment can be extensive and visually striking. Some of the trees are deep and some are wide. Here’s a wide one that starts with Oliver Woodford’s exportfig.

With our link data, we can figure out which is the single most influential file in terms of link count. And the answer is… one of the oldest files on the site from Christophe Couvreur. Take a look at how wide this one is.

The largest fully-connected subnet of files has over 200 files in it, all of them connected to each other through a series of links, each link a thoughtful token of respect.

If we go from a file-centric to an author-centric view, we can determine the most influential authors. Who has the most acknowledgments considering all the files they submitted (we ignore self-linking for this analysis). In other words, which authors are the most influential?

1. John D’Errico
2. Christophe COUVREUR
3. Yair Altman
4. Jos
5. Scott Hirsch
6. Yi Cao
7. Malcolm Lidierth
8. Robert Bemis
9. Douglas Schwarz
10. Dirk-Jan Kroon

These people are heroes of the MATLAB Community. They have freely given, and we have all benefited.

I want to see your name on this list!

Let’s say you agree. Let’s say you’re ready to let people start building on your code. How do you start? Paradoxically, a good way to start is by borrowing someone else’s code. You’ll learn good coding practices by borrowing from the best, and you’re likelier to build something worth borrowing by standing on the shoulders of others. You’ll also be bound more deeply into the community. And soon you’ll see how satisfying it is to have someone use your code as a springboard for something wonderful that never would have occurred to you. So give by taking. You’ll be surprised how far it gets you.